asianpricer.cpp

Quantitative Finance C++ source code, asian option pricer

#include "stdafx.h"
#include <complex>
#include <iostream>

using namespace std;

const double M_PI=3.141592653589793238462643383280;

extern void HyperGeomFun(double a_real,double a_imag,double b_real,double b_imag,
                    double z_real,double z_imag,double* res_real,double* res_imag);
extern void cgamma(double arg_real,double arg_imag,double* res_real,
						double* res_imag,int logFlag);


complex<double> tau_fn(complex<double> T,complex<double> t,complex<double> sigma)
{
  return complex<double>(double(1)/double(4),0)*sigma*sigma*(T-t);
}

complex<double> v_fn(complex<double> r,complex<double> q,complex<double> sigma)
{
  return complex<double>(2,0)*(r-q)/(sigma*sigma) - complex<double>(1,0);
}

complex<double> alpha_fn(complex<double> S,complex<double> ES,complex<double> K,
                         complex<double> sigma, complex<double> T,complex<double> t,
                         complex<double> t0)
{
  return sigma*sigma*(K*(T-t0)-(t-t0)*ES)/(complex<double>(4,0)*S);
}

complex<double> mu_fn(complex<double> n,complex<double> p)
{
  return pow(n*n+complex<double>(2,0)*p,0.5);
}


complex<double> LaplaceTransform(complex<double> p,complex<double> mu,
                                 complex<double> v,complex<double> alpha)
{
  complex<double> res;
  complex<double> a;
  complex<double> b;
  complex<double> z;
  double res_real,res_imag;
  a=(mu+v+complex<double>(4,0))/complex<double>(2,0);
  b=mu+complex<double>(1,0);
  z=complex<double>(1,0)/(complex<double>(2,0)*alpha);
  HyperGeomFun(a.real(),a.imag(),
    b.real(),b.imag(),
    z.real(),z.imag(),
    &res_real,&res_imag);
  
  complex<double> tmpcgamma1,tmphypergeom,tmpcgamma2,tmpdenom,tmp5;
  complex<double> ba=(mu+v+complex<double>(4,0))/complex<double>(2,0);
  double res_real2;
  double res_imag2;
  cgamma(ba.real(),ba.imag(),&res_real2,&res_imag2,1);
  tmpcgamma1=complex<double>(res_real2,res_imag2);			
  
  
  tmphypergeom=complex<double>(res_real,res_imag);
  
  double res_real3;
  double res_imag3;
  complex<double> tmp4in=mu+complex<double>(1,0);
  cgamma(tmp4in.real(),tmp4in.imag(),&res_real3,&res_imag3,1);
  tmpcgamma2=complex<double>(res_real3,res_imag3);
  
  tmp5=-z+((v+complex<double>(2,0)-mu)/complex<double>(2,0))*log(complex<double>(2,0)*alpha);
  tmpdenom=log(p*(p-complex<double>(2,0)*(v+complex<double>(1,0))));
  res=tmpcgamma1+tmphypergeom-tmpcgamma2+tmp5-tmpdenom;
  return res;
}



class CParams
{
public:
  complex<double> contour;
  complex<double> v;
  complex<double> alpha;
  complex<double> ti;
};

CParams paramsObj;

double KrondonIntegrand(double x)
{
  complex<double> contour=paramsObj.contour;
  complex<double> p=complex<double>(contour.real(),x);
  complex<double> v=paramsObj.v;
  complex<double> mu=mu_fn(v,p);
  complex<double> alpha=paramsObj.alpha;
  complex<double> ti=paramsObj.ti;
  complex<double> zi(0,p.imag());
  complex<double> g2 = LaplaceTransform(p,mu,v,alpha);
  complex<double> tmptest=exp(zi*ti + g2);
  complex<double> res=
    exp(zi*ti + g2) / (complex<double>(M_PI*2,0))  ;
  return res.real();
}


typedef double ( * lpFunc)(double);
int functionEvaluations_=0;

double GaussKronrod( lpFunc f,
                    const double a, const double b,
                    const double tolerance)  {
  
    static const double g7w[] = { 0.417959183673469,
                                  0.381830050505119, 
                                  0.279705391489277, 
                                  0.129484966168870 };
    static const double k15w[] = { 0.209482141084728, 
                                   0.204432940075298,
                                   0.190350578064785, 
                                   0.169004726639267,
                                   0.140653259715525, 
                                   0.104790010322250,
                                   0.063092092629979, 
                                   0.022935322010529 };
    static const double k15t[] = { 0.000000000000000, 
                                   0.207784955007898,
                                   0.405845151377397, 
                                   0.586087235467691,
                                   0.741531185599394, 
                                   0.864864423359769,
                                   0.949107912342758, 
                                   0.991455371120813 };

    const double halflength = (b - a) / 2;
    const double center = (a + b) / 2;

    double g7; 
    double k15; 

    double t, fsum; 
    double fc = f(center);
    g7 = fc * g7w[0];
    k15 = fc * k15w[0];
    int j, j2;
    for (j = 1, j2 = 2; j < 4; j++, j2 += 2) {
        t = halflength * k15t[j2];
        fsum = f(center - t) + f(center + t);
        g7  += fsum * g7w[j];
        k15 += fsum * k15w[j2];
    }
    for (j2 = 1; j2 < 8; j2 += 2) {
        t = halflength * k15t[j2];
        fsum = f(center - t) + f(center + t);
        k15 += fsum * k15w[j2];
    }
    g7 = halflength * g7;
    k15 = halflength * k15;
    functionEvaluations_ += 15;

    if (abs(k15 - g7) < tolerance) {
        return k15;
    } else {
        cout << "maximum number of function evaluations exceeded" << endl;
        return GaussKronrod(f, a, center, tolerance/2)
             + GaussKronrod(f, center, b, tolerance/2);
    }
}


double AsianPrice(double S,double ES,double K,double r,double q,
        double sigma,double T,double t,double t0,double truncation)
{
  complex<double> ti=tau_fn(complex<double>(T,0),complex<double>(t,0),complex<double>(sigma,0));
  complex<double> v=v_fn(complex<double>(r,0),complex<double>(q,0),complex<double>(sigma,0));
  complex<double> n=v;
  complex<double> alpha=alpha_fn(complex<double>(S,0),complex<double>(ES,0),complex<double>(K,0),complex<double>(sigma,0),complex<double>(T,0),complex<double>(t,0),complex<double>(t0,0));
  double contour=((complex<double>(2,0)*n+complex<double>(2,0))/complex<double>(sigma*sigma,0)).real();
  paramsObj.contour=complex<double>(contour,0);
  paramsObj.v=v;
  paramsObj.alpha=alpha;
  paramsObj.ti=ti;
  
  double valsum=0;
  for (int tmpi=0;tmpi<int(truncation);tmpi=tmpi+10)
  {
    valsum+=GaussKronrod(KrondonIntegrand,tmpi,tmpi+10,1e-20);
  }
  double h=ti.real();
  complex<double> alpha2=(complex<double>(2,0)*v+complex<double>(2,0))/(complex<double>(sigma,0)*complex<double>(sigma,0));
  double asainPrice=exp(-r*T)*(S/h)*exp(alpha2.real()*h)*2*valsum;
  
  return asainPrice;
}



int main(int argc, char* argv[])
{
  double asianPrice=AsianPrice(2,0,2,0.05,0,0.5,1,0,0,1000);
  cout << "price=" << asianPrice << endl;
  return 0;
}