Actual source code: blackscholes.c

  1: /**********************************************************************
  2:     American Put Options Pricing using the Black-Scholes Equation

  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.

  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:

 15:        dV/dt + 0.5(sigma**2)(S**alpha)(d2V/dS2) + (r - D)S(dV/dS) - rV = 0.

 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.

 21:      To completely specify the problem, we need to impose some boundary
 22:      conditions.  These are as follows:

 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

 28:      where T is the exercise time time and E the strike price (price paid
 29:      for the contract).

 31:      An explicit formula for the value of an European option can be
 32:      found.  See the references for examples.

 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.

 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)

 46:      where the variables are the same as before and we have the same boundary
 47:      conditions.

 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.

 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.

 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: ***************************************************************************/

 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: */

 73: #include <petscdmda.h>
 74: #include <petsctao.h>

 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";

 90: /*
 91:   User-defined application context - contains data needed by the
 92:   application-provided call-back routines, FormFunction(), and FormJacobian().
 93: */

 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) */

100:   PetscReal rate;                /* Interest rate */
101:   PetscReal sigma, alpha, delta; /* Underlying asset properties */
102:   PetscReal strike, expiry;      /* Option contract properties */

104:   PetscReal es;     /* Finite value used for maximum asset value */
105:   PetscReal ds, dt; /* Discretization properties */
106:   PetscInt  ms, mt; /* Number of elements */

108:   DM dm;
109: } AppCtx;

111: /* -------- User-defined Routines --------- */

113: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
114: PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
115: PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void *);

117: int main(int argc, char **argv)
118: {
119:   Vec        x;    /* solution vector */
120:   Vec        c;    /* Constraints function vector */
121:   Mat        J;    /* jacobian matrix */
122:   PetscBool  flg;  /* A return variable when checking for user options */
123:   Tao        tao;  /* Tao solver context */
124:   AppCtx     user; /* user-defined work context */
125:   PetscInt   i, j;
126:   PetscInt   xs, xm, gxs, gxm;
127:   PetscReal  sval = 0;
128:   PetscReal *x_array;
129:   Vec        localX;

131:   /* Initialize PETSc, TAO */
133:   PetscInitialize(&argc, &argv, (char *)0, help);

135:   /*
136:      Initialize the user-defined application context with reasonable
137:      values for the American option to price
138:   */
139:   user.rate   = 0.04;
140:   user.sigma  = 0.40;
141:   user.alpha  = 2.00;
142:   user.delta  = 0.01;
143:   user.strike = 10.0;
144:   user.expiry = 1.0;
145:   user.mt     = 10;
146:   user.ms     = 150;
147:   user.es     = 100.0;

149:   /* Read in alternative values for the American option to price */
150:   PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg);
151:   PetscOptionsGetReal(NULL, NULL, "-delta", &user.delta, &flg);
152:   PetscOptionsGetReal(NULL, NULL, "-es", &user.es, &flg);
153:   PetscOptionsGetReal(NULL, NULL, "-expiry", &user.expiry, &flg);
154:   PetscOptionsGetInt(NULL, NULL, "-ms", &user.ms, &flg);
155:   PetscOptionsGetInt(NULL, NULL, "-mt", &user.mt, &flg);
156:   PetscOptionsGetReal(NULL, NULL, "-rate", &user.rate, &flg);
157:   PetscOptionsGetReal(NULL, NULL, "-sigma", &user.sigma, &flg);
158:   PetscOptionsGetReal(NULL, NULL, "-strike", &user.strike, &flg);

160:   /* Check that the options set are allowable (needs to be done) */

162:   user.ms++;
163:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.ms, 1, 1, NULL, &user.dm);
164:   DMSetFromOptions(user.dm);
165:   DMSetUp(user.dm);
166:   /* Create appropriate vectors and matrices */

168:   DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL);
169:   DMDAGetGhostCorners(user.dm, &gxs, NULL, NULL, &gxm, NULL, NULL);

171:   DMCreateGlobalVector(user.dm, &x);
172:   /*
173:      Finish filling in the user-defined context with the values for
174:      dS, dt, and allocating space for the constants
175:   */
176:   user.ds = user.es / (user.ms - 1);
177:   user.dt = user.expiry / user.mt;

179:   PetscMalloc1(gxm, &(user.Vt1));
180:   PetscMalloc1(gxm, &(user.c));
181:   PetscMalloc1(gxm, &(user.d));

