Actual source code: blackscholes.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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 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.

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

    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:   PetscReal *Vt1;                /* Value of the option at time T + dt */
    111:   PetscReal *c;                  /* Constant -- (r - D)S */
    112:   PetscReal *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:   PetscInt  ms, mt;               /* Number of elements */

    122:   DM        dm;
    123: } AppCtx;

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

    127: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
    128: PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
    129: PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void*);

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

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

    149:   /*
    150:      Initialize the user-defined application context with reasonable
    151:      values for the American option to price
    152:   */
    153:   user.rate = 0.04;
    154:   user.sigma = 0.40;
    155:   user.alpha = 2.00;
    156:   user.delta = 0.01;
    157:   user.strike = 10.0;
    158:   user.expiry = 1.0;
    159:   user.mt = 10;
    160:   user.ms = 150;
    161:   user.es = 100.0;

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

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

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

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

    185:   DMCreateGlobalVector(user.dm,&x);
    186:   /*
    187:      Finish filling in the user-defined context with the values for
    188:      dS, dt, and allocating space for the constants
    189:   */
    190:   user.ds = user.es / (user.ms-1);
    191:   user.dt = user.expiry / user.mt;

    193:   PetscMalloc1(gxm,&(user.Vt1));
    194:   PetscMalloc1(gxm,&(user.c));
    195:   PetscMalloc1(gxm,&(user.d));

    197:   /*
    198:      Calculate the values for the constant.  Vt1 begins with the ending
    199:      boundary condition.
    200:   */
    201:   for (i=0; i<gxm; i++) {
    202:     sval = (gxs+i)*user.ds;
    203:     user.Vt1[i] = PetscMax(user.strike - sval, 0);
    204:     user.c[i] = (user.delta - user.rate)*sval;
    205:     user.d[i] = -0.5*user.sigma*user.sigma*PetscPowReal(sval, user.alpha);
    206:   }
    207:   if (gxs+gxm==user.ms){
    208:     user.Vt1[gxm-1] = 0;
    209:   }
    210:   VecDuplicate(x, &c);

    212:   /*
    213:      Allocate the matrix used by TAO for the Jacobian.  Each row of
    214:      the Jacobian matrix will have at most three elements.
    215:   */
    216:   DMCreateMatrix(user.dm,&J);

    218:   /* The TAO code begins here */

    220:   /* Create TAO solver and set desired solution method  */
    221:   TaoCreate(PETSC_COMM_WORLD, &tao);
    222:   TaoSetType(tao,TAOSSILS);

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

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

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

    239:   /* Set routines for function and Jacobian evaluation */
    240:   TaoSetFromOptions(tao);

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

    245:     /* Solve the current version */
    246:     TaoSolve(tao);

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

    260:   /* Free TAO data structures */
    261:   TaoDestroy(&tao);

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

    273:   PetscFinalize();
    274:   return ierr;
    275: }

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

    287:   /* Set the variable bounds */
    288:   VecSet(xu, ub);
    289:   DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL);

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

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

    309:   return 0;
    310: }
    311: /* -------------------------------------------------------------------- */

    313: /*
    314:     FormFunction - Evaluates gradient of f.

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

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

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

    349:   /* Get pointers to the vector data */
    350:   VecGetArray(localX, &x);
    351:   VecGetArray(localF, &f);

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

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

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

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

    388: /* ------------------------------------------------------------------- */
    389: /*
    390:    FormJacobian - Evaluates Jacobian matrix.

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

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

    414:   /* Set various matrix options */
    415:   MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
    416:   MatAssembled(J,&assembled);
    417:   if (assembled){MatZeroEntries(J); }

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

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

    443:   /* Assemble the Jacobian matrix */
    444:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
    445:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
    446:   PetscLogFlops(18*(gxm)+5);
    447:   return 0;
    448: }

    452: /*TEST

    454:    build:
    455:       requires: !complex

    457:    test:
    458:       args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
    459:       requires: !single

    461:    test:
    462:       suffix: 2
    463:       args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
    464:       requires: !single

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

    471:    test:
    472:       suffix: 4
    473:       args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
    474:       requires: !single

    476:    test:
    477:       suffix: 5
    478:       args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5
    479:       requires: !single

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

    486:    test:
    487:       suffix: 7
    488:       args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5
    489:       requires: !single

    491: TEST*/