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 /* 92 User-defined application context - contains data needed by the 93 application-provided call-back routines, FormFunction(), and FormJacobian(). 94 */ 95 96 typedef struct { 97 PetscReal *Vt1; /* Value of the option at time T + dt */ 98 PetscReal *c; /* Constant -- (r - D)S */ 99 PetscReal *d; /* Constant -- -0.5(sigma**2)(S**alpha) */ 100 101 PetscReal rate; /* Interest rate */ 102 PetscReal sigma, alpha, delta; /* Underlying asset properties */ 103 PetscReal strike, expiry; /* Option contract properties */ 104 105 PetscReal es; /* Finite value used for maximum asset value */ 106 PetscReal ds, dt; /* Discretization properties */ 107 PetscInt ms, mt; /* Number of elements */ 108 109 DM dm; 110 } AppCtx; 111 112 /* -------- User-defined Routines --------- */ 113 114 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *); 115 PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *); 116 PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void*); 117 118 int main(int argc, char **argv) 119 { 120 Vec x; /* solution vector */ 121 Vec c; /* Constraints function vector */ 122 Mat J; /* jacobian matrix */ 123 PetscBool flg; /* A return variable when checking for user options */ 124 Tao tao; /* Tao solver context */ 125 AppCtx user; /* user-defined work context */ 126 PetscInt i, j; 127 PetscInt xs,xm,gxs,gxm; 128 PetscReal sval = 0; 129 PetscReal *x_array; 130 Vec localX; 131 132 /* Initialize PETSc, TAO */ 133 PetscFunctionBeginUser; 134 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 135 136 /* 137 Initialize the user-defined application context with reasonable 138 values for the American option to price 139 */ 140 user.rate = 0.04; 141 user.sigma = 0.40; 142 user.alpha = 2.00; 143 user.delta = 0.01; 144 user.strike = 10.0; 145 user.expiry = 1.0; 146 user.mt = 10; 147 user.ms = 150; 148 user.es = 100.0; 149 150 /* Read in alternative values for the American option to price */ 151 PetscCall(PetscOptionsGetReal(NULL,NULL, "-alpha", &user.alpha, &flg)); 152 PetscCall(PetscOptionsGetReal(NULL,NULL, "-delta", &user.delta, &flg)); 153 PetscCall(PetscOptionsGetReal(NULL,NULL, "-es", &user.es, &flg)); 154 PetscCall(PetscOptionsGetReal(NULL,NULL, "-expiry", &user.expiry, &flg)); 155 PetscCall(PetscOptionsGetInt(NULL,NULL, "-ms", &user.ms, &flg)); 156 PetscCall(PetscOptionsGetInt(NULL,NULL, "-mt", &user.mt, &flg)); 157 PetscCall(PetscOptionsGetReal(NULL,NULL, "-rate", &user.rate, &flg)); 158 PetscCall(PetscOptionsGetReal(NULL,NULL, "-sigma", &user.sigma, &flg)); 159 PetscCall(PetscOptionsGetReal(NULL,NULL, "-strike", &user.strike, &flg)); 160 161 /* Check that the options set are allowable (needs to be done) */ 162 163 user.ms++; 164 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.ms,1,1,NULL,&user.dm)); 165 PetscCall(DMSetFromOptions(user.dm)); 166 PetscCall(DMSetUp(user.dm)); 167 /* Create appropriate vectors and matrices */ 168 169 PetscCall(DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL)); 170 PetscCall(DMDAGetGhostCorners(user.dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 171 172 PetscCall(DMCreateGlobalVector(user.dm,&x)); 173 /* 174 Finish filling in the user-defined context with the values for 175 dS, dt, and allocating space for the constants 176 */ 177 user.ds = user.es / (user.ms-1); 178 user.dt = user.expiry / user.mt; 179 180 PetscCall(PetscMalloc1(gxm,&(user.Vt1))); 181 PetscCall(PetscMalloc1(gxm,&(user.c))); 182 PetscCall(PetscMalloc1(gxm,&(user.d))); 183 184 /* 185 Calculate the values for the constant. Vt1 begins with the ending 186 boundary condition. 187 */ 188 for (i=0; i<gxm; i++) { 189 sval = (gxs+i)*user.ds; 190 user.Vt1[i] = PetscMax(user.strike - sval, 0); 191 user.c[i] = (user.delta - user.rate)*sval; 192 user.d[i] = -0.5*user.sigma*user.sigma*PetscPowReal(sval, user.alpha); 193 } 194 if (gxs+gxm==user.ms) { 195 user.Vt1[gxm-1] = 0; 196 } 197 PetscCall(VecDuplicate(x, &c)); 198 199 /* 200 Allocate the matrix used by TAO for the Jacobian. Each row of 201 the Jacobian matrix will have at most three elements. 202 */ 203 PetscCall(DMCreateMatrix(user.dm,&J)); 204 205 /* The TAO code begins here */ 206 207 /* Create TAO solver and set desired solution method */ 208 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 209 PetscCall(TaoSetType(tao,TAOSSILS)); 210 211 /* Set routines for constraints function and Jacobian evaluation */ 212 PetscCall(TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user)); 213 PetscCall(TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user)); 214 215 /* Set the variable bounds */ 216 PetscCall(TaoSetVariableBoundsRoutine(tao,ComputeVariableBounds,(void*)&user)); 217 218 /* Set initial solution guess */ 219 PetscCall(VecGetArray(x,&x_array)); 220 for (i=0; i< xm; i++) 221 x_array[i] = user.Vt1[i-gxs+xs]; 222 PetscCall(VecRestoreArray(x,&x_array)); 223 /* Set data structure */ 224 PetscCall(TaoSetSolution(tao, x)); 225 226 /* Set routines for function and Jacobian evaluation */ 227 PetscCall(TaoSetFromOptions(tao)); 228 229 /* Iteratively solve the linear complementarity problems */ 230 for (i = 1; i < user.mt; i++) { 231 232 /* Solve the current version */ 233 PetscCall(TaoSolve(tao)); 234 235 /* Update Vt1 with the solution */ 236 PetscCall(DMGetLocalVector(user.dm,&localX)); 237 PetscCall(DMGlobalToLocalBegin(user.dm,x,INSERT_VALUES,localX)); 238 PetscCall(DMGlobalToLocalEnd(user.dm,x,INSERT_VALUES,localX)); 239 PetscCall(VecGetArray(localX,&x_array)); 240 for (j = 0; j < gxm; j++) { 241 user.Vt1[j] = x_array[j]; 242 } 243 PetscCall(VecRestoreArray(x,&x_array)); 244 PetscCall(DMRestoreLocalVector(user.dm,&localX)); 245 } 246 247 /* Free TAO data structures */ 248 PetscCall(TaoDestroy(&tao)); 249 250 /* Free PETSc data structures */ 251 PetscCall(VecDestroy(&x)); 252 PetscCall(VecDestroy(&c)); 253 PetscCall(MatDestroy(&J)); 254 PetscCall(DMDestroy(&user.dm)); 255 /* Free user-defined workspace */ 256 PetscCall(PetscFree(user.Vt1)); 257 PetscCall(PetscFree(user.c)); 258 PetscCall(PetscFree(user.d)); 259 260 PetscCall(PetscFinalize()); 261 return 0; 262 } 263 264 /* -------------------------------------------------------------------- */ 265 PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, void*ctx) 266 { 267 AppCtx *user = (AppCtx *) ctx; 268 PetscInt i; 269 PetscInt xs,xm; 270 PetscInt ms = user->ms; 271 PetscReal sval=0.0,*xl_array,ub= PETSC_INFINITY; 272 273 /* Set the variable bounds */ 274 PetscCall(VecSet(xu, ub)); 275 PetscCall(DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL)); 276 277 PetscCall(VecGetArray(xl,&xl_array)); 278 for (i = 0; i < xm; i++) { 279 sval = (xs+i)*user->ds; 280 xl_array[i] = PetscMax(user->strike - sval, 0); 281 } 282 PetscCall(VecRestoreArray(xl,&xl_array)); 283 284 if (xs==0) { 285 PetscCall(VecGetArray(xu,&xl_array)); 286 xl_array[0] = PetscMax(user->strike, 0); 287 PetscCall(VecRestoreArray(xu,&xl_array)); 288 } 289 if (xs+xm==ms) { 290 PetscCall(VecGetArray(xu,&xl_array)); 291 xl_array[xm-1] = 0; 292 PetscCall(VecRestoreArray(xu,&xl_array)); 293 } 294 295 return 0; 296 } 297 /* -------------------------------------------------------------------- */ 298 299 /* 300 FormFunction - Evaluates gradient of f. 301 302 Input Parameters: 303 . tao - the Tao context 304 . X - input vector 305 . ptr - optional user-defined context, as set by TaoAppSetConstraintRoutine() 306 307 Output Parameters: 308 . F - vector containing the newly evaluated gradient 309 */ 310 PetscErrorCode FormConstraints(Tao tao, Vec X, Vec F, void *ptr) 311 { 312 AppCtx *user = (AppCtx *) ptr; 313 PetscReal *x, *f; 314 PetscReal *Vt1 = user->Vt1, *c = user->c, *d = user->d; 315 PetscReal rate = user->rate; 316 PetscReal dt = user->dt, ds = user->ds; 317 PetscInt ms = user->ms; 318 PetscInt i, xs,xm,gxs,gxm; 319 Vec localX,localF; 320 PetscReal zero=0.0; 321 322 PetscCall(DMGetLocalVector(user->dm,&localX)); 323 PetscCall(DMGetLocalVector(user->dm,&localF)); 324 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 325 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 326 PetscCall(DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL)); 327 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 328 PetscCall(VecSet(F, zero)); 329 /* 330 The problem size is smaller than the discretization because of the 331 two fixed elements (V(0,T) = E and V(Send,T) = 0. 332 */ 333 334 /* Get pointers to the vector data */ 335 PetscCall(VecGetArray(localX, &x)); 336 PetscCall(VecGetArray(localF, &f)); 337 338 /* Left Boundary */ 339 if (gxs==0) { 340 f[0] = x[0]-user->strike; 341 } else { 342 f[0] = 0; 343 } 344 345 /* All points in the interior */ 346 /* for (i=gxs+1;i<gxm-1;i++) { */ 347 for (i=1;i<gxm-1;i++) { 348 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]) + 349 (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] + Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]); 350 } 351 352 /* Right boundary */ 353 if (gxs+gxm==ms) { 354 f[xm-1]=x[gxm-1]; 355 } else { 356 f[xm-1]=0; 357 } 358 359 /* Restore vectors */ 360 PetscCall(VecRestoreArray(localX, &x)); 361 PetscCall(VecRestoreArray(localF, &f)); 362 PetscCall(DMLocalToGlobalBegin(user->dm,localF,INSERT_VALUES,F)); 363 PetscCall(DMLocalToGlobalEnd(user->dm,localF,INSERT_VALUES,F)); 364 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 365 PetscCall(DMRestoreLocalVector(user->dm,&localF)); 366 PetscCall(PetscLogFlops(24.0*(gxm-2))); 367 /* 368 info=VecView(F,PETSC_VIEWER_STDOUT_WORLD); 369 */ 370 return 0; 371 } 372 373 /* ------------------------------------------------------------------- */ 374 /* 375 FormJacobian - Evaluates Jacobian matrix. 376 377 Input Parameters: 378 . tao - the Tao context 379 . X - input vector 380 . ptr - optional user-defined context, as set by TaoSetJacobian() 381 382 Output Parameters: 383 . J - Jacobian matrix 384 */ 385 PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr) 386 { 387 AppCtx *user = (AppCtx *) ptr; 388 PetscReal *c = user->c, *d = user->d; 389 PetscReal rate = user->rate; 390 PetscReal dt = user->dt, ds = user->ds; 391 PetscInt ms = user->ms; 392 PetscReal val[3]; 393 PetscInt col[3]; 394 PetscInt i; 395 PetscInt gxs,gxm; 396 PetscBool assembled; 397 398 /* Set various matrix options */ 399 PetscCall(MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE)); 400 PetscCall(MatAssembled(J,&assembled)); 401 if (assembled) PetscCall(MatZeroEntries(J)); 402 403 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL)); 404 405 if (gxs==0) { 406 i = 0; 407 col[0] = 0; 408 val[0]=1.0; 409 PetscCall(MatSetValues(J,1,&i,1,col,val,INSERT_VALUES)); 410 } 411 for (i=1; i < gxm-1; i++) { 412 col[0] = gxs + i - 1; 413 col[1] = gxs + i; 414 col[2] = gxs + i + 1; 415 val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds); 416 val[1] = 1.0/dt + rate - d[i]/(ds*ds); 417 val[2] = c[i]/(4*ds) + d[i]/(2*ds*ds); 418 PetscCall(MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES)); 419 } 420 if (gxs+gxm==ms) { 421 i = ms-1; 422 col[0] = i; 423 val[0]=1.0; 424 PetscCall(MatSetValues(J,1,&i,1,col,val,INSERT_VALUES)); 425 } 426 427 /* Assemble the Jacobian matrix */ 428 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 429 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 430 PetscCall(PetscLogFlops(18.0*(gxm)+5)); 431 return 0; 432 } 433 434 /*TEST 435 436 build: 437 requires: !complex 438 439 test: 440 args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5 441 requires: !single 442 443 test: 444 suffix: 2 445 args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5 446 requires: !single 447 448 test: 449 suffix: 3 450 args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5 451 requires: !single 452 453 test: 454 suffix: 4 455 args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5 456 requires: !single 457 458 test: 459 suffix: 5 460 args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5 461 requires: !single 462 463 test: 464 suffix: 6 465 args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5 466 requires: !single 467 468 test: 469 suffix: 7 470 args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5 471 requires: !single 472 473 TEST*/ 474