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[] =
 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";

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

101: /*
102:   User-defined application context - contains data needed by the
103:   application-provided call-back routines, FormFunction(), and FormJacobian().
104: */

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

111:   PetscReal rate;                /* Interest rate */
112:   PetscReal sigma, alpha, delta; /* Underlying asset properties */
113:   PetscReal strike, expiry;      /* Option contract properties */

115:   PetscReal es;                  /* Finite value used for maximum asset value */
116:   PetscReal ds, dt;              /* Discretization properties */
117:   PetscInt  ms, mt;               /* Number of elements */

119:   DM        dm;
120: } AppCtx;

122: /* -------- User-defined Routines --------- */

124: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
125: PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
126: PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void*);

128: int main(int argc, char **argv)
129: {
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;

143:   /* Initialize PETSc, TAO */
144:   PetscInitialize(&argc, &argv, (char *)0, help);if (ierr) return ierr;

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;

160:   /* Read in alternative values for the American option to price */
161:   PetscOptionsGetReal(NULL,NULL, "-alpha", &user.alpha, &flg);
162:   PetscOptionsGetReal(NULL,NULL, "-delta", &user.delta, &flg);
163:   PetscOptionsGetReal(NULL,NULL, "-es", &user.es, &flg);
164:   PetscOptionsGetReal(NULL,NULL, "-expiry", &user.expiry, &flg);
165:   PetscOptionsGetInt(NULL,NULL, "-ms", &user.ms, &flg);
166:   PetscOptionsGetInt(NULL,NULL, "-mt", &user.mt, &flg);
167:   PetscOptionsGetReal(NULL,NULL, "-rate", &user.rate, &flg);
168:   PetscOptionsGetReal(NULL,NULL, "-sigma", &user.sigma, &flg);
169:   PetscOptionsGetReal(NULL,NULL, "-strike", &user.strike, &flg);

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

173:   user.ms++;
174:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.ms,1,1,NULL,&user.dm);
175:   DMSetFromOptions(user.dm);
176:   DMSetUp(user.dm);
177:   /* Create appropriate vectors and matrices */

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

182:   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;

190:   PetscMalloc1(gxm,&(user.Vt1));
191:   PetscMalloc1(gxm,&(user.c));
192:   PetscMalloc1(gxm,&(user.d));

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:   VecDuplicate(x, &c);

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:   DMCreateMatrix(user.dm,&J);

215:   /* The TAO code begins here */

217:   /* Create TAO solver and set desired solution method  */
218:   TaoCreate(PETSC_COMM_WORLD, &tao);
219:   TaoSetType(tao,TAOSSILS);

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

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

228:   /* Set initial solution guess */
229:   VecGetArray(x,&x_array);
230:   for (i=0; i< xm; i++)
231:     x_array[i] = user.Vt1[i-gxs+xs];
232:   VecRestoreArray(x,&x_array);
233:   /* Set data structure */
234:   TaoSetInitialVector(tao, x);

236:   /* Set routines for function and Jacobian evaluation */
237:   TaoSetFromOptions(tao);

239:   /* Iteratively solve the linear complementarity problems  */
240:   for (i = 1; i < user.mt; i++) {

242:     /* Solve the current version */
243:     TaoSolve(tao);

245:     /* Update Vt1 with the solution */
246:     DMGetLocalVector(user.dm,&localX);
247:     DMGlobalToLocalBegin(user.dm,x,INSERT_VALUES,localX);
248:     DMGlobalToLocalEnd(user.dm,x,INSERT_VALUES,localX);
249:     VecGetArray(localX,&x_array);
250:     for (j = 0; j < gxm; j++) {
251:       user.Vt1[j] = x_array[j];
252:     }
253:     VecRestoreArray(x,&x_array);
254:     DMRestoreLocalVector(user.dm,&localX);
255:   }

257:   /* Free TAO data structures */
258:   TaoDestroy(&tao);

260:   /* Free PETSc data structures */
261:   VecDestroy(&x);
262:   VecDestroy(&c);
263:   MatDestroy(&J);
264:   DMDestroy(&user.dm);
265:   /* Free user-defined workspace */
266:   PetscFree(user.Vt1);
267:   PetscFree(user.c);
268:   PetscFree(user.d);

270:   PetscFinalize();
271:   return ierr;
272: }

274: /* -------------------------------------------------------------------- */
275: PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, void*ctx)
276: {
277:   AppCtx         *user = (AppCtx *) ctx;
279:   PetscInt       i;
280:   PetscInt       xs,xm;
281:   PetscInt       ms = user->ms;
282:   PetscReal      sval=0.0,*xl_array,ub= PETSC_INFINITY;

284:   /* Set the variable bounds */
285:   VecSet(xu, ub);
286:   DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL);

288:   VecGetArray(xl,&xl_array);
289:   for (i = 0; i < xm; i++) {
290:     sval = (xs+i)*user->ds;
291:     xl_array[i] = PetscMax(user->strike - sval, 0);
292:   }
293:   VecRestoreArray(xl,&xl_array);

295:   if (xs==0) {
296:     VecGetArray(xu,&xl_array);
297:     xl_array[0] = PetscMax(user->strike, 0);
298:     VecRestoreArray(xu,&xl_array);
299:   }
300:   if (xs+xm==ms) {
301:     VecGetArray(xu,&xl_array);
302:     xl_array[xm-1] = 0;
303:     VecRestoreArray(xu,&xl_array);
304:   }

