// 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 <stdio.h>
#include <math.h>
#include <time.h>

#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;
}