Actual source code: blackscholes.c

tao-2.1-p0 2012-07-24
```  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: ***************************************************************************/

67: /*
68:   Include "tao.h" so we can use TAO solvers.
69:   Include "petscdm.h" so that we can use distributed meshes (DMs) for managing
70:   the parallel mesh.
71: */

73: #include "petscdm.h"
74: #include "tao.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: TaoInitialize(); TaoFinalize();
94:    Routines: TaoCreate(); TaoDestroy();
95:    Routines: TaoSetJacobianRoutine(); TaoAppSetConstraintRoutine();
96:    Routines: TaoSetFromOptions();
97:    Routines: TaoSolveApplication();
98:    Routines: TaoSetVariableBoundsRoutine(); TaoSetInitialSolutionVec();
99:    Processors: 1
100: T*/

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

108: typedef struct {
109:   PetscReal *Vt1;                /* Value of the option at time T + dt */
110:   PetscReal *c;                  /* Constant -- (r - D)S */
111:   PetscReal *d;                         /* Constant -- -0.5(sigma**2)(S**alpha) */

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

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

121:   DM dm;
122: } AppCtx;

124: /* -------- User-defined Routines --------- */

126: PetscErrorCode FormConstraints(TaoSolver, Vec, Vec, void *);
127: PetscErrorCode FormJacobian(TaoSolver, Vec, Mat *, Mat*, MatStructure*, void *);
128: PetscErrorCode ComputeVariableBounds(TaoSolver, Vec, Vec, void*);

132: int main(int argc, char **argv)
133: {
135:   Vec      x;             /* solution vector */
136:   Vec      c;             /* Constraints function vector */
137:   Mat J;                  /* jacobian matrix */
138:   PetscBool  flg;          /* A return variable when checking for user options */
139:   TaoSolver tao;          /* TaoSolver solver context */
140:   AppCtx user;                  /* user-defined work context */
141:   PetscInt i, j;
142:   PetscInt    xs,xm,gxs,gxm;
143:   PetscReal sval = 0;
144:   PetscReal *x_array;
145:   Vec    localX;

147:   /* Initialize PETSc, TAO */
148:   PetscInitialize(&argc, &argv, (char *)0, help);
149:   TaoInitialize(&argc, &argv, (char *)0, help);

151:   /*
152:      Initialize the user-defined application context with reasonable
153:      values for the American option to price
154:   */
155:   user.rate = 0.04;
156:   user.sigma = 0.40;
157:   user.alpha = 2.00;
158:   user.delta = 0.01;
159:   user.strike = 10.0;
160:   user.expiry = 1.0;
161:   user.mt = 10;
162:   user.ms = 150;
163:   user.es = 100.0;
164:
165:   /* Read in alternative values for the American option to price */
166:   PetscOptionsGetReal(PETSC_NULL, "-alpha", &user.alpha, &flg);
167:
168:   PetscOptionsGetReal(PETSC_NULL, "-delta", &user.delta, &flg);
169:
170:   PetscOptionsGetReal(PETSC_NULL, "-es", &user.es, &flg);
171:
172:   PetscOptionsGetReal(PETSC_NULL, "-expiry", &user.expiry, &flg);
173:
174:   PetscOptionsGetInt(PETSC_NULL, "-ms", &user.ms, &flg);
175:
176:   PetscOptionsGetInt(PETSC_NULL, "-mt", &user.mt, &flg);
177:
178:   PetscOptionsGetReal(PETSC_NULL, "-rate", &user.rate, &flg);
179:
180:   PetscOptionsGetReal(PETSC_NULL, "-sigma", &user.sigma, &flg);
181:
182:   PetscOptionsGetReal(PETSC_NULL, "-strike", &user.strike, &flg);
183:

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

187:   user.ms++;
188:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,user.ms,1,1,
189:                     PETSC_NULL,&user.dm);
190:   /* Create appropriate vectors and matrices */

192:   DMDAGetCorners(user.dm,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
193:   DMDAGetGhostCorners(user.dm,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL);

195:   DMCreateGlobalVector(user.dm,&x);
196:   /*
197:      Finish filling in the user-defined context with the values for
198:      dS, dt, and allocating space for the constants
199:   */
200:   user.ds = user.es / (user.ms-1);
201:   user.dt = user.expiry / user.mt;

203:   PetscMalloc((gxm)*sizeof(double),&(user.Vt1));
204:   PetscMalloc((gxm)*sizeof(double),&(user.c));
205:   PetscMalloc((gxm)*sizeof(double),&(user.d));

207:   /*
208:      Calculate the values for the constant.  Vt1 begins with the ending
209:      boundary condition.
210:   */
211:   for (i=0; i<gxm; i++) {
212:     sval = (gxs+i)*user.ds;
213:     user.Vt1[i] = PetscMax(user.strike - sval, 0);
214:     user.c[i] = (user.delta - user.rate)*sval;
215:     user.d[i] = -0.5*user.sigma*user.sigma*pow(sval, user.alpha);
216:   }
217:   if (gxs+gxm==user.ms){
218:     user.Vt1[gxm-1] = 0;
219:   }

221:   VecDuplicate(x, &c);

223:   /*
224:      Allocate the matrix used by TAO for the Jacobian.  Each row of
225:      the Jacobian matrix will have at most three elements.
226:   */
227:   DMCreateMatrix(user.dm,MATAIJ,&J);
228:
229:   /* The TAO code begins here */

231:   /* Create TAO solver and set desired solution method  */
232:   TaoCreate(PETSC_COMM_WORLD, &tao);
233:   TaoSetType(tao,"tao_ssils");

235:   /* Set routines for constraints function and Jacobian evaluation */
236:   TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
237:
238:   TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);
239:
240:   /* Set the variable bounds */
241:   TaoSetVariableBoundsRoutine(tao,ComputeVariableBounds,(void*)&user);
242:

244:   /* Set initial solution guess */
245:   VecGetArray(x,&x_array);
246:   for (i=0; i< xm; i++)
247:     x_array[i] = user.Vt1[i-gxs+xs];
248:   VecRestoreArray(x,&x_array);
249:   /* Set data structure */
250:   TaoSetInitialVector(tao, x);

252:   /* Set routines for function and Jacobian evaluation */
253:   TaoSetFromOptions(tao);

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

258:     /* Solve the current version */
259:     TaoSolve(tao);

261:     /* Update Vt1 with the solution */
262:     DMGetLocalVector(user.dm,&localX);
263:     DMGlobalToLocalBegin(user.dm,x,INSERT_VALUES,localX);
264:     DMGlobalToLocalEnd(user.dm,x,INSERT_VALUES,localX);
265:     VecGetArray(localX,&x_array);
266:     for (j = 0; j < gxm; j++) {
267:       user.Vt1[j] = x_array[j];
268:     }
269:     VecRestoreArray(x,&x_array);
270:     DMRestoreLocalVector(user.dm,&localX);
271:   }

273:   /* Free TAO data structures */
274:   TaoDestroy(&tao);

276:   /* Free PETSc data structures */
277:   VecDestroy(&x);
278:   VecDestroy(&c);
279:   MatDestroy(&J);
280:   DMDestroy(&user.dm);
281:   /* Free user-defined workspace */
282:   PetscFree(user.Vt1);
283:   PetscFree(user.c);
284:   PetscFree(user.d);

286:   /* Finalize TAO and PETSc */
287:   PetscFinalize();
288:   TaoFinalize();

290:   return 0;
291: }

293: /* -------------------------------------------------------------------- */
296: PetscErrorCode ComputeVariableBounds(TaoSolver tao, Vec xl, Vec xu, void*ctx)
297: {
298:   AppCtx *user = (AppCtx *) ctx;
300:   PetscInt  i;
301:   PetscInt  xs,xm;
302:   PetscInt  ms = user->ms;
303:   PetscReal sval=0.0,*xl_array,ub= TAO_INFINITY;

305:   /* Set the variable bounds */
306:   VecSet(xu, ub);
307:   DMDAGetCorners(user->dm,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

309:   VecGetArray(xl,&xl_array);
310:   for (i = 0; i < xm; i++){
311:     sval = (xs+i)*user->ds;
312:     xl_array[i] = PetscMax(user->strike - sval, 0);
313:   }
314:   VecRestoreArray(xl,&xl_array);

316:   if (xs==0){
317:     VecGetArray(xu,&xl_array);
318:     xl_array[0] = PetscMax(user->strike, 0);
319:     VecRestoreArray(xu,&xl_array);
320:   }
321:   if (xs+xm==ms){
322:     VecGetArray(xu,&xl_array);
323:     xl_array[xm-1] = 0;
324:     VecRestoreArray(xu,&xl_array);
325:   }

327:   return 0;
328: }
329: /* -------------------------------------------------------------------- */

333: /*
334:     FormFunction - Evaluates gradient of f.

336:     Input Parameters:
337: .   tao  - the TaoSolver context
338: .   X    - input vector
339: .   ptr  - optional user-defined context, as set by TaoAppSetConstraintRoutine()
340:
341:     Output Parameters:
342: .   F - vector containing the newly evaluated gradient
343: */
344: PetscErrorCode FormConstraints(TaoSolver tao, Vec X, Vec F, void *ptr)
345: {
346:   AppCtx *user = (AppCtx *) ptr;
347:   PetscReal *x, *f;
348:   PetscReal *Vt1 = user->Vt1, *c = user->c, *d = user->d;
349:   PetscReal rate = user->rate;
350:   PetscReal dt = user->dt, ds = user->ds;
351:   PetscInt ms = user->ms;
352:   PetscErrorCode    ierr;
353:   PetscInt i, xs,xm,gxs,gxm;
354:   Vec    localX,localF;
355:   PetscReal zero=0.0;

357:   DMGetLocalVector(user->dm,&localX);
358:   DMGetLocalVector(user->dm,&localF);
359:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
360:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
361:   DMDAGetCorners(user->dm,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
362:   DMDAGetGhostCorners(user->dm,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL);
363:   VecSet(F, zero);
364:   /*
365:      The problem size is smaller than the discretization because of the
366:      two fixed elements (V(0,T) = E and V(Send,T) = 0.
367:   */

369:   /* Get pointers to the vector data */
370:   VecGetArray(localX, &x);
371:   VecGetArray(localF, &f);
372:
373:   /* Left Boundary */
374:   if (gxs==0){
375:     f[0] = x[0]-user->strike;
376:   } else {
377:     f[0] = 0;
378:   }

380:   /* All points in the interior */
381:   /*  for (i=gxs+1;i<gxm-1;i++){ */
382:   for (i=1;i<gxm-1;i++){
383:     f[i] = (1.0/dt + rate)*x[i] - Vt1[i]/dt +
384:       (c[i]/(4*ds))*(x[i+1] - x[i-1] + Vt1[i+1] - Vt1[i-1]) +
385:       (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] +
386:                         Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]);
387:   }

389:   /* Right boundary */
390:   if (gxs+gxm==ms){
391:     f[xm-1]=x[gxm-1];
392:   } else {
393:     f[xm-1]=0;
394:   }

396:   /* Restore vectors */
397:   VecRestoreArray(localX, &x);
398:   VecRestoreArray(localF, &f);
399:   DMLocalToGlobalBegin(user->dm,localF,INSERT_VALUES,F);
400:   DMLocalToGlobalEnd(user->dm,localF,INSERT_VALUES,F);
401:   DMRestoreLocalVector(user->dm,&localX);
402:   DMRestoreLocalVector(user->dm,&localF);
403:   PetscLogFlops(24*(gxm-2));
404:   /*
405:   info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
406:   */
407:   return 0;
408: }

410: /* ------------------------------------------------------------------- */
413: /*
414:    FormJacobian - Evaluates Jacobian matrix.

416:    Input Parameters:
417: .  tao  - the TaoSolver context
418: .  X    - input vector
419: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

421:    Output Parameters:
422: .  J    - Jacobian matrix
423: */
424: PetscErrorCode FormJacobian(TaoSolver tao, Vec X, Mat *tJ, Mat *tJPre, MatStructure *flag, void *ptr)
425: {
426:   AppCtx *user = (AppCtx *) ptr;
427:   Mat J = *tJ;
428:   PetscReal *c = user->c, *d = user->d;
429:   PetscReal rate = user->rate;
430:   PetscReal dt = user->dt, ds = user->ds;
431:   PetscInt ms = user->ms;
432:   PetscReal val[3];
434:   PetscInt col[3];
435:   PetscInt i;
436:   PetscInt gxs,gxm;
437:   PetscBool  assembled;

439:   /* Set various matrix options */
440:   *flag=SAME_NONZERO_PATTERN;
441:   MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
442:   MatAssembled(J,&assembled);
443:   if (assembled){MatZeroEntries(J);  }

446:   DMDAGetGhostCorners(user->dm,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL);

448:   if (gxs==0){
449:     i = 0;
450:     col[0] = 0;
451:     val[0]=1.0;
452:     MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
453:   }
454:   for (i=1; i < gxm-1; i++) {
455:     col[0] = gxs + i - 1;
456:     col[1] = gxs + i;
457:     col[2] = gxs + i + 1;
458:     val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds);
459:     val[1] = 1.0/dt + rate - d[i]/(ds*ds);
460:     val[2] =  c[i]/(4*ds) + d[i]/(2*ds*ds);
461:     MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES);
462:   }
463:   if (gxs+gxm==ms){
464:     i = ms-1;
465:     col[0] = i;
466:     val[0]=1.0;
467:     MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
468:   }

470:   /* Assemble the Jacobian matrix */
471:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
472:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
473:   PetscLogFlops(18*(gxm)+5);
474:   return 0;
475: }

```