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