//        bsfindiffvsanalytic.cpp

// this code runs test output for vanilla put
// comparing the Greeks between the analytical model and the finite-difference method
// 
// finite difference method is theta on log grid where pure diffusion is solved
 
#include <cmath>
#include <iostream>
#include <fstream>
#include "nrutil_nr.h" //allow access to public domain Numerical Recipes matrix routines

using namespace std;

const double phinorm = 0.398942280401433; 

// pdf
double phi(double x)
{
    return phinorm*exp(-x*x/2);
}

/* The following is the compact code provided by Marsaglia in his 2004 paper! 
*/
double PhiMarsaglia(double x)

{long double s=x,t=0,b=x,q=x*x,i=1;
while(s!=t) s=(t=s)+(b*=q/(i+=2));
return 0.5+s*exp(-0.5*q-0.91893853320467274178L)
    ;
}

double BSPut(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d1 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd;
    double d2 = d1-sd;
    return K*exp(-r*t)*PhiMarsaglia(-d2)-S*exp(-q*t)*PhiMarsaglia(-d1);
}

double BSPutDelta(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d1 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd;
    return exp(-q*t)*(PhiMarsaglia(d1)-1.0);
}

double BSPutGamma(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d1 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd;
    return exp(-q*t)*phi(d1)/(S*sd);
}

double BSPutVega(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d1 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd;
    return S*exp(-q*t)*phi(d1)*sqrt(t);
}

double BSPutTheta(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d1 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd;
    double d2 = d1-sd;
    return -S*phi(d1)*exp(-q*t)*sigma/(2*sqrt(t))-q*S*exp(-q*t)*PhiMarsaglia(-d1)+r*K*exp(-r*t)*PhiMarsaglia(-d2);
}

double BSPutRho(double S, double K, double sigma, double r, double q, double t)
{
    double sd = sigma*sqrt(t);
    double d2 =(log(S/K)+(r-q)*t+0.5*sd*sd)/sd-sd;
    return -K*t*exp(-r*t)*PhiMarsaglia(-d2);
}

void tridiagsolve(NRVec<double> &a, NRVec<double> &b, NRVec<double> &c, NRVec<double> &r, NRVec<double> &x)
{
	int j;
	int P=r.size();
	NRVec<double> u(P);
	NRVec<double> d(P);
	NRVec<double> l(P);
	NRVec<double> z(P);
	u[0]=c[0];
	d[0]=b[0];
	z[0]=r[0];
	// Do LU decomposition and forward subsititution as you go
	for (j=1;j<P;j++) {
		u[j]=c[j];
		l[j]=a[j]/d[j-1];
		d[j]=b[j]-l[j]*u[j-1];
		z[j] = r[j]-l[j]*z[j-1];
        }
	//follow with back substitution
	x[P-1] = z[P-1]/d[P-1];
	for (j=(P-2);j>=0;j--){
        x[j] = (z[j] - u[j]*x[j+1])/d[j];
        }
}

double NonDimExpiry(double T, double sigma)
{
       return sigma*sigma*T/2;
}

double kone(double r, double sigma)
{
       return 2.0*r/(sigma*sigma);      
}

double ktwo(double r, double q, double sigma)
{
       return 2.0*(r-q)/(sigma*sigma);      
}

double valmul(double K, double r, double q, double x, double tau, double sigma)
{
       double k1=kone(r,sigma);
       double k2=ktwo(r,q,sigma);
       return K*exp(-(k2-1)*x/2.0 - (k1 + (k2-1)*(k2-1)/4.0)*tau);
}

double upperputbc(double x, double tau, double r, double q, double sigma)
{
       return 0.0;
}

double lowerputbc(double x, double tau, double r, double q, double sigma)
{
       double k2=ktwo(r,q,sigma);
       return exp((k2-1)*x/2 + (k2-1)*(k2-1)*tau/4)-exp((k2+1)*x/2 + (k2+1)*(k2+1)*tau/4) ;
}

double putexercise(double x, double r, double q, double sigma)
{
       double k2=ktwo(r,q,sigma);
       double rez;
       rez = exp((k2-1)*x/2)*(1-exp(x));
       if (rez < 0.0) rez = 0.0; else rez=rez;
       return rez; 
}