306:   return 0;
307: }
308: /* -------------------------------------------------------------------- */

310: /*
311:     FormFunction - Evaluates gradient of f.

313:     Input Parameters:
314: .   tao  - the Tao context
315: .   X    - input vector
316: .   ptr  - optional user-defined context, as set by TaoAppSetConstraintRoutine()

318:     Output Parameters:
319: .   F - vector containing the newly evaluated gradient
320: */
321: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec F, void *ptr)
322: {
323:   AppCtx         *user = (AppCtx *) ptr;
324:   PetscReal      *x, *f;
325:   PetscReal      *Vt1 = user->Vt1, *c = user->c, *d = user->d;
326:   PetscReal      rate = user->rate;
327:   PetscReal      dt = user->dt, ds = user->ds;
328:   PetscInt       ms = user->ms;
330:   PetscInt       i, xs,xm,gxs,gxm;
331:   Vec            localX,localF;
332:   PetscReal      zero=0.0;

334:   DMGetLocalVector(user->dm,&localX);
335:   DMGetLocalVector(user->dm,&localF);
336:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
337:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
338:   DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL);
339:   DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL);
340:   VecSet(F, zero);
341:   /*
342:      The problem size is smaller than the discretization because of the
343:      two fixed elements (V(0,T) = E and V(Send,T) = 0.
344:   */

346:   /* Get pointers to the vector data */
347:   VecGetArray(localX, &x);
348:   VecGetArray(localF, &f);

350:   /* Left Boundary */
351:   if (gxs==0) {
352:     f[0] = x[0]-user->strike;
353:   } else {
354:     f[0] = 0;
355:   }

357:   /* All points in the interior */
358:   /*  for (i=gxs+1;i<gxm-1;i++) { */
359:   for (i=1;i<gxm-1;i++) {
360:     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]) +
361:            (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] + Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]);
362:   }

364:   /* Right boundary */
365:   if (gxs+gxm==ms) {
366:     f[xm-1]=x[gxm-1];
367:   } else {
368:     f[xm-1]=0;
369:   }

371:   /* Restore vectors */
372:   VecRestoreArray(localX, &x);
373:   VecRestoreArray(localF, &f);
374:   DMLocalToGlobalBegin(user->dm,localF,INSERT_VALUES,F);
375:   DMLocalToGlobalEnd(user->dm,localF,INSERT_VALUES,F);
376:   DMRestoreLocalVector(user->dm,&localX);
377:   DMRestoreLocalVector(user->dm,&localF);
378:   PetscLogFlops(24.0*(gxm-2));
379:   /*
380:   info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
381:   */
382:   return 0;
383: }

385: /* ------------------------------------------------------------------- */
386: /*
387:    FormJacobian - Evaluates Jacobian matrix.

389:    Input Parameters:
390: .  tao  - the Tao context
391: .  X    - input vector
392: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

394:    Output Parameters:
395: .  J    - Jacobian matrix
396: */
397: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr)
398: {
399:   AppCtx         *user = (AppCtx *) ptr;
400:   PetscReal      *c = user->c, *d = user->d;
401:   PetscReal      rate = user->rate;
402:   PetscReal      dt = user->dt, ds = user->ds;
403:   PetscInt       ms = user->ms;
404:   PetscReal      val[3];
406:   PetscInt       col[3];
407:   PetscInt       i;
408:   PetscInt       gxs,gxm;
409:   PetscBool      assembled;

411:   /* Set various matrix options */
412:   MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
413:   MatAssembled(J,&assembled);
414:   if (assembled) {MatZeroEntries(J);}

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

418:   if (gxs==0) {
419:     i = 0;
420:     col[0] = 0;
421:     val[0]=1.0;
422:     MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
423:   }
424:   for (i=1; i < gxm-1; i++) {
425:     col[0] = gxs + i - 1;
426:     col[1] = gxs + i;
427:     col[2] = gxs + i + 1;
428:     val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds);
429:     val[1] = 1.0/dt + rate - d[i]/(ds*ds);
430:     val[2] =  c[i]/(4*ds) + d[i]/(2*ds*ds);
431:     MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES);
432:   }
433:   if (gxs+gxm==ms) {
434:     i = ms-1;
435:     col[0] = i;
436:     val[0]=1.0;
437:     MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
438:   }

440:   /* Assemble the Jacobian matrix */
441:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
442:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
443:   PetscLogFlops(18.0*(gxm)+5);
444:   return 0;
445: }

447: /*TEST

449:    build:
450:       requires: !complex

452:    test:
453:       args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
454:       requires: !single

456:    test:
457:       suffix: 2
458:       args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
459:       requires: !single

461:    test:
462:       suffix: 3
463:       args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5
464:       requires: !single

466:    test:
467:       suffix: 4
468:       args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
469:       requires: !single

471:    test:
472:       suffix: 5
473:       args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5
474:       requires: !single

476:    test:
477:       suffix: 6
478:       args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5
479:       requires: !single

481:    test:
482:       suffix: 7
483:       args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5
484:       requires: !single

486: TEST*/