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 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 */
    132:   PetscFunctionBeginUser;
    133:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg));
    151:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta", &user.delta, &flg));
    152:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-es", &user.es, &flg));
    153:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-expiry", &user.expiry, &flg));
    154:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ms", &user.ms, &flg));
    155:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mt", &user.mt, &flg));
    156:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-rate", &user.rate, &flg));
    157:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma", &user.sigma, &flg));
    158:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-strike", &user.strike, &flg));

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

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

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

    171:   PetscCall(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:   PetscCall(PetscMalloc1(gxm, &user.Vt1));
    180:   PetscCall(PetscMalloc1(gxm, &user.c));
    181:   PetscCall(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:   PetscCall(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:   PetscCall(DMCreateMatrix(user.dm, &J));

    202:   /* The TAO code begins here */

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

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

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

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

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

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

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

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

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

    253:   PetscCall(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:   PetscFunctionBeginUser;
    267:   /* Set the variable bounds */
    268:   PetscCall(VecSet(xu, ub));
    269:   PetscCall(DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL));

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

    278:   if (xs == 0) {
    279:     PetscCall(VecGetArray(xu, &xl_array));
    280:     xl_array[0] = PetscMax(user->strike, 0);
    281:     PetscCall(VecRestoreArray(xu, &xl_array));
    282:   }
    283:   if (xs + xm == ms) {
    284:     PetscCall(VecGetArray(xu, &xl_array));
    285:     xl_array[xm - 1] = 0;
    286:     PetscCall(VecRestoreArray(xu, &xl_array));
    287:   }
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    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:   PetscFunctionBeginUser;
    316:   PetscCall(DMGetLocalVector(user->dm, &localX));
    317:   PetscCall(DMGetLocalVector(user->dm, &localF));
    318:   PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
    319:   PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
    320:   PetscCall(DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL));
    321:   PetscCall(DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL));
    322:   PetscCall(VecSet(F, zero));
    323:   /*
    324:      The problem size is smaller than the discretization because of the
    325:      two fixed elements (V(0,T) = E and V(Send,T) = 0.
    326:   */

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

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

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

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

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

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

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

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

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

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

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

    421:   /* Assemble the Jacobian matrix */
    422:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    423:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    424:   PetscCall(PetscLogFlops(18.0 * (gxm) + 5));
    425:   PetscFunctionReturn(PETSC_SUCCESS);
    426: }

    428: /*TEST

    430:    build:
    431:       requires: !complex

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

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

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

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

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

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

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

    467: TEST*/