void ThetaImplicitDiffusion(int N, int M, double theta, double dt, double dx, double S, double K, double sigma, double r, double q, double t)
{
       int i, j;
       const double PI = 3.141592653589793;
       NRMat<double> solvals(0.0,M+1,2*N+1); //NR dynamic matrix
       double alpha = dt/(dx*dx); // key diffusion parameter
       NRVec<double> rhs(0.0,2*N-1); // create vectors for implicit solver
       NRVec<double> x(0.0,2*N-1);
       NRVec<double> a(0.0,2*N-1);
       NRVec<double> b(0.0,2*N-1);
       NRVec<double> c(0.0,2*N-1);
       NRVec<double> vals(0.0,2*N+1);
       NRVec<double> prevals(0.0,2*N+1);
       NRVec<double> deltas(0.0,2*N+1);
       NRVec<double> gammas(0.0,2*N+1);
       NRVec<double> thetas(0.0,2*N+1);
       
       // These are the Greeks that are obviously reported for one call to algorithm
       // In European case would use rho-delta and vega-gamma identities to do rho and vega
       
       // output alpha
       cout << "alpha = " << alpha << "\n";
       
       // initialize initital data with scaled exercise function
       
       for (j = 0; j <= 2*N; j++)
       {
           solvals[0][j] = putexercise((j-N)*dx,r,q,sigma);
       }
             
       // initialize boundary conditions
       
       for (i = 1; i<= M; i++)
       {
           solvals[i][0] = lowerputbc(-N*dx, i*dt, r, q, sigma);
           solvals[i][2*N] = upperputbc(N*dx, i*dt, r, q, sigma);
           }
       // Run the theta implicit algorithm
       
       // First fill the arrays to feed to implicit solver
       // Note form for theta implicit scheme
       
       for (j=0; j<=2*N-2;j++)
       {
           a[j] = -alpha*theta;
           b[j] = 1.0+2.0*alpha*theta;
           c[j] = -alpha*theta;
       }
       
       for (i= 1; i<= M; i++)
       {
           for (j = 1; j<= 2*N-1; j++)
           {
               rhs[j-1]=(1.0-2.0*alpha*(1-theta))*solvals[i-1][j]+alpha*(1-theta)*(solvals[i-1][j+1]+solvals[i-1][j-1]);
           }
           
           // Now correct the outermost entries of r for the boundary conditions in the difference scheme
           // do so in manner independent of option type;
           
           rhs[0]+=alpha*theta*solvals[i][0];
           rhs[2*N-2]+= alpha*theta*solvals[i][2*N];
 
           tridiagsolve(a, b, c, rhs, x);
           // Note that for a grid with constant parameters as here
           // It would be more efficient to make one call to get L and U
           // before entering the evolution and then
           // only call the fwd and back substitution at each step
           // The version here would be more appropriate for the case when
           // the matrix varied across the entire grid.
           for (j = 1; j<= 2*N-1; j++)
           {
           solvals[i][j] = x[j-1]; 
           }   
       }
    
    for (j=0; j<= 2*N; j++)
        // load up values array with solution at terminal point
        {vals[j] = solvals[M][j];
        // Now convery to real prices by undoing the transformation
        vals[j]=vals[j]*valmul(K, r, q, (j-N)*dx, M*dt, sigma);
        //Repeat at previous time step
        prevals[j] = solvals[M-1][j];
        prevals[j]=prevals[j]*valmul(K, r, q, (j-N)*dx, (M-1)*dt, sigma);
        //Load theta array using bery simple difference based on Euler scheme
        thetas[j] = (prevals[j]-vals[j])/dt;
        // Convert to derivative w.r.t real time
        thetas[j] = thetas[j]*sigma*sigma/2;
        }
       
    ofstream out("bsfdvsanal.txt");
    for (j=0; j<= 2*N; j++)
    {   out.precision(15);
        out << K*exp((j-N)*dx) << " " << vals[j] << " " << BSPut(K*exp((j-N)*dx), K, sigma, r, q, t) << " " << thetas[j] << " " << BSPutTheta(K*exp((j-N)*dx), K, sigma, r, q, t) << "\n";
        }
       return ;
    
}
       


using namespace std;
int main()
{
    double result;
    char name[5];
    int N;
    int M;
    int j;
    double theta;
    double dx;
    double dtau;
    double k1;
    double k2;
    double s;
    double k;
    double vol;
    double r;
    double q;
    double t;
    double bsp;
    double bspd;
    double bspg;
    double bspt;
    double bspv;
    double bspr;

    cout << "Test computations of Black-Scholes pricing and Greeks for Put \n ";
    cout << "Analytical formulae vs finite-difference estimates  \n";
    cout << "Enter number of positive and negative sample points (N) for spatial grid \n ";
    cin >> N;
    cout << "Enter number of sample points (M) for temporal grid \n ";
    cin >> M;
    cout << "Enter dx \n ";
    cin >> dx;
    cout << "Enter dtau \n ";
    cin >> dtau;
    cout << "Enter strike \n ";
    cin >> k;
    cout << "Enter volatility \n ";
    cin >> vol;
    cout << "Enter risk-free rate (continuously compounded) \n ";
    cin >> r;
    cout << "Enter continuous yield (continuously compounded) \n ";
    cin >> q;
    cout << "Enter time remaining to maturity in years \n ";
    cin >> t;
    cout << "Enter theta parameter \n ";
    cout << "0.0 for explicit, 1.0 for full implicit, 0.5 for CN, -1 to trigger Douglas \n ";
    cin >> theta;
    if (theta < 0.0) theta = 0.5-dx*dx/(12.0*dtau); else theta=theta;
    
       
    cout << "Calculating...\n";
    cout << " Results for value and theta being written to bsfdvsanal.txt \n ";
    
    ThetaImplicitDiffusion(N, M, theta, dtau, dx, s, k, vol, r, q, t);

    cout << "Hit any key+<RET> to finish \n ";
    cin >> name;

    return(0);
}