Actual source code: blackscholes.c

      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 American 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 with PETSc support.
     69:   Include "petscda.h" so that we can use distributed arrays (DAs) for managing
     70:   the parallel mesh.
     71: */

     73: #include "petscda.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: TaoApplicationCreate(); TaoAppDestroy();
     96:    Routines: TaoAppSetJacobianRoutine(); TaoAppSetConstraintRoutine();
     97:    Routines: TaoAppSetJacobianMat(); TaoSetOptions();
     98:    Routines: TaoSolveApplication();
     99:    Routines: TaoAppSetVariableBoundsRoutine(); TaoAppSetInitialSolutionVec();
    100:    Processors: 1
    101: T*/

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

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

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

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

    122:   DA da;
    123: } AppCtx;

    125: /* -------- User-defined Routines --------- */

    127: int FormConstraints(TAO_APPLICATION, Vec, Vec, void *);
    128: int FormJacobian(TAO_APPLICATION, Vec, Mat *, Mat*, MatStructure*, void *);
    129: int ComputeVariableBounds(TAO_APPLICATION, Vec, Vec, void*);

    133: int main(int argc, char **argv)
    134: {
    135:   int      info;          /* used to check for functions returning nonzeros */
    136:   Vec      x;             /* solution vector */
    137:   Vec      c;             /* Constraints function vector */
    138:   Mat J;                  /* jacobian matrix */
    139:   PetscTruth flg;          /* A return variable when checking for user options */
    140:   TAO_SOLVER tao;          /* TAO_SOLVER solver context */
    141:   TAO_APPLICATION my_app; /* The PETSc application */
    142:   AppCtx user;                  /* user-defined work context */
    143:   int i, j;
    144:   int    xs,xm,gxs,gxm;
    145:   double sval = 0;
    146:   PetscScalar *x_array;
    147:   Vec    localX;

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

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

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

    189:   user.ms++;
    190:   info = DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,user.ms,1,1,
    191:                     PETSC_NULL,&user.da); CHKERRQ(info);
    192:   /* Create appropriate vectors and matrices */

    194:   info = DAGetCorners(user.da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);
    195:   info = DAGetGhostCorners(user.da,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);

    197:   info = DACreateGlobalVector(user.da,&x);CHKERRQ(info);
    198:   /*
    199:      Finish filling in the user-defined context with the values for
    200:      dS, dt, and allocating space for the constants
    201:   */
    202:   user.ds = user.es / (user.ms-1);
    203:   user.dt = user.expiry / user.mt;

    205:   info = PetscMalloc((gxm)*sizeof(double),&(user.Vt1)); CHKERRQ(info);
    206:   info = PetscMalloc((gxm)*sizeof(double),&(user.c)); CHKERRQ(info);
    207:   info = PetscMalloc((gxm)*sizeof(double),&(user.d)); CHKERRQ(info);

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

    223:   info = VecDuplicate(x, &c); CHKERRQ(info);

    225:   /*
    226:      Allocate the matrix used by TAO for the Jacobian.  Each row of
    227:      the Jacobian matrix will have at most three elements.
    228:   */
    229:   info = DAGetMatrix(user.da,MATAIJ,&J);CHKERRQ(info);
    230:
    231:   /* The TAO code begins here */

    233:   /* Create TAO solver and set desired solution method  */
    234:   info = TaoCreate(PETSC_COMM_WORLD, "tao_ssils", &tao); CHKERRQ(info);
    235:   info = TaoApplicationCreate(PETSC_COMM_WORLD, &my_app); CHKERRQ(info);

    237:   /* Set routines for constraints function and Jacobian evaluation */
    238:   info = TaoAppSetConstraintRoutine(my_app, FormConstraints, (void *)&user);
    239:   CHKERRQ(info);
    240:   info = TaoAppSetJacobianRoutine(my_app, FormJacobian, (void *)&user); CHKERRQ(info);
    241:   info = TaoAppSetJacobianMat(my_app, J, J); CHKERRQ(info);
    242:
    243:   /* Set the variable bounds */
    244:   info = TaoAppSetVariableBoundsRoutine(my_app,ComputeVariableBounds,(void*)&user);
    245:   CHKERRQ(info);

    247:   /* Set initial solution guess */
    248:   info = VecGetArray(x,&x_array); CHKERRQ(info);
    249:   for (i=0; i< xm; i++)
    250:     x_array[i] = user.Vt1[i-gxs+xs];
    251:   info = VecRestoreArray(x,&x_array); CHKERRQ(info);
    252:   /* Set data structure */
    253:   info = TaoAppSetInitialSolutionVec(my_app, x); CHKERRQ(info);

    255:   /* Set routines for function and Jacobian evaluation */
    256:   info = TaoSetOptions(my_app,tao); CHKERRQ(info);

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

    261:     /* Solve the current version */
    262:     info = TaoSolveApplication(my_app,tao);  CHKERRQ(info);

    264:     /* Update Vt1 with the solution */
    265:     info = DAGetLocalVector(user.da,&localX);CHKERRQ(info);
    266:     info = DAGlobalToLocalBegin(user.da,x,INSERT_VALUES,localX); CHKERRQ(info);
    267:     info = DAGlobalToLocalEnd(user.da,x,INSERT_VALUES,localX); CHKERRQ(info);
    268:     info = VecGetArray(localX,&x_array); CHKERRQ(info);
    269:     for (j = 0; j < gxm; j++) {
    270:       user.Vt1[j] = x_array[j];
    271:     }
    272:     info = VecRestoreArray(x,&x_array); CHKERRQ(info);
    273:     info = DARestoreLocalVector(user.da,&localX); CHKERRQ(info);
    274:   }

    276:   /* Free TAO data structures */

    278:   info = TaoDestroy(tao); CHKERRQ(info);

    280:   info = TaoAppDestroy(my_app);  CHKERRQ(info);

    282:   /* Free PETSc data structures */
    283:   info = VecDestroy(x); CHKERRQ(info);
    284:   info = VecDestroy(c); CHKERRQ(info);
    285:   info = MatDestroy(J); CHKERRQ(info);
    286:   info = DADestroy(user.da); CHKERRQ(info);
    287:   /* Free user-defined workspace */
    288:   PetscFree(user.Vt1);
    289:   PetscFree(user.c);
    290:   PetscFree(user.d);

    292:   /* Finalize TAO and PETSc */
    293:   PetscFinalize();
    294:   TaoFinalize();

    296:   return 0;
    297: }

    299: /* -------------------------------------------------------------------- */
    302: int ComputeVariableBounds(TAO_APPLICATION tao, Vec xl, Vec xu, void*ctx)
    303: {
    304:   AppCtx *user = (AppCtx *) ctx;
    305:   int i,info;
    306:   int    xs,xm;
    307:   int ms = user->ms;
    308:   PetscScalar sval=0.0,*xl_array,ub= TAO_INFINITY;

    310:   /* Set the variable bounds */
    311:   info = VecSet(xu, ub); CHKERRQ(info);
    312:   info = DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);

    314:   info = VecGetArray(xl,&xl_array); CHKERRQ(info);
    315:   for (i = 0; i < xm; i++){
    316:     sval = (xs+i)*user->ds;
    317:     xl_array[i] = PetscMax(user->strike - sval, 0);
    318:   }
    319:   info = VecRestoreArray(xl,&xl_array); CHKERRQ(info);

    321:   if (xs==0){
    322:     info = VecGetArray(xu,&xl_array); CHKERRQ(info);
    323:     xl_array[0] = PetscMax(user->strike, 0);
    324:     info = VecRestoreArray(xu,&xl_array); CHKERRQ(info);
    325:   }
    326:   if (xs+xm==ms){
    327:     info = VecGetArray(xu,&xl_array); CHKERRQ(info);
    328:     xl_array[xm-1] = 0;
    329:     info = VecRestoreArray(xu,&xl_array); CHKERRQ(info);
    330:   }

    332:   return 0;
    333: }
    334: /* -------------------------------------------------------------------- */

    338: /*
    339:     FormFunction - Evaluates gradient of f.

    341:     Input Parameters:
    342: .   tao  - the TAO_SOLVER context
    343: .   X    - input vector
    344: .   ptr  - optional user-defined context, as set by TaoAppSetConstraintRoutine()
    345:
    346:     Output Parameters:
    347: .   F - vector containing the newly evaluated gradient
    348: */
    349: int FormConstraints(TAO_APPLICATION tao, Vec X, Vec F, void *ptr)
    350: {
    351:   AppCtx *user = (AppCtx *) ptr;
    352:   PetscScalar *x, *f;
    353:   double *Vt1 = user->Vt1, *c = user->c, *d = user->d;
    354:   double rate = user->rate;
    355:   double dt = user->dt, ds = user->ds;
    356:   int ms = user->ms;
    357:   int i, info;
    358:   int    xs,xm,gxs,gxm;
    359:   Vec    localX,localF;
    360:   double zero=0.0;

    362:   info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);
    363:   info = DAGetLocalVector(user->da,&localF);CHKERRQ(info);
    364:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    365:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    366:   info = DAGetCorners(user->da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);
    367:   info = DAGetGhostCorners(user->da,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);
    368:   info = VecSet(F, zero);CHKERRQ(info);
    369:   /*
    370:      The problem size is smaller than the discretization because of the
    371:      two fixed elements (V(0,T) = E and V(Send,T) = 0.
    372:   */

    374:   /* Get pointers to the vector data */
    375:   info = VecGetArray(localX, &x); CHKERRQ(info);
    376:   info = VecGetArray(localF, &f); CHKERRQ(info);
    377:
    378:   /* Left Boundary */
    379:   if (gxs==0){
    380:     f[0] = x[0]-user->strike;
    381:   } else {
    382:     f[0] = 0;
    383:   }

    385:   /* All points in the interior */
    386:   /*  for (i=gxs+1;i<gxm-1;i++){ */
    387:   for (i=1;i<gxm-1;i++){
    388:     f[i] = (1.0/dt + rate)*x[i] - Vt1[i]/dt +
    389:       (c[i]/(4*ds))*(x[i+1] - x[i-1] + Vt1[i+1] - Vt1[i-1]) +
    390:       (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] +
    391:                         Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]);
    392:   }

    394:   /* Right boundary */
    395:   if (gxs+gxm==ms){
    396:     f[xm-1]=x[gxm-1];
    397:   } else {
    398:     f[xm-1]=0;
    399:   }

    401:   /* Restore vectors */
    402:   info = VecRestoreArray(localX, &x); CHKERRQ(info);
    403:   info = VecRestoreArray(localF, &f); CHKERRQ(info);
    404:   info = DALocalToGlobalBegin(user->da,localF,F); CHKERRQ(info);
    405:   info = DALocalToGlobalEnd(user->da,localF,F); CHKERRQ(info);
    406:   info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);
    407:   info = DARestoreLocalVector(user->da,&localF); CHKERRQ(info);
    408:   info = PetscLogFlops(24*(gxm-2)); CHKERRQ(info);
    409:   /*
    410:   info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
    411:   */
    412:   return 0;
    413: }

    415: /* ------------------------------------------------------------------- */
    418: /*
    419:    FormJacobian - Evaluates Jacobian matrix.

    421:    Input Parameters:
    422: .  tao  - the TAO_SOLVER context
    423: .  X    - input vector
    424: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

    426:    Output Parameters:
    427: .  J    - Jacobian matrix
    428: */
    429: int FormJacobian(TAO_APPLICATION tao, Vec X, Mat *tJ, Mat *tJPre, MatStructure *flag, void *ptr)
    430: {
    431:   AppCtx *user = (AppCtx *) ptr;
    432:   Mat J = *tJ;
    433:   double *c = user->c, *d = user->d;
    434:   double rate = user->rate;
    435:   double dt = user->dt, ds = user->ds;
    436:   int ms = user->ms;
    437:   PetscScalar val[3];
    438:   int col[3];
    439:   int i, info;
    440:   int gxs,gxm;
    441:   PetscTruth assembled;

    443:   /* Set various matrix options */
    444:   *flag=SAME_NONZERO_PATTERN;
    445:   info = MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info);
    446:   info = MatSetOption(J,MAT_COLUMNS_SORTED); CHKERRQ(info);
    447:   info = MatSetOption(J,MAT_ROWS_SORTED); CHKERRQ(info);
    448:   info = MatAssembled(J,&assembled); CHKERRQ(info);
    449:   if (assembled){info = MatZeroEntries(J);  CHKERRQ(info);}

    452:   info = DAGetGhostCorners(user->da,&gxs,PETSC_NULL,PETSC_NULL,&gxm,PETSC_NULL,PETSC_NULL); CHKERRQ(info);

    454:   if (gxs==0){
    455:     i = 0;
    456:     col[0] = 0;
    457:     val[0]=1.0;
    458:     info = MatSetValues(J,1,&i,1,col,val,INSERT_VALUES); CHKERRQ(info);
    459:   }
    460:   for (i=1; i < gxm-1; i++) {
    461:     col[0] = gxs + i - 1;
    462:     col[1] = gxs + i;
    463:     col[2] = gxs + i + 1;
    464:     val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds);
    465:     val[1] = 1.0/dt + rate - d[i]/(ds*ds);
    466:     val[2] =  c[i]/(4*ds) + d[i]/(2*ds*ds);
    467:     info = MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES); CHKERRQ(info);
    468:   }
    469:   if (gxs+gxm==ms){
    470:     i = ms-1;
    471:     col[0] = i;
    472:     val[0]=1.0;
    473:     info = MatSetValues(J,1,&i,1,col,val,INSERT_VALUES); CHKERRQ(info);
    474:   }

    476:   /* Assemble the Jacobian matrix */
    477:   info = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    478:   info = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    479:   info = PetscLogFlops(18*(gxm)+5); CHKERRQ(info);
    480:   return 0;
    481: }