/* A program to compute the early exercise boundary for American style stock put options. Written by Jonathan Goodman, Courant Institute of Mathematical Sciences, NYU goodman@cims.nyu.edu http://www.math.nyu.edu/faculty/goodman Dan Ostrov, Mathematics Department, Santa Clara University dostrov@scu.edu written summer 2000. This was written to produce the pictures in our paper on the early exercise boundary for American style puts. It is not particularly efficient and is not "state of the art" computational technique. */ #include #include /* Will use the log or exp functions somewhere */ #include /* Use the C++ file communication system for dumping computed results. */ #include /* To be able to print more decimal digits in the output*/ #define NX_MAX 16002 /* A dimension bigger than we will probably ever use to simplify storage allocation */ #define NUMAVG 500 #define T_STEP_RATIO .9 /* Thr fraction of max possible CFL ratio to use. */ #define MACH_EPSILON 1.e-15 /* A small tolerence related to the roundoff level for double precision computer arithmetic. Used in time_step in taking the max of the computed solution and the computed intrinsic value. */ // A program to solve the black scholes problem for an american style put // in log variables. //*************************** time_step ************************************** /* Take a single time step for the american put problem. First take a Forward Euler step, using precomputed coefficients a, b, and c. Then take the maximium of that with the intrinsic value. Report an error if the early exercise boundary is in the first grid interval. Apply dirichlet boundary conditions at the right boundary. The strike price is taken to be S=1, which corresponds to x=0 in log variables. u is the solution being computed. Its first and last entries hold boundary values and are not updated. The first entry is u[0]. The last is u[nx]. The u array has nx+1 entries, of which nx-1 are updated. The old values are overwritten by new values. The routine returns a number, x_crit, which is the x value of the first grid point where the constraint u>1-S = 1-exp(x) is not binding. */ static int time_step( double u[], double uintr[], double x_crit[], double x_min, double dx, double dt, double sigma, double r, int nx) { double* u_new= new double[NX_MAX]; // Temporary storage for the time step routine. // Coefficients for the finite difference method. double a = (dt*sigma*sigma)/(2*dx*dx) - ( ( r - sigma*sigma/2 )*dt)/(2*dx); double b = 1 - ( dt*sigma*sigma)/(dx*dx) - r*dt; double c = (dt*sigma*sigma)/(2*dx*dx) + ( ( r - sigma*sigma/2 )*dt)/(2*dx); // A safety precaution, probably unnecessary. if ( ( a < 0 ) || ( b < 0 ) || ( c < 0 ) ) { cout << "in time_step, bailing because a coefficient is negative." << " a, b, c, = " << a << " " << b << " " << c << endl; delete u_new; return 1; } // Take a time step for the PDE int ix; for ( ix=1; ix uintr[ix]*( 1 + MACH_EPSILON) ) // Allow for roundoff in the u_new computation. break; // Stop if you've left the intrinsic value curve . . . u_new[ix] = uintr[ix]; // Otherwise make correct the numerical solution. } // end of "for (int ix=1; ix= NX_MAX ) { cout << "In main program, not big enough arrays. \n" << " asked for nx = " << nx << ". Had NX_MAX = " << NX_MAX << endl; return(3); } double* u= new double[NX_MAX]; double* uintr= new double[NX_MAX]; double x_crit; double xcritical[NUMAVG]; double time[NUMAVG]; //Initialize the solution double x = x_min; double x_cutoff = -1.e-3; // Purely technical: change how the // intrinsic value is computed. int ix; for ( ix = 0; ix <= nx; ix++) { // This loop includes both end points. if ( x > 0 ) u[ix] = 0.; else if ( x < x_cutoff ) // The argument is large enough that the formula u[ix] = 1 - exp(x) ; // is accurate enough. else u[ix] = - x * ( 1 + x/2 * ( 1 + x/3 * ( 1 + x/4 ) ) ); // This is a factored form of the Taylor series for 1-exp(x), // which has no cancellation and therefore should be more accurate // when x is small. It should give full double precision accuracy. // when |x| < 10^-4. x = x + dx; // The first trip through the loop has ix=0, so the // increment goes at the end. } // end of "for ( ix = 0; ix <= nx; ix++) ... " for ( ix = 0; ix <= nx; ix++) uintr[ix]=u[ix]; //End of initialization double t = 0; int it, error_code; double vsum = 0; // Used in comparing computed x_crit to theory. for ( it = 0; t < T; it++ ) { t += dt; error_code = time_step( u, uintr, &x_crit, x_min, dx, dt, sigma, r, nx); xcritical[it % NUMAVG] = x_crit; time[it % NUMAVG]=t; if ( error_code != 0 ) { cout << "Main returning because of time step error code " << error_code << endl; cout << " it = " << it << " dt = " << dt; return(1); } } // This is the fun part, seeing how the data fit the theory . . . int count; for (count = 0; count < NUMAVG; count++) vsum += exp(-xcritical[count]*xcritical[count]/(sigma*sigma*time[count]))/time[count]; crit_price_output << t << " " << x_crit << " " << vsum/numavgd << " " << it << endl; cout << t <<" " << it << endl; } crit_price_output.close(); cout << "done! " << endl; return 0; }