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