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