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