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 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: DMGetMatrix(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: }