//        blackscholesgreeks.cpp

// code to work out the standard model using
// Marsaglia's normal function
// extended to include Greeks
// this code runs test output for vanilla put
 
#include <cmath>
#include <iostream>
#include <fstream>

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 BSCall(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*exp(-q*t)*PhiMarsaglia(d1)-K*exp(-r*t)*PhiMarsaglia(d2);
}

double BSCallDelta(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);
}

double BSCallGamma(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 BSCallVega(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 BSCallTheta(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 BSCallRho(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);
}

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

using namespace std;
int main()
{
    double result;
    char name[5];
    int N;
    int j;
    double range;
    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 << "Marsaglia's method for Phi \n";
    cout << "Enter number of positive and negative sample points (N) \n ";
    cin >> N;
    cout << "Enter range parameter (range) for prices above and below strike \n ";
    cin >> range;
    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 to maturity in years \n ";
    cin >> t;
    
    
    
    cout << " Result being written to bsgreeksoutput.txt \n " << endl;

    ofstream out("bsgreeksoutput.txt");
    for (j=0; j<= 2*N; j++)
    {  
        s = k -range + j*range/N;
        bsp = BSPut(s,k,vol,r,q,t);
        bspd = BSPutDelta(s,k,vol,r,q,t);
        bspg = BSPutGamma(s,k,vol,r,q,t);
        bspt = BSPutTheta(s,k,vol,r,q,t);
        bspv = BSPutVega(s,k,vol,r,q,t);
        bspr = BSPutRho(s,k,vol,r,q,t);
        out.precision(15);
        out << s << " " << bsp << " " << bspd << " " << bspg << " " << bspt << " " << bspv << " "<< bspr << "\n";
        }
    cout << "Hit any key+<RET> to finish \n ";
    cin >> name;

    return(0);
}