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[] = "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 */
133: PetscInitialize(&argc, &argv, (char *)0, 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: PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg);
151: PetscOptionsGetReal(NULL, NULL, "-delta", &user.delta, &flg);
152: PetscOptionsGetReal(NULL, NULL, "-es", &user.es, &flg);
153: PetscOptionsGetReal(NULL, NULL, "-expiry", &user.expiry, &flg);
154: PetscOptionsGetInt(NULL, NULL, "-ms", &user.ms, &flg);
155: PetscOptionsGetInt(NULL, NULL, "-mt", &user.mt, &flg);
156: PetscOptionsGetReal(NULL, NULL, "-rate", &user.rate, &flg);
157: PetscOptionsGetReal(NULL, NULL, "-sigma", &user.sigma, &flg);
158: PetscOptionsGetReal(NULL, NULL, "-strike", &user.strike, &flg);
160: /* Check that the options set are allowable (needs to be done) */
162: user.ms++;
163: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.ms, 1, 1, NULL, &user.dm);
164: DMSetFromOptions(user.dm);
165: DMSetUp(user.dm);
166: /* Create appropriate vectors and matrices */
168: DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL);
169: DMDAGetGhostCorners(user.dm, &gxs, NULL, NULL, &gxm, NULL, NULL);
171: 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: PetscMalloc1(gxm, &(user.Vt1));
180: PetscMalloc1(gxm, &(user.c));
181: 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: 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: DMCreateMatrix(user.dm, &J);
202: /* The TAO code begins here */
204: /* Create TAO solver and set desired solution method */
205: TaoCreate(PETSC_COMM_WORLD, &tao);
206: TaoSetType(tao, TAOSSILS);
208: /* Set routines for constraints function and Jacobian evaluation */
209: TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
210: TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);
212: /* Set the variable bounds */
213: TaoSetVariableBoundsRoutine(tao, ComputeVariableBounds, (void *)&user);
215: /* Set initial solution guess */
216: VecGetArray(x, &x_array);
217: for (i = 0; i < xm; i++) x_array[i] = user.Vt1[i - gxs + xs];
218: VecRestoreArray(x, &x_array);
219: /* Set data structure */
220: TaoSetSolution(tao, x);
222: /* Set routines for function and Jacobian evaluation */
223: TaoSetFromOptions(tao);
225: /* Iteratively solve the linear complementarity problems */
226: for (i = 1; i < user.mt; i++) {
227: /* Solve the current version */
228: TaoSolve(tao);
230: /* Update Vt1 with the solution */
231: DMGetLocalVector(user.dm, &localX);
232: DMGlobalToLocalBegin(user.dm, x, INSERT_VALUES, localX);
233: DMGlobalToLocalEnd(user.dm, x, INSERT_VALUES, localX);
234: VecGetArray(localX, &x_array);
235: for (j = 0; j < gxm; j++) user.Vt1[j] = x_array[j];
236: VecRestoreArray(x, &x_array);
237: DMRestoreLocalVector(user.dm, &localX);
238: }
240: /* Free TAO data structures */
241: TaoDestroy(&tao);
243: /* Free PETSc data structures */
244: VecDestroy(&x);
245: VecDestroy(&c);
246: MatDestroy(&J);
247: DMDestroy(&user.dm);
248: /* Free user-defined workspace */
249: PetscFree(user.Vt1);
250: PetscFree(user.c);
251: PetscFree(user.d);
253: 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: /* Set the variable bounds */
267: VecSet(xu, ub);
268: DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL);
270: VecGetArray(xl, &xl_array);
271: for (i = 0; i < xm; i++) {
272: sval = (xs + i) * user->ds;
273: xl_array[i] = PetscMax(user->strike - sval, 0);
274: }
275: VecRestoreArray(xl, &xl_array);
277: if (xs == 0) {
278: VecGetArray(xu, &xl_array);
279: xl_array[0] = PetscMax(user->strike, 0);
280: VecRestoreArray(xu, &xl_array);
281: }
282: if (xs + xm == ms) {
283: VecGetArray(xu, &xl_array);
284: xl_array[xm - 1] = 0;
285: VecRestoreArray(xu, &xl_array);
286: }
288: return 0;
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: DMGetLocalVector(user->dm, &localX);
316: DMGetLocalVector(user->dm, &localF);
317: DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
318: DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
319: DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL);
320: DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL);
321: VecSet(F, zero);
322: /*
323: The problem size is smaller than the discretization because of the
324: two fixed elements (V(0,T) = E and V(Send,T) = 0.
325: */
327: /* Get pointers to the vector data */
328: VecGetArray(localX, &x);
329: VecGetArray(localF, &f);
331: /* Left Boundary */
332: if (gxs == 0) {
333: f[0] = x[0] - user->strike;
334: } else {
335: f[0] = 0;
336: }
338: /* All points in the interior */
339: /* for (i=gxs+1;i<gxm-1;i++) { */
340: for (i = 1; i < gxm - 1; i++) {
341: 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]);
342: }
344: /* Right boundary */
345: if (gxs + gxm == ms) {
346: f[xm - 1] = x[gxm - 1];
347: } else {
348: f[xm - 1] = 0;
349: }
351: /* Restore vectors */
352: VecRestoreArray(localX, &x);
353: VecRestoreArray(localF, &f);
354: DMLocalToGlobalBegin(user->dm, localF, INSERT_VALUES, F);
355: DMLocalToGlobalEnd(user->dm, localF, INSERT_VALUES, F);
356: DMRestoreLocalVector(user->dm, &localX);
357: DMRestoreLocalVector(user->dm, &localF);
358: PetscLogFlops(24.0 * (gxm - 2));
359: /*
360: info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
361: */
362: return 0;
363: }
365: /* ------------------------------------------------------------------- */
366: /*
367: FormJacobian - Evaluates Jacobian matrix.
369: Input Parameters:
370: . tao - the Tao context
371: . X - input vector
372: . ptr - optional user-defined context, as set by TaoSetJacobian()
374: Output Parameters:
375: . J - Jacobian matrix
376: */
377: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr)
378: {
379: AppCtx *user = (AppCtx *)ptr;
380: PetscReal *c = user->c, *d = user->d;
381: PetscReal rate = user->rate;
382: PetscReal dt = user->dt, ds = user->ds;
383: PetscInt ms = user->ms;
384: PetscReal val[3];
385: PetscInt col[3];
386: PetscInt i;
387: PetscInt gxs, gxm;
388: PetscBool assembled;
390: /* Set various matrix options */
391: MatSetOption(J, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
392: MatAssembled(J, &assembled);
393: if (assembled) MatZeroEntries(J);
395: DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL);
397: if (gxs == 0) {
398: i = 0;
399: col[0] = 0;
400: val[0] = 1.0;
401: MatSetValues(J, 1, &i, 1, col, val, INSERT_VALUES);
402: }
403: for (i = 1; i < gxm - 1; i++) {
404: col[0] = gxs + i - 1;
405: col[1] = gxs + i;
406: col[2] = gxs + i + 1;
407: val[0] = -c[i] / (4 * ds) + d[i] / (2 * ds * ds);
408: val[1] = 1.0 / dt + rate - d[i] / (ds * ds);
409: val[2] = c[i] / (4 * ds) + d[i] / (2 * ds * ds);
410: MatSetValues(J, 1, &col[1], 3, col, val, INSERT_VALUES);
411: }
412: if (gxs + gxm == ms) {
413: i = ms - 1;
414: col[0] = i;
415: val[0] = 1.0;
416: MatSetValues(J, 1, &i, 1, col, val, INSERT_VALUES);
417: }
419: /* Assemble the Jacobian matrix */
420: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
421: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
422: PetscLogFlops(18.0 * (gxm) + 5);
423: return 0;
424: }
426: /*TEST
428: build:
429: requires: !complex
431: test:
432: args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
433: requires: !single
435: test:
436: suffix: 2
437: args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
438: requires: !single
440: test:
441: suffix: 3
442: args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5
443: requires: !single
445: test:
446: suffix: 4
447: args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
448: requires: !single
450: test:
451: suffix: 5
452: args: -tao_monitor -tao_type asils -tao_subset_type matrixfree -pc_type jacobi -tao_max_it 6 -tao_gttol 1.e-5
453: requires: !single
455: test:
456: suffix: 6
457: args: -tao_monitor -tao_type asfls -tao_subset_type subvec -tao_max_it 10 -tao_gttol 1.e-5
458: requires: !single
460: test:
461: suffix: 7
462: args: -tao_monitor -tao_type asfls -tao_subset_type mask -tao_max_it 10 -tao_gttol 1.e-5
463: requires: !single
465: TEST*/