1 /********************************************************************** 2 American Put Options Pricing using the Black-Scholes Equation 3 4 Background (European Options): 5 The standard European option is a contract where the holder has the right 6 to either buy (call option) or sell (put option) an underlying asset at 7 a designated future time and price. 8 9 The classic Black-Scholes model begins with an assumption that the 10 price of the underlying asset behaves as a lognormal random walk. 11 Using this assumption and a no-arbitrage argument, the following 12 linear parabolic partial differential equation for the value of the 13 option results: 14 15 dV/dt + 0.5(sigma**2)(S**alpha)(d2V/dS2) + (r - D)S(dV/dS) - rV = 0. 16 17 Here, sigma is the volatility of the underling asset, alpha is a 18 measure of elasticity (typically two), D measures the dividend payments 19 on the underling asset, and r is the interest rate. 20 21 To completely specify the problem, we need to impose some boundary 22 conditions. These are as follows: 23 24 V(S, T) = max(E - S, 0) 25 V(0, t) = E for all 0 <= t <= T 26 V(s, t) = 0 for all 0 <= t <= T and s->infinity 27 28 where T is the exercise time time and E the strike price (price paid 29 for the contract). 30 31 An explicit formula for the value of an European option can be 32 found. See the references for examples. 33 34 Background (American Options): 35 The American option is similar to its European counterpart. The 36 difference is that the holder of the American option can exercise 37 their right to buy or sell the asset at any time prior to the 38 expiration. This additional ability introduce a free boundary into 39 the Black-Scholes equation which can be modeled as a linear 40 complementarity problem. 41 42 0 <= -(dV/dt + 0.5(sigma**2)(S**alpha)(d2V/dS2) + (r - D)S(dV/dS) - rV) 43 complements 44 V(S,T) >= max(E-S,0) 45 46 where the variables are the same as before and we have the same boundary 47 conditions. 48 49 There is not explicit formula for calculating the value of an American 50 option. Therefore, we discretize the above problem and solve the 51 resulting linear complementarity problem. 52 53 We will use backward differences for the time variables and central 54 differences for the space variables. Crank-Nicholson averaging will 55 also be used in the discretization. The algorithm used by the code 56 solves for V(S,t) for a fixed t and then uses this value in the 57 calculation of V(S,t - dt). The method stops when V(S,0) has been 58 found. 59 60 References: 61 + * - Huang and Pang, "Options Pricing and Linear Complementarity," 62 Journal of Computational Finance, volume 2, number 3, 1998. 63 - * - Wilmott, "Derivatives: The Theory and Practice of Financial Engineering," 64 John Wiley and Sons, New York, 1998. 65 ***************************************************************************/ 66 67 /* 68 Include "petsctao.h" so we can use TAO solvers. 69 Include "petscdmda.h" so that we can use distributed meshes (DMs) for managing 70 the parallel mesh. 71 */ 72 73 #include <petscdmda.h> 74 #include <petsctao.h> 75 76 static char help[] = 77 "This example demonstrates use of the TAO package to\n\ 78 solve a linear complementarity problem for pricing American put options.\n\ 79 The code uses backward differences in time and central differences in\n\ 80 space. The command line options are:\n\ 81 -rate <r>, where <r> = interest rate\n\ 82 -sigma <s>, where <s> = volatility of the underlying\n\ 83 -alpha <a>, where <a> = elasticity of the underlying\n\ 84 -delta <d>, where <d> = dividend rate\n\ 85 -strike <e>, where <e> = strike price\n\ 86 -expiry <t>, where <t> = the expiration date\n\ 87 -mt <tg>, where <tg> = number of grid points in time\n\ 88 -ms <sg>, where <sg> = number of grid points in space\n\ 89 -es <se>, where <se> = ending point of the space discretization\n\n"; 90 91 /*T 92 Concepts: TAO^Solving a complementarity problem 93 Routines: TaoCreate(); TaoDestroy(); 94 Routines: TaoSetJacobianRoutine(); TaoAppSetConstraintRoutine(); 95 Routines: TaoSetFromOptions(); 96 Routines: TaoSolveApplication(); 97 Routines: TaoSetVariableBoundsRoutine(); TaoSetInitialSolutionVec(); 98 Processors: 1 99 T*/ 100 101 /* 102 User-defined application context - contains data needed by the 103 application-provided call-back routines, FormFunction(), and FormJacobian(). 104 */ 105 106 typedef struct { 107 PetscReal *Vt1; /* Value of the option at time T + dt */ 108 PetscReal *c; /* Constant -- (r - D)S */ 109 PetscReal *d; /* Constant -- -0.5(sigma**2)(S**alpha) */ 110 111 PetscReal rate; /* Interest rate */ 112 PetscReal sigma, alpha, delta; /* Underlying asset properties */ 113 PetscReal strike, expiry; /* Option contract properties */ 114 115 PetscReal es; /* Finite value used for maximum asset value */ 116 PetscReal ds, dt; /* Discretization properties */ 117 PetscInt ms, mt; /* Number of elements */ 118 119 DM dm; 120 } AppCtx; 121 122 /* -------- User-defined Routines --------- */ 123 124 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 125 PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *); 126 PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void*); 127 128 int main(int argc, char **argv) 129 { 130 Vec x; /* solution vector */ 131 Vec c; /* Constraints function vector */ 132 Mat J; /* jacobian matrix */ 133 PetscBool flg; /* A return variable when checking for user options */ 134 Tao tao; /* Tao solver context */ 135 AppCtx user; /* user-defined work context */ 136 PetscInt i, j; 137 PetscInt xs,xm,gxs,gxm; 138 PetscReal sval = 0; 139 PetscReal *x_array; 140 Vec localX; 141 142 /* Initialize PETSc, TAO */ 143 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 144 145 /* 146 Initialize the user-defined application context with reasonable 147 values for the American option to price 148 */ 149 user.rate = 0.04; 150 user.sigma = 0.40; 151 user.alpha = 2.00; 152 user.delta = 0.01; 153 user.strike = 10.0; 154 user.expiry = 1.0; 155 user.mt = 10; 156 user.ms = 150; 157 user.es = 100.0; 158 159 /* Read in alternative values for the American option to price */ 160 PetscCall(PetscOptionsGetReal(NULL,NULL, "-alpha", &user.alpha, &flg)); 161 PetscCall(PetscOptionsGetReal(NULL,NULL, "-delta", &user.delta, &flg)); 162 PetscCall(PetscOptionsGetReal(NULL,NULL, "-es", &user.es, &flg)); 163 PetscCall(PetscOptionsGetReal(NULL,NULL, "-expiry", &user.expiry, &flg)); 164 PetscCall(PetscOptionsGetInt(NULL,NULL, "-ms", &user.ms, &flg)); 165 PetscCall(PetscOptionsGetInt(NULL,NULL, "-mt", &user.mt, &flg)); 166 PetscCall(PetscOptionsGetReal(NULL,NULL, "-rate", &user.rate, &flg)); 167 PetscCall(PetscOptionsGetReal(NULL,NULL, "-sigma", &user.sigma, &flg)); 168 PetscCall(PetscOptionsGetReal(NULL,NULL, "-strike", &user.strike, &flg)); 169 170 /* Check that the options set are allowable (needs to be done) */ 171 172 user.ms++; 173 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.ms,1,1,NULL,&user.dm)); 174 PetscCall(DMSetFromOptions(user.dm)); 175 PetscCall(DMSetUp(user.dm)); 176 /* Create appropriate vectors and matrices */ 177 178 PetscCall(DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL)); 179 PetscCall(DMDAGetGhostCorners(user.dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 180 181 PetscCall(DMCreateGlobalVector(user.dm,&x)); 182 /* 183 Finish filling in the user-defined context with the values for 184 dS, dt, and allocating space for the constants 185 */ 186 user.ds = user.es / (user.ms-1); 187 user.dt = user.expiry / user.mt; 188 189 PetscCall(PetscMalloc1(gxm,&(user.Vt1))); 190 PetscCall(PetscMalloc1(gxm,&(user.c))); 191 PetscCall(PetscMalloc1(gxm,&(user.d))); 192 193 /* 194 Calculate the values for the constant. Vt1 begins with the ending 195 boundary condition. 196 */ 197 for (i=0; i<gxm; i++) { 198 sval = (gxs+i)*user.ds; 199 user.Vt1[i] = PetscMax(user.strike - sval, 0); 200 user.c[i] = (user.delta - user.rate)*sval; 201 user.d[i] = -0.5*user.sigma*user.sigma*PetscPowReal(sval, user.alpha); 202 } 203 if (gxs+gxm==user.ms) { 204 user.Vt1[gxm-1] = 0; 205 } 206 PetscCall(VecDuplicate(x, &c)); 207 208 /* 209 Allocate the matrix used by TAO for the Jacobian. Each row of 210 the Jacobian matrix will have at most three elements. 211 */ 212 PetscCall(DMCreateMatrix(user.dm,&J)); 213 214 /* The TAO code begins here */ 215 216 /* Create TAO solver and set desired solution method */ 217 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 218 PetscCall(TaoSetType(tao,TAOSSILS)); 219 220 /* Set routines for constraints function and Jacobian evaluation */ 221 PetscCall(TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user)); 222 PetscCall(TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user)); 223 224 /* Set the variable bounds */ 225 PetscCall(TaoSetVariableBoundsRoutine(tao,ComputeVariableBounds,(void*)&user)); 226 227 /* Set initial solution guess */ 228 PetscCall(VecGetArray(x,&x_array)); 229 for (i=0; i< xm; i++) 230 x_array[i] = user.Vt1[i-gxs+xs]; 231 PetscCall(VecRestoreArray(x,&x_array)); 232 /* Set data structure */ 233 PetscCall(TaoSetSolution(tao, x)); 234 235 /* Set routines for function and Jacobian evaluation */ 236 PetscCall(TaoSetFromOptions(tao)); 237 238 /* Iteratively solve the linear complementarity problems */ 239 for (i = 1; i < user.mt; i++) { 240 241 /* Solve the current version */ 242 PetscCall(TaoSolve(tao)); 243 244 /* Update Vt1 with the solution */ 245 PetscCall(DMGetLocalVector(user.dm,&localX)); 246 PetscCall(DMGlobalToLocalBegin(user.dm,x,INSERT_VALUES,localX)); 247 PetscCall(DMGlobalToLocalEnd(user.dm,x,INSERT_VALUES,localX)); 248 PetscCall(VecGetArray(localX,&x_array)); 249 for (j = 0; j < gxm; j++) { 250 user.Vt1[j] = x_array[j]; 251 } 252 PetscCall(VecRestoreArray(x,&x_array)); 253 PetscCall(DMRestoreLocalVector(user.dm,&localX)); 254 } 255 256 /* Free TAO data structures */ 257 PetscCall(TaoDestroy(&tao)); 258 259 /* Free PETSc data structures */ 260 PetscCall(VecDestroy(&x)); 261 PetscCall(VecDestroy(&c)); 262 PetscCall(MatDestroy(&J)); 263 PetscCall(DMDestroy(&user.dm)); 264 /* Free user-defined workspace */ 265 PetscCall(PetscFree(user.Vt1)); 266 PetscCall(PetscFree(user.c)); 267 PetscCall(PetscFree(user.d)); 268 269 PetscCall(PetscFinalize()); 270 return 0; 271 } 272 273 /* -------------------------------------------------------------------- */ 274 PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, void*ctx) 275 { 276 AppCtx *user = (AppCtx *) ctx; 277 PetscInt i; 278 PetscInt xs,xm; 279 PetscInt ms = user->ms; 280 PetscReal sval=0.0,*xl_array,ub= PETSC_INFINITY; 281 282 /* Set the variable bounds */ 283 PetscCall(VecSet(xu, ub)); 284 PetscCall(DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL)); 285 286 PetscCall(VecGetArray(xl,&xl_array)); 287 for (i = 0; i < xm; i++) { 288 sval = (xs+i)*user->ds; 289 xl_array[i] = PetscMax(user->strike - sval, 0); 290 } 291 PetscCall(VecRestoreArray(xl,&xl_array)); 292 293 if (xs==0) { 294 PetscCall(VecGetArray(xu,&xl_array)); 295 xl_array[0] = PetscMax(user->strike, 0); 296 PetscCall(VecRestoreArray(xu,&xl_array)); 297 } 298 if (xs+xm==ms) { 299 PetscCall(VecGetArray(xu,&xl_array)); 300 xl_array[xm-1] = 0; 301 PetscCall(VecRestoreArray(xu,&xl_array)); 302 } 303 304 return 0; 305 } 306 /* -------------------------------------------------------------------- */ 307 308 /* 309 FormFunction - Evaluates gradient of f. 310 311 Input Parameters: 312 . tao - the Tao context 313 . X - input vector 314 . ptr - optional user-defined context, as set by TaoAppSetConstraintRoutine() 315 316 Output Parameters: 317 . F - vector containing the newly evaluated gradient 318 */ 319 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec F, void *ptr) 320 { 321 AppCtx *user = (AppCtx *) ptr; 322 PetscReal *x, *f; 323 PetscReal *Vt1 = user->Vt1, *c = user->c, *d = user->d; 324 PetscReal rate = user->rate; 325 PetscReal dt = user->dt, ds = user->ds; 326 PetscInt ms = user->ms; 327 PetscInt i, xs,xm,gxs,gxm; 328 Vec localX,localF; 329 PetscReal zero=0.0; 330 331 PetscCall(DMGetLocalVector(user->dm,&localX)); 332 PetscCall(DMGetLocalVector(user->dm,&localF)); 333 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 334 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 335 PetscCall(DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL)); 336 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 337 PetscCall(VecSet(F, zero)); 338 /* 339 The problem size is smaller than the discretization because of the 340 two fixed elements (V(0,T) = E and V(Send,T) = 0. 341 */ 342 343 /* Get pointers to the vector data */ 344 PetscCall(VecGetArray(localX, &x)); 345 PetscCall(VecGetArray(localF, &f)); 346 347 /* Left Boundary */ 348 if (gxs==0) { 349 f[0] = x[0]-user->strike; 350 } else { 351 f[0] = 0; 352 } 353 354 /* All points in the interior */ 355 /* for (i=gxs+1;i<gxm-1;i++) { */ 356 for (i=1;i<gxm-1;i++) { 357 f[i] = (1.0/dt + rate)*x[i] - Vt1[i]/dt + (c[i]/(4*ds))*(x[i+1] - x[i-1] + Vt1[i+1] - Vt1[i-1]) + 358 (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] + Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]); 359 } 360 361 /* Right boundary */ 362 if (gxs+gxm==ms) { 363 f[xm-1]=x[gxm-1]; 364 } else { 365 f[xm-1]=0; 366 } 367 368 /* Restore vectors */ 369 PetscCall(VecRestoreArray(localX, &x)); 370 PetscCall(VecRestoreArray(localF, &f)); 371 PetscCall(DMLocalToGlobalBegin(user->dm,localF,INSERT_VALUES,F)); 372 PetscCall(DMLocalToGlobalEnd(user->dm,localF,INSERT_VALUES,F)); 373 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 374 PetscCall(DMRestoreLocalVector(user->dm,&localF)); 375 PetscCall(PetscLogFlops(24.0*(gxm-2))); 376 /* 377 info=VecView(F,PETSC_VIEWER_STDOUT_WORLD); 378 */ 379 return 0; 380 } 381 382 /* ------------------------------------------------------------------- */ 383 /* 384 FormJacobian - Evaluates Jacobian matrix. 385 386 Input Parameters: 387 . tao - the Tao context 388 . X - input vector 389 . ptr - optional user-defined context, as set by TaoSetJacobian() 390 391 Output Parameters: 392 . J - Jacobian matrix 393 */ 394 PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr) 395 { 396 AppCtx *user = (AppCtx *) ptr; 397 PetscReal *c = user->c, *d = user->d; 398 PetscReal rate = user->rate; 399 PetscReal dt = user->dt, ds = user->ds; 400 PetscInt ms = user->ms; 401 PetscReal val[3]; 402 PetscInt col[3]; 403 PetscInt i; 404 PetscInt gxs,gxm; 405 PetscBool assembled; 406 407 /* Set various matrix options */ 408 PetscCall(MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE)); 409 PetscCall(MatAssembled(J,&assembled)); 410 if (assembled) PetscCall(MatZeroEntries(J)); 411 412 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 413 414 if (gxs==0) { 415 i = 0; 416 col[0] = 0; 417 val[0]=1.0; 418 PetscCall(MatSetValues(J,1,&i,1,col,val,INSERT_VALUES)); 419 } 420 for (i=1; i < gxm-1; i++) { 421 col[0] = gxs + i - 1; 422 col[1] = gxs + i; 423 col[2] = gxs + i + 1; 424 val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds); 425 val[1] = 1.0/dt + rate - d[i]/(ds*ds); 426 val[2] = c[i]/(4*ds) + d[i]/(2*ds*ds); 427 PetscCall(MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES)); 428 } 429 if (gxs+gxm==ms) { 430 i = ms-1; 431 col[0] = i; 432 val[0]=1.0; 433 PetscCall(MatSetValues(J,1,&i,1,col,val,INSERT_VALUES)); 434 } 435 436 /* Assemble the Jacobian matrix */ 437 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 438 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 439 PetscCall(PetscLogFlops(18.0*(gxm)+5)); 440 return 0; 441 } 442 443 /*TEST 444 445 build: 446 requires: !complex 447 448 test: 449 args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5 450 requires: !single 451 452 test: 453 suffix: 2 454 args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5 455 requires: !single 456 457 test: 458 suffix: 3 459 args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5 460 requires: !single 461 462 test: 463 suffix: 4 464 args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5 465 requires: !single 466 467 test: 468 suffix: 5 469 args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5 470 requires: !single 471 472 test: 473 suffix: 6 474 args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5 475 requires: !single 476 477 test: 478 suffix: 7 479 args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5 480 requires: !single 481 482 TEST*/ 483