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 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[] =
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*/
101: /*
102: User-defined application context - contains data needed by the
103: application-provided call-back routines, FormFunction(), and FormJacobian().
104: */
106: typedef struct {
107: PetscReal *Vt1; /* Value of the option at time T + dt */
108: PetscReal *c; /* Constant -- (r - D)S */
109: PetscReal *d; /* Constant -- -0.5(sigma**2)(S**alpha) */
111: PetscReal rate; /* Interest rate */
112: PetscReal sigma, alpha, delta; /* Underlying asset properties */
113: PetscReal strike, expiry; /* Option contract properties */
115: PetscReal es; /* Finite value used for maximum asset value */
116: PetscReal ds, dt; /* Discretization properties */
117: PetscInt ms, mt; /* Number of elements */
119: DM dm;
120: } AppCtx;
122: /* -------- User-defined Routines --------- */
124: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
125: PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
126: PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void*);
128: int main(int argc, char **argv)
129: {
131: Vec x; /* solution vector */
132: Vec c; /* Constraints function vector */
133: Mat J; /* jacobian matrix */
134: PetscBool flg; /* A return variable when checking for user options */
135: Tao tao; /* Tao solver context */
136: AppCtx user; /* user-defined work context */
137: PetscInt i, j;
138: PetscInt xs,xm,gxs,gxm;
139: PetscReal sval = 0;
140: PetscReal *x_array;
141: Vec localX;
143: /* Initialize PETSc, TAO */
144: PetscInitialize(&argc, &argv, (char *)0, help);if (ierr) return ierr;
146: /*
147: Initialize the user-defined application context with reasonable
148: values for the American option to price
149: */
150: user.rate = 0.04;
151: user.sigma = 0.40;
152: user.alpha = 2.00;
153: user.delta = 0.01;
154: user.strike = 10.0;
155: user.expiry = 1.0;
156: user.mt = 10;
157: user.ms = 150;
158: user.es = 100.0;
160: /* Read in alternative values for the American option to price */
161: PetscOptionsGetReal(NULL,NULL, "-alpha", &user.alpha, &flg);
162: PetscOptionsGetReal(NULL,NULL, "-delta", &user.delta, &flg);
163: PetscOptionsGetReal(NULL,NULL, "-es", &user.es, &flg);
164: PetscOptionsGetReal(NULL,NULL, "-expiry", &user.expiry, &flg);
165: PetscOptionsGetInt(NULL,NULL, "-ms", &user.ms, &flg);
166: PetscOptionsGetInt(NULL,NULL, "-mt", &user.mt, &flg);
167: PetscOptionsGetReal(NULL,NULL, "-rate", &user.rate, &flg);
168: PetscOptionsGetReal(NULL,NULL, "-sigma", &user.sigma, &flg);
169: PetscOptionsGetReal(NULL,NULL, "-strike", &user.strike, &flg);
171: /* Check that the options set are allowable (needs to be done) */
173: user.ms++;
174: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.ms,1,1,NULL,&user.dm);
175: DMSetFromOptions(user.dm);
176: DMSetUp(user.dm);
177: /* Create appropriate vectors and matrices */
179: DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL);
180: DMDAGetGhostCorners(user.dm,&gxs,NULL,NULL,&gxm,NULL,NULL);
182: DMCreateGlobalVector(user.dm,&x);
183: /*
184: Finish filling in the user-defined context with the values for
185: dS, dt, and allocating space for the constants
186: */
187: user.ds = user.es / (user.ms-1);
188: user.dt = user.expiry / user.mt;
190: PetscMalloc1(gxm,&(user.Vt1));
191: PetscMalloc1(gxm,&(user.c));
192: PetscMalloc1(gxm,&(user.d));
194: /*
195: Calculate the values for the constant. Vt1 begins with the ending
196: boundary condition.
197: */
198: for (i=0; i<gxm; i++) {
199: sval = (gxs+i)*user.ds;
200: user.Vt1[i] = PetscMax(user.strike - sval, 0);
201: user.c[i] = (user.delta - user.rate)*sval;
202: user.d[i] = -0.5*user.sigma*user.sigma*PetscPowReal(sval, user.alpha);
203: }
204: if (gxs+gxm==user.ms) {
205: user.Vt1[gxm-1] = 0;
206: }
207: VecDuplicate(x, &c);
209: /*
210: Allocate the matrix used by TAO for the Jacobian. Each row of
211: the Jacobian matrix will have at most three elements.
212: */
213: DMCreateMatrix(user.dm,&J);
215: /* The TAO code begins here */
217: /* Create TAO solver and set desired solution method */
218: TaoCreate(PETSC_COMM_WORLD, &tao);
219: TaoSetType(tao,TAOSSILS);
221: /* Set routines for constraints function and Jacobian evaluation */
222: TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
223: TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);
225: /* Set the variable bounds */
226: TaoSetVariableBoundsRoutine(tao,ComputeVariableBounds,(void*)&user);
228: /* Set initial solution guess */
229: VecGetArray(x,&x_array);
230: for (i=0; i< xm; i++)
231: x_array[i] = user.Vt1[i-gxs+xs];
232: VecRestoreArray(x,&x_array);
233: /* Set data structure */
234: TaoSetInitialVector(tao, x);
236: /* Set routines for function and Jacobian evaluation */
237: TaoSetFromOptions(tao);
239: /* Iteratively solve the linear complementarity problems */
240: for (i = 1; i < user.mt; i++) {
242: /* Solve the current version */
243: TaoSolve(tao);
245: /* Update Vt1 with the solution */
246: DMGetLocalVector(user.dm,&localX);
247: DMGlobalToLocalBegin(user.dm,x,INSERT_VALUES,localX);
248: DMGlobalToLocalEnd(user.dm,x,INSERT_VALUES,localX);
249: VecGetArray(localX,&x_array);
250: for (j = 0; j < gxm; j++) {
251: user.Vt1[j] = x_array[j];
252: }
253: VecRestoreArray(x,&x_array);
254: DMRestoreLocalVector(user.dm,&localX);
255: }
257: /* Free TAO data structures */
258: TaoDestroy(&tao);
260: /* Free PETSc data structures */
261: VecDestroy(&x);
262: VecDestroy(&c);
263: MatDestroy(&J);
264: DMDestroy(&user.dm);
265: /* Free user-defined workspace */
266: PetscFree(user.Vt1);
267: PetscFree(user.c);
268: PetscFree(user.d);
270: PetscFinalize();
271: return ierr;
272: }
274: /* -------------------------------------------------------------------- */
275: PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, void*ctx)
276: {
277: AppCtx *user = (AppCtx *) ctx;
279: PetscInt i;
280: PetscInt xs,xm;
281: PetscInt ms = user->ms;
282: PetscReal sval=0.0,*xl_array,ub= PETSC_INFINITY;
284: /* Set the variable bounds */
285: VecSet(xu, ub);
286: DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL);
288: VecGetArray(xl,&xl_array);
289: for (i = 0; i < xm; i++) {
290: sval = (xs+i)*user->ds;
291: xl_array[i] = PetscMax(user->strike - sval, 0);
292: }
293: VecRestoreArray(xl,&xl_array);
295: if (xs==0) {
296: VecGetArray(xu,&xl_array);
297: xl_array[0] = PetscMax(user->strike, 0);
298: VecRestoreArray(xu,&xl_array);
299: }
300: if (xs+xm==ms) {
301: VecGetArray(xu,&xl_array);
302: xl_array[xm-1] = 0;
303: VecRestoreArray(xu,&xl_array);
304: }
306: return 0;
307: }
308: /* -------------------------------------------------------------------- */
310: /*
311: FormFunction - Evaluates gradient of f.
313: Input Parameters:
314: . tao - the Tao context
315: . X - input vector
316: . ptr - optional user-defined context, as set by TaoAppSetConstraintRoutine()
318: Output Parameters:
319: . F - vector containing the newly evaluated gradient
320: */
321: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec F, void *ptr)
322: {
323: AppCtx *user = (AppCtx *) ptr;
324: PetscReal *x, *f;
325: PetscReal *Vt1 = user->Vt1, *c = user->c, *d = user->d;
326: PetscReal rate = user->rate;
327: PetscReal dt = user->dt, ds = user->ds;
328: PetscInt ms = user->ms;
330: PetscInt i, xs,xm,gxs,gxm;
331: Vec localX,localF;
332: PetscReal zero=0.0;
334: DMGetLocalVector(user->dm,&localX);
335: DMGetLocalVector(user->dm,&localF);
336: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
337: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
338: DMDAGetCorners(user->dm,&xs,NULL,NULL,&xm,NULL,NULL);
339: DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL);
340: VecSet(F, zero);
341: /*
342: The problem size is smaller than the discretization because of the
343: two fixed elements (V(0,T) = E and V(Send,T) = 0.
344: */
346: /* Get pointers to the vector data */
347: VecGetArray(localX, &x);
348: VecGetArray(localF, &f);
350: /* Left Boundary */
351: if (gxs==0) {
352: f[0] = x[0]-user->strike;
353: } else {
354: f[0] = 0;
355: }
357: /* All points in the interior */
358: /* for (i=gxs+1;i<gxm-1;i++) { */
359: for (i=1;i<gxm-1;i++) {
360: 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]) +
361: (d[i]/(2*ds*ds))*(x[i+1] -2*x[i] + x[i-1] + Vt1[i+1] - 2*Vt1[i] + Vt1[i-1]);
362: }
364: /* Right boundary */
365: if (gxs+gxm==ms) {
366: f[xm-1]=x[gxm-1];
367: } else {
368: f[xm-1]=0;
369: }
371: /* Restore vectors */
372: VecRestoreArray(localX, &x);
373: VecRestoreArray(localF, &f);
374: DMLocalToGlobalBegin(user->dm,localF,INSERT_VALUES,F);
375: DMLocalToGlobalEnd(user->dm,localF,INSERT_VALUES,F);
376: DMRestoreLocalVector(user->dm,&localX);
377: DMRestoreLocalVector(user->dm,&localF);
378: PetscLogFlops(24.0*(gxm-2));
379: /*
380: info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
381: */
382: return 0;
383: }
385: /* ------------------------------------------------------------------- */
386: /*
387: FormJacobian - Evaluates Jacobian matrix.
389: Input Parameters:
390: . tao - the Tao context
391: . X - input vector
392: . ptr - optional user-defined context, as set by TaoSetJacobian()
394: Output Parameters:
395: . J - Jacobian matrix
396: */
397: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr)
398: {
399: AppCtx *user = (AppCtx *) ptr;
400: PetscReal *c = user->c, *d = user->d;
401: PetscReal rate = user->rate;
402: PetscReal dt = user->dt, ds = user->ds;
403: PetscInt ms = user->ms;
404: PetscReal val[3];
406: PetscInt col[3];
407: PetscInt i;
408: PetscInt gxs,gxm;
409: PetscBool assembled;
411: /* Set various matrix options */
412: MatSetOption(J,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
413: MatAssembled(J,&assembled);
414: if (assembled) {MatZeroEntries(J);}
416: DMDAGetGhostCorners(user->dm,&gxs,NULL,NULL,&gxm,NULL,NULL);
418: if (gxs==0) {
419: i = 0;
420: col[0] = 0;
421: val[0]=1.0;
422: MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
423: }
424: for (i=1; i < gxm-1; i++) {
425: col[0] = gxs + i - 1;
426: col[1] = gxs + i;
427: col[2] = gxs + i + 1;
428: val[0] = -c[i]/(4*ds) + d[i]/(2*ds*ds);
429: val[1] = 1.0/dt + rate - d[i]/(ds*ds);
430: val[2] = c[i]/(4*ds) + d[i]/(2*ds*ds);
431: MatSetValues(J,1,&col[1],3,col,val,INSERT_VALUES);
432: }
433: if (gxs+gxm==ms) {
434: i = ms-1;
435: col[0] = i;
436: val[0]=1.0;
437: MatSetValues(J,1,&i,1,col,val,INSERT_VALUES);
438: }
440: /* Assemble the Jacobian matrix */
441: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
442: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
443: PetscLogFlops(18.0*(gxm)+5);
444: return 0;
445: }
447: /*TEST
449: build:
450: requires: !complex
452: test:
453: args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
454: requires: !single
456: test:
457: suffix: 2
458: args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
459: requires: !single
461: test:
462: suffix: 3
463: args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5
464: requires: !single
466: test:
467: suffix: 4
468: args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
469: requires: !single
471: test:
472: suffix: 5
473: args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5
474: requires: !single
476: test:
477: suffix: 6
478: args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5
479: requires: !single
481: test:
482: suffix: 7
483: args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5
484: requires: !single
486: TEST*/