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 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 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.
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 */
66
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 */
72
73 #include <petscdmda.h>
74 #include <petsctao.h>
75
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";
89
90 /*
91 User-defined application context - contains data needed by the
92 application-provided call-back routines, FormFunction(), and FormJacobian().
93 */
94
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) */
99
100 PetscReal rate; /* Interest rate */
101 PetscReal sigma, alpha, delta; /* Underlying asset properties */
102 PetscReal strike, expiry; /* Option contract properties */
103
104 PetscReal es; /* Finite value used for maximum asset value */
105 PetscReal ds, dt; /* Discretization properties */
106 PetscInt ms, mt; /* Number of elements */
107
108 DM dm;
109 } AppCtx;
110
111 /* -------- User-defined Routines --------- */
112
113 PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
114 PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
115 PetscErrorCode ComputeVariableBounds(Tao, Vec, Vec, void *);
116
main(int argc,char ** argv)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;
130
131 /* Initialize PETSc, TAO */
132 PetscFunctionBeginUser;
133 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
134
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;
148
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));
159
160 /* Check that the options set are allowable (needs to be done) */
161
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 */
167
168 PetscCall(DMDAGetCorners(user.dm, &xs, NULL, NULL, &xm, NULL, NULL));
169 PetscCall(DMDAGetGhostCorners(user.dm, &gxs, NULL, NULL, &gxm, NULL, NULL));
170
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;
178
179 PetscCall(PetscMalloc1(gxm, &user.Vt1));
180 PetscCall(PetscMalloc1(gxm, &user.c));
181 PetscCall(PetscMalloc1(gxm, &user.d));
182
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));
195
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));
201
202 /* The TAO code begins here */
203
204 /* Create TAO solver and set desired solution method */
205 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
206 PetscCall(TaoSetType(tao, TAOSSILS));
207
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));
211
212 /* Set the variable bounds */
213 PetscCall(TaoSetVariableBoundsRoutine(tao, ComputeVariableBounds, (void *)&user));
214
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));
221
222 /* Set routines for function and Jacobian evaluation */
223 PetscCall(TaoSetFromOptions(tao));
224
225 /* Iteratively solve the linear complementarity problems */
226 for (i = 1; i < user.mt; i++) {
227 /* Solve the current version */
228 PetscCall(TaoSolve(tao));
229
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 }
239
240 /* Free TAO data structures */
241 PetscCall(TaoDestroy(&tao));
242
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));
252
253 PetscCall(PetscFinalize());
254 return 0;
255 }
256
257 /* -------------------------------------------------------------------- */
ComputeVariableBounds(Tao tao,Vec xl,Vec xu,PetscCtx ctx)258 PetscErrorCode ComputeVariableBounds(Tao tao, Vec xl, Vec xu, PetscCtx 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;
265
266 PetscFunctionBeginUser;
267 /* Set the variable bounds */
268 PetscCall(VecSet(xu, ub));
269 PetscCall(DMDAGetCorners(user->dm, &xs, NULL, NULL, &xm, NULL, NULL));
270
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));
277
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 /* -------------------------------------------------------------------- */
291
292 /*
293 FormFunction - Evaluates gradient of f.
294
295 Input Parameters:
296 . tao - the Tao context
297 . X - input vector
298 . ptr - optional user-defined context, as set by TaoAppSetConstraintRoutine()
299
300 Output Parameters:
301 . F - vector containing the newly evaluated gradient
302 */
FormConstraints(Tao tao,Vec X,Vec F,void * ptr)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;
314
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 */
327
328 /* Get pointers to the vector data */
329 PetscCall(VecGetArray(localX, &x));
330 PetscCall(VecGetArray(localF, &f));
331
332 /* Left Boundary */
333 if (gxs == 0) {
334 f[0] = x[0] - user->strike;
335 } else {
336 f[0] = 0;
337 }
338
339 /* All points in the interior */
340 /* for (i=gxs+1;i<gxm-1;i++) { */
341 for (i = 1; i < gxm - 1; i++) 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
343 /* Right boundary */
344 if (gxs + gxm == ms) {
345 f[xm - 1] = x[gxm - 1];
346 } else {
347 f[xm - 1] = 0;
348 }
349
350 /* Restore vectors */
351 PetscCall(VecRestoreArray(localX, &x));
352 PetscCall(VecRestoreArray(localF, &f));
353 PetscCall(DMLocalToGlobalBegin(user->dm, localF, INSERT_VALUES, F));
354 PetscCall(DMLocalToGlobalEnd(user->dm, localF, INSERT_VALUES, F));
355 PetscCall(DMRestoreLocalVector(user->dm, &localX));
356 PetscCall(DMRestoreLocalVector(user->dm, &localF));
357 PetscCall(PetscLogFlops(24.0 * (gxm - 2)));
358 /*
359 info=VecView(F,PETSC_VIEWER_STDOUT_WORLD);
360 */
361 PetscFunctionReturn(PETSC_SUCCESS);
362 }
363
364 /* ------------------------------------------------------------------- */
365 /*
366 FormJacobian - Evaluates Jacobian matrix.
367
368 Input Parameters:
369 . tao - the Tao context
370 . X - input vector
371 . ptr - optional user-defined context, as set by TaoSetJacobian()
372
373 Output Parameters:
374 . J - Jacobian matrix
375 */
FormJacobian(Tao tao,Vec X,Mat J,Mat tJPre,void * ptr)376 PetscErrorCode FormJacobian(Tao tao, Vec X, Mat J, Mat tJPre, void *ptr)
377 {
378 AppCtx *user = (AppCtx *)ptr;
379 PetscReal *c = user->c, *d = user->d;
380 PetscReal rate = user->rate;
381 PetscReal dt = user->dt, ds = user->ds;
382 PetscInt ms = user->ms;
383 PetscReal val[3];
384 PetscInt col[3];
385 PetscInt i;
386 PetscInt gxs, gxm;
387 PetscBool assembled;
388
389 PetscFunctionBeginUser;
390 /* Set various matrix options */
391 PetscCall(MatSetOption(J, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
392 PetscCall(MatAssembled(J, &assembled));
393 if (assembled) PetscCall(MatZeroEntries(J));
394
395 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, NULL, NULL, &gxm, NULL, NULL));
396
397 if (gxs == 0) {
398 i = 0;
399 col[0] = 0;
400 val[0] = 1.0;
401 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(J, 1, &i, 1, col, val, INSERT_VALUES));
417 }
418
419 /* Assemble the Jacobian matrix */
420 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
421 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
422 PetscCall(PetscLogFlops(18.0 * (gxm) + 5));
423 PetscFunctionReturn(PETSC_SUCCESS);
424 }
425
426 /*TEST
427
428 build:
429 requires: !complex
430
431 test:
432 args: -tao_monitor -tao_type ssils -tao_gttol 1.e-5
433 requires: !single
434
435 test:
436 suffix: 2
437 args: -tao_monitor -tao_type ssfls -tao_max_it 10 -tao_gttol 1.e-5
438 requires: !single
439
440 test:
441 suffix: 3
442 args: -tao_monitor -tao_type asils -tao_subset_type subvec -tao_gttol 1.e-5
443 requires: !single
444
445 test:
446 suffix: 4
447 args: -tao_monitor -tao_type asils -tao_subset_type mask -tao_gttol 1.e-5
448 requires: !single
449
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
454
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
459
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
464
465 TEST*/
466