double x0, x1, f0, f1, xk; //Change variables double tau = 0.5*sigmaSq*time; // f(x,t) -> double dtau= 0.5*sigmaSq*timeStep; // fm(xk,tau) double xk0 = Math.log(x[1]/strike); // lognormal mesh double xkn = Math.log(x[n]/strike); double dxk =(xkn-xk0)/(n-1); double k1 = 2*rate/sigmaSq; double k2 = 2*(rate-divid)/sigmaSq; double k2m1= k2-1.; double k2p1= k2+1.; int j,k; //=== Interpolate only once from (x,t) to lognormal variables (xk,tau) if (time<=timeStep) { if(isPut) { //Initialize fm[] directly as put-option for (k=1; k<=n; k++) { xk=xk0+(k-1)*dxk; fm[k]=Math.max(0., Math.exp(0.5*k2m1*xk) - Math.exp(0.5*k2p1*xk) ); } } else if (isCall) { //Left as an exercise } else { //Interpolate fm[] from IC in f[] j=1; ; x0=xk0; f0=f[1]/strike*Math.exp(0.5*k2m1*xk0); x1=x0; f1=f0; for (k=1; k<n; k++) { //Loop over lognormal mesh index xk=xk0+(k-1)*dxk; // given xk, find x0,x1 | x[j] < xk < x[j+1] while (xk>=x1) { j++; x0=x1; f0=f1; x1=Math.log(x[j]/strike); } f1=f[j]/strike*Math.exp(0.5*k2m1*x1); fm[k]= f0 + (xk-x0)*(f1-f0)/(x1-x0); } fm[n]= fm[n-1] + dxk*(fm[n-1]-fm[n-2]); } } else { //Retrieve fm[] from previous time step } //===== Solve diffusion equation with an explicit 2 levels scheme double D = dtau/(dxk*dxk); for (j=2; j<n; j++) f[j]= fm[j] + D*(fm[j+1]-2.*fm[j]+fm[j-1]); if (isPut) { f[1]= Math.exp(0.5*k2m1*xk0+0.25*k2m1*k2m1*tau); f[n]= f[n-1]; fp[0]=strike*Math.exp(-rate*time); // } else if (isCall) { //Left as an exercise } else { f[1]= f[2]; f[n]= f[n-1]; fp[0]=fp[1]; } //===== Interpolate rest from lognormal to financial mesh variables k=1; x0=x[0]; x1=x0; f0=fp[0]; xk=xk0; f1=f[1]*strike*Math.exp(-0.5*k2m1*xk-(0.25*k2m1*k2m1+k1)*tau); for (j=1; j<n; j++) { //Loop over financial mesh index while (x[j]>=x1){ k++;x0=x1;f0=f1;xk=xk0+(k-1)*dxk;x1=strike*Math.exp(xk);} f1=f[k]*strike*Math.exp(-0.5*k2m1*xk-(0.25*k2m1*k2m1+k1)*tau); fp[j]= f0 +(x[j]-x0)/(x1-x0)*(f1-f0); //Lin interpol in x } if (isPut) { fp[n]=f[n]*strike*Math.exp(-0.5*k2m1*xkn-(0.25*k2m1*k2m1+k1)*tau); // } else if (isCall) { //Left as an exercise } else { fp[0]=fp[1]; fp[n]=fp[n-1]; }
SYLLABUS Previous: 4.4.1 Naive implementation using Up: 4.4 Methods for European Next: 4.5 Methods for European