// LevyLattice.cpp : Defines the entry point for the console application. // Copyright 2002 by Alan L. Lewis // Supplemental code for my article: // American Options under Jump-diffusions: an Introduction // Wilmott Magazine March 2003, 46-51 #include "stdafx.h" #include #include #include #define max(a, b) (((a) > (b)) ? (a) : (b)) #define NJ 2 /* number of Point Jumps */ #define INFINITY 1.0e100 double K = 100.0; double S0 = 100.0; double xj[2] = {log(1.25), log(0.50)}; /* note ordering -- up jump, then down */ double xmin, xmax, dx; #define NX_MAX 256000 /* these array indices will run from 1 to NX_MAX */ double Vold[NX_MAX+1]; double Vnew[NX_MAX+1]; double ExValue[NX_MAX+1]; int NX; /* assumes i runs from 1 to NX */ double x(int i) {return (xmin + (i - 1) * dx);} double Vf(int i) { if( i < 1) return (K - exp(x(i))); else if (i > NX) return 0; else return (Vold[i]); } /* MS replacement for main */ int _tmain(int argc, _TCHAR* argv[]) { float T, q1, r, sig, lambda, cutoff; /* note double doesn't work */ float fl; char lastline[2]; int n, iK, i, j, k, iup, idown; double q[NJ], c[NJ], d[NJ], w[NJ]; int nj[NJ]; double discount, mu, p, dt, ContinuationVal, ContinuationVallast, putval; double slope, scrit; double a1, a2; double tmp; clock_t clock1, clock2; top: printf("Enter Years To Expiration:"); scanf("%f",&T); printf("Enter Number of Time Steps :"); scanf("%i",&n); printf("Enter Up Jump Probability (q1):"); scanf("%f",&q1); printf("Enter r sigma lambda:"); scanf("%f %f %f",&r,&sig,&lambda); printf("Enter xmax = log(Smax/S0) cutoff:"); scanf("%f",&cutoff); printf("You entered T=%8.2f Steps=%d q1=%8.2f r=%8.6f sig=%8.6f lambda=%8.4f \n", T,n,q1,r, sig,lambda); q[0] = q1; /* note q[0] is the up jump */ q[1] = 1 - q1; mu = r - 0.50 * sig * sig; tmp = 0; for(j = 0; j < NJ; j++) tmp+= (exp(xj[j])-1.0)* q[j]; mu-= lambda * tmp; dt = T/n; dx = sig * sqrt(dt); p = 0.5 + 0.5 * mu / sig * sqrt(dt); printf("dt=%10.6g dx=%10.6g p=%10.6g\n",dt,dx,p); /* approximate values -- to be adjusted */ xmin = log(S0) - cutoff; xmax = log(S0) + cutoff; NX = 1 + ceil((xmax - xmin)/dx); if( NX % 2 == 0) NX = NX + 1; /* now odd */ xmin = log(K) - 0.5 * (NX - 1) * dx; xmax = log(K) + 0.5 * (NX - 1) * dx; printf("Tsteps=%d Xsteps=%d Cutoff=%8.2f Smin=%10.6g Smax=%10.6g\n",n,NX,cutoff, exp(xmin),exp(xmax)); if( NX > NX_MAX) { printf("Can't run: NX needs to be less than %d\n", NX_MAX); goto end; } clock1 = clock(); iK = (NX + 1)/2; /* strike index */ for(j = 0; j < NJ; j++) { fl = floor(xj[j]/dx); w[j] = (xj[j] - fl * dx) / dx; /* weight on high point */ nj[j] = fl; /* index of low point */ } /* init tables */ for( i = 1; i <= NX; i++) { Vnew[i] = 0; ExValue[i] = Vold[i] = max(K - exp(x(i)), 0.0); } /* repetitive constants */ discount = exp(-r * dt); a1 = discount* (1.0 - lambda*dt) * p; a2 = discount * (1.0 - lambda*dt)*(1 - p); for (j=0; j < NJ; j++) { c[j] = discount * lambda * dt * w[j] * q[j]; d[j] = discount * lambda * dt * (1.0 - w[j]) * q[j]; } /* time step outer loop */ for (k = 1; k <= n; k++) { /* iterate from strike to upper boundary -- early exercise impossible */ iup = NX; for(i = iK; i <= NX; i++) { tmp = 0.0; for(j = 0; j < NJ; j++) tmp+= d[j] * Vf(i + nj[j]) + c[j] * Vf(i + nj[j] + 1); Vnew[i] = a1 * Vf(i+1) + a2 * Vf(i-1) + tmp; if(Vnew[i] == 0.0) { iup = i - 1; break; } } /* iterate from strike to lower boundary */ idown = 1; for(i = iK - 1; i >= 1; i--) { tmp = 0.0; for(j = 0; j < NJ; j++) tmp+= d[j] * Vf(i + nj[j]) + c[j] * Vf(i + nj[j] + 1); ContinuationVal = a1 * Vf(i+1) + a2 * Vf(i-1) + tmp; if(ContinuationVal < ExValue[i]) { idown = i + 1; break; } else { Vnew[i] = ContinuationVal; ContinuationVallast = ContinuationVal; } } /* only remap entries that have changed */ for (i = idown; i <= iup; i++) Vold[i] = Vnew[i]; } /* end of time step k loop */ putval = Vold[iK]; slope = (ContinuationVallast - ContinuationVal)/(exp(x(idown))- exp(x(idown-1))); scrit = (K - ContinuationVal + exp(x(idown-1))* slope)/(1.0 + slope); printf("At T = %8.2f Scritical = %10.6f ATM Putval = %10.6f\n",T,scrit,putval); clock2 = clock(); printf("\n( LevyLattice.cpp Runtime=%8.2f secs)\n",(double)(clock2-clock1)/CLOCKS_PER_SEC); printf("\n"); end: printf("Enter 'c'(Continue ) or 'x'(Exit)\n"); scanf("%s",lastline); if(lastline[0]== 'x') return 0; else goto top; }