183:   /*
184:      Calculate the values for the constant.  Vt1 begins with the ending
185:      boundary condition.
186:   */
187:   for (i = 0; i < gxm; i++) {
188:     sval        = (gxs + i) * user.ds;
189:     user.Vt1[i] = PetscMax(user.strike - sval, 0);
190:     user.c[i]   = (user.delta - user.rate) * sval;
191:     user.d[i]   = -0.5 * user.sigma * user.sigma * PetscPowReal(sval, user.alpha);
192:   }
193:   if (gxs + gxm == user.ms) user.Vt1[gxm - 1] = 0;
194:   VecDuplicate(x, &c);

196:   /*
197:      Allocate the matrix used by TAO for the Jacobian.  Each row of
198:      the Jacobian matrix will have at most three elements.
199:   */
200:   DMCreateMatrix(user.dm, &J);

202:   /* The TAO code begins here */

204:   /* Create TAO solver and set desired solution method  */
205:   TaoCreate(PETSC_COMM_WORLD, &tao);
206:   TaoSetType(tao, TAOSSILS);

208:   /* Set routines for constraints function and Jacobian evaluation */
209:   TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
210:   TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);

212:   /* Set the variable bounds */
213:   TaoSetVariableBoundsRoutine(tao, ComputeVariableBounds, (void *)&user);

215:   /* Set initial solution guess */
216:   VecGetArray(x, &x_array);
217:   for (i = 0; i < xm; i++) x_array[i] = user.Vt1[i - gxs + xs];
218:   VecRestoreArray(x, &x_array);
219:   /* Set data structure */
220:   TaoSetSolution(tao, x);

222:   /* Set routines for function and Jacobian evaluation */
223:   TaoSetFromOptions(tao);

225:   /* Iteratively solve the linear complementarity problems  */
226:   for (i = 1; i < user.mt; i++) {
227:     /* Solve the current version */
228:     TaoSolve(tao);

230:     /* Update Vt1 with the solution */
231:     DMGetLocalVector(user.dm, &localX);
232:     DMGlobalToLocalBegin(user.dm, x, INSERT_VALUES, localX);
233:     DMGlobalToLocalEnd(user.dm, x, INSERT_VALUES, localX);
234:     VecGetArray(localX, &x_array);
235:     for (j = 0; j < gxm; j++) user.Vt1[j] = x_array[j];
236:     VecRestoreArray(x, &x_array);
237:     DMRestoreLocalVector(user.dm, &localX);
238:   }

240:   /* Free TAO data structures */
241:   TaoDestroy(&tao);

243:   /* Free PETSc data structures */
244:   VecDestroy(&x);
245:   VecDestroy(&c);
246:   MatDestroy(&J);
247:   DMDestroy(&user.dm);
248:   /* Free user-defined workspace */
249:   PetscFree(user.Vt1);
250:   PetscFree(user.c);
251:   PetscFree(user.d);

253:   PetscFinalize();
254:   return 0;
255: }

257: /* -------------------------------------------------------------------- */
258: PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, void *ctx)
259: {
260:   AppCtx   *user = (AppCtx *)ctx;
261:   PetscInt  i;
262:   PetscInt  xs, xm;
263:   PetscInt  ms   = user->ms;
264:   PetscReal sval = 0.0, *xl_array, ub = PETSC_INFINITY;

266:   /* Set the variable bounds */
267:   VecSet(xu, ub);
268:   DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL);

270:   VecGetArray(xl, &xl_array);
271:   for (i = 0; i < xm; i++) {
272:     sval        = (xs + i) * user->ds;
273:     xl_array[i] = PetscMax(user->strike - sval, 0);
274:   }
275:   VecRestoreArray(xl, &xl_array);

277:   if (xs == 0) {
278:     VecGetArray(xu, &xl_array);
279:     xl_array[0] = PetscMax(user->strike, 0);
280:     VecRestoreArray(xu, &xl_array);
281:   }
282:   if (xs + xm == ms) {
283:     VecGetArray(xu, &xl_array);
284:     xl_array[xm - 1] = 0;
285:     VecRestoreArray(xu, &xl_array);
286:   }

288:   return 0;
289: }
290: /* -------------------------------------------------------------------- */

292: /*
293:     FormFunction - Evaluates gradient of f.

295:     Input Parameters:
296: .   tao  - the Tao context
297: .   X    - input vector
298: .   ptr  - optional user-defined context, as set by TaoAppSetConstraintRoutine()

300:     Output Parameters:
301: .   F - vector containing the newly evaluated gradient
302: */
303: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec F, void *ptr)
304: {
305:   AppCtx    *user = (AppCtx *)ptr;
306:   PetscReal *x, *f;
307:   PetscReal *Vt1 = user->Vt1, *c = user->c, *d = user->d;
308:   PetscReal  rate = user->rate;
309:   PetscReal  dt = user->dt, ds = user->ds;
310:   PetscInt   ms = user->ms;
311:   PetscInt   i, xs, xm, gxs, gxm;
312:   Vec        localX, localF;
313:   PetscReal  zero = 0.0;

315:   DMGetLocalVector(user->dm, &localX);
316:   DMGetLocalVector(user->dm, &localF);
317:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
318:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
319:   DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL);
320:   DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL);
321:   VecSet(F, zero);
322:   /*
323:      The problem size is smaller than the discretization because of the
324:      two fixed elements (V(0,T) = E and V(Send,T) = 0.
325:   */

327:   /* Get pointers to the vector data */
328:   VecGetArray(localX, &x);
329:   VecGetArray(localF, &f);

331:   /* Left Boundary */
332:   if (gxs == 0) {
333:     f[0] = x[0] - user->strike;
334:   } else {
335:     f[0] = 0;
336:   }

338:   /* All points in the interior */
339:   /*  for (i=gxs+1;i<gxm-1;i++) { */
340:   for (i = 1; i < gxm - 1; i++) {
341:     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]);
342:   }

344:   /* Right boundary */
345:   if (gxs + gxm == ms) {
346:     f[xm - 1] = x[gxm - 1];
347:   } else {
348:     f[xm - 1] = 0;
349:   }

351:   /* Restore vectors */
352:   VecRestoreArray(localX, &x);
353:   VecRestoreArray(localF, &f);
354:   DMLocalToGlobalBegin(user->dm, localF, INSERT_VALUES, F);
355:   DMLocalToGlobalEnd(user->dm, localF, INSERT_VALUES, F);
356:   DMRestoreLocalVector(user->dm, &localX);
357:   DMRestoreLocalVector(user->dm, &localF);
358:   PetscLogFlops(24.0 * (gxm - 2));
359:   /*
360:   info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
361:   */
362:   return 0;
363: }

365: /* ------------------------------------------------------------------- */
366: /*
367:    FormJacobian - Evaluates Jacobian matrix.

369:    Input Parameters:
370: .  tao  - the Tao context
371: .  X    - input vector
372: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

374:    Output Parameters:
375: .  J    - Jacobian matrix
376: */
377: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr)
378: {
379:   AppCtx    *user = (AppCtx *)ptr;
380:   PetscReal *c = user->c, *d = user->d;
381:   PetscReal  rate = user->rate;
382:   PetscReal  dt = user->dt, ds = user->ds;
383:   PetscInt   ms = user->ms;
384:   PetscReal  val[3];
385:   PetscInt   col[3];
386:   PetscInt   i;
387:   PetscInt   gxs, gxm;
388:   PetscBool  assembled;

390:   /* Set various matrix options */
391:   MatSetOption(J, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
392:   MatAssembled(J, &assembled);
393:   if (assembled) MatZeroEntries(J);

395:   DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL);

397:   if (gxs == 0) {
398:     i      = 0;
399:     col[0] = 0;
400:     val[0] = 1.0;
401:     MatSetValues(J, 1, &i, 1, col, val, INSERT_VALUES);
402:   }
403:   for (i = 1; i < gxm - 1; i++) {
404:     col[0] = gxs + i - 1;
405:     col[1] = gxs + i;
406:     col[2] = gxs + i + 1;
407:     val[0] = -c[i] / (4 * ds) + d[i] / (2 * ds * ds);
408:     val[1] = 1.0 / dt + rate - d[i] / (ds * ds);
409:     val[2] = c[i] / (4 * ds) + d[i] / (2 * ds * ds);
410:     MatSetValues(J, 1, &col[1], 3, col, val, INSERT_VALUES);
411:   }
412:   if (gxs + gxm == ms) {
413:     i      = ms - 1;
414:     col[0] = i;
415:     val[0] = 1.0;
416:     MatSetValues(J, 1, &i, 1, col, val, INSERT_VALUES);
417:   }

419:   /* Assemble the Jacobian matrix */
420:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
421:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
422:   PetscLogFlops(18.0 * (gxm) + 5);
423:   return 0;
424: }

426: /*TEST

428:    build:
429:       requires: !complex

431:    test:
432:       args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
433:       requires: !single

435:    test:
436:       suffix: 2
437:       args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
438:       requires: !single

440:    test:
441:       suffix: 3
442:       args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5
443:       requires: !single

445:    test:
446:       suffix: 4
447:       args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
448:       requires: !single

450:    test:
451:       suffix: 5
452:       args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5
453:       requires: !single

455:    test:
456:       suffix: 6
457:       args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5
458:       requires: !single

460:    test:
461:       suffix: 7
462:       args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5
463:       requires: !single

465: TEST*/