MathFinance - Financial Functions Library in C++

MathFinance - Financial Functions Library in C++

Here you find C++ sample source code and an executalbe file to price spread options. There are several auxiliary functions involving the normal distribution such as
  1. cumulative standard normal distribution functions in one and two dimensions
  2. invers of the cumulative standard normal distribution function.
  3. Guass-Legendre numerical integration
  4. Black-Scholes Formula

Valuation of spread options follwing the article
Low-Fat Spreads by K. Ravindran, RISK, Oct 1993.
download LFspread.exe
view and download spreadinput.txt: This contains the input parameters and the output
The project LFspread is designed to be a Win32 application.
It reads the input data from spreadinput.txt and writes the results into the same file
as well as on the screen.
There are several source files to include in your C/C++ project:

// spread_testenv.cpp
#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include "spread_pricing.h"



void
main(int argc,char**argv)
{
        int greek;
        double s1;
        double s2;
        double r1;
        double r2;
        double rd;
        double v1;
        double v2;
        double rho;
        double T;
        double a;
        double b;
        double K;
        int phi;

        double result;
        
        cout << "SpreadOption" << endl;
        ifstream inputdata("spreadinput.txt");

        inputdata >> s1;
        inputdata >> s2;
        inputdata >> r1;
        inputdata >> r2;
        inputdata >> rd;
        inputdata >> v1;
        inputdata >> v2;
        inputdata >> rho;
        inputdata >> T;
        inputdata >> a;
        inputdata >> b;
        inputdata >> K;
        inputdata >> phi;


        cout << " s1   = " << s1 << endl;
        cout << " s2   = " << s2 << endl;
        cout << " r1   = " << r1 << endl;
        cout << " r2   = " << r2 << endl;
        cout << " rd   = " << rd << endl;
        cout << " v1   = " << v1 << endl;
        cout << " v2   = " << v2 << endl;
        cout << " rho  = " << rho << endl;
        cout << " T    = " << T << endl;
        cout << " a    = " << a << endl;
        cout << " b    = " << b << endl;
        cout << " K    = " << K << endl;
        cout << " phi  = " << phi << endl;


        cout << "spread option value :" << endl;

        ofstream outputdata("spreadinput.txt");

        outputdata <<  s1 << endl;
        outputdata <<  s2 << endl;
        outputdata <<  r1 << endl;
        outputdata <<  r2 << endl;
        outputdata <<  rd << endl;
        outputdata <<  v1 << endl;
        outputdata <<  v2 << endl;
        outputdata <<  rho << endl;
        outputdata <<  T << endl;
        outputdata <<  a << endl;
        outputdata <<  b << endl;
        outputdata <<  K << endl;
        outputdata <<  phi << endl;
        

        for(greek=10;greek<11;greek++)
        { 
                result = LFspread(s1, s2, r1, r2, rd, v1, v2, rho, T, a, b, K, phi, greek);
                cout << result << endl;
                outputdata << result << endl;
        }
        cout << "enter any number to quit ";
        cin >> result;
        
}

//******************************************************* //    file    : spread_pricing.cpp //    author  : Uwe Wystup, http://www.mathfinance.de // //    date    : February 2000 // //******************************************************* #include <math.h> #include "spread_pricing.h" /*_________________________________________________________________________________________*/ /* function :         spread options.                                                      */ /* author:            Uwe Wystup                                                           */ /* contact info:      http://www.mathfinance.de                                            */ /* language:          C++                                                                  */ /* references:        (1) RAVINDRAN, K. Low-Fat Spreads                                    */ /*                        RISK, Oct 1993                                                   */ /*                    (2) WYSTUP, U (2000). The MathFinance Formula Catalogue              */ /*                        http://www.mathfinance.de                                        */ /*_________________________________________________________________________________________*/ double LFspread(double s1, double s2, double r1, double r2, double rd, double \                     v1, double v2, double rho, double T, double \                     a, double b, double K, int phi, int greek) { double result; int k; const int n=41; double x[n]; //these are the abscissas  double w[n]; //these are the weights   result = 0; GaussLegendre(-5.,5., x, w, n-1); // use Gauss-Legendre integration for (k=1;k<n;k++) {     result += w[k] * \     Vanilla(a * s1 * exp(rho * v1 * sqrt(T) * x[k] - 0.5 * v1 * v1 * rho * rho * T), \             b * s2 * exp(v2 * sqrt(T) * x[k] + (rd - r2 - 0.5 * v2 * v2) * T) + K, \             v1 * sqrt(1 - rho * rho), rd, r1, T, phi, 10) * ndf(x[k]); } return result; }
/* file    : vanilla.cpp author  : Uwe Wystup, http://www.mathfinance.de date    : February 2000 */ #include <math.h> #include "normal.h" // function : Black-Scholes formula for European put and call options (vanilla) //            value and Greeks double Vanilla(double spot, double strike, double vol, double rd, double rf, double tau, int \                  phi, int greek) { double result; double d1, d2;          d1 = (log(spot / strike) + (rd - rf + vol * vol * 0.5) * tau) \                 / (vol * sqrt(tau));     d2 = d1 - vol * sqrt(tau);       switch(greek)         {     case 10: // value         result = phi * (spot * exp(-rf * tau) * nc(phi * d1) \                                 - strike * exp(-rd * tau) * nc(phi * d2));                 break;         } return result; }
// GaussLegendre.cpp
#include<math.h>

#define PI 3.141592654
#define EPS 0.00000000000003



// Gauss-Legendre integration
// see Numerical Recipies in C, pg. 152
// by Press, Teukoslky, Vetterling and Flannery
// Cambridge Universtiy Press 1992
//
// Given the lower and uppe limits of integration x1 and x2, and given n,
// this routine returns array x[1..n] and w[1..n] of length n,
// containing the abcissas and weight of the Gauss-Legendre n-point
// quadrature formula.
 
void GaussLegendre(double x1, double x2, double* x, double* w, int n) 
{
  int m,j,i;
  double z1,z,xm,xl,pp,p3,p2,p1;

  m  = (n+1)/2;
  xm = 0.5 * (x2+x1);
  xl = 0.5 * (x2-x1);

  for (i=1;i<=m;i++)
  {
    z = cos(PI*(i-0.25)/(n+0.5));

    do
    {
      p1 = 1.0;
      p2 = 0.0;
      for (j=1;j<=n;j++)
      {
        p3 = p2;
        p2 = p1;
        p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
      }
      pp = n * (z*p1-p2) / (z*z-1.0);
      z1 = z;
      z  = z1 - p1 / pp;
    } while(fabs(z-z1)>EPS);

    x[i]     = xm - xl* z;
    x[n+1-i] = xm + xl* z;

    w[i]     = 2.0 * xl / ((1.0-z*z)*pp*pp);
    w[n+1-i] = w[i];
  }
}





//*******************************************************************
//    program : normal.cpp
//    author  : Uwe Wystup, http://www.mathfinance.de
//
//    created : January 2000
//
//    usage   : auxiliary functions around the normal distribution
//*******************************************************************

#include <math.h>


#define PI 3.141592654


double max(double a, double b)
{
if (a>b)
        return a;
else
        return b;
}


double min(double a, double b)
{
if (a>b)
        return b;
else
        return a;
}

// replicates Sgn as in visual basic, the signum of a real number
double sgn(double a)
{
if (a>0)
        return 1.;
else if (a<0)
        return -1.;
else
        return 0.;
}


// standard normal density function
double ndf(double t)
{
return 0.398942280401433*exp(-t*t/2);
}

// standard normal cumulative distribution function
double nc(double x)
{
double result;
if (x<-7.)
        result = ndf(x)/sqrt(1.+x*x);
else if (x>7.)
        result = 1. - nc(-x);
else
{
result = 0.2316419;
static double a[5] = {0.31938153,-0.356563782,1.781477937,-1.821255978,1.330274429};
result=1./(1+result*fabs(x));
result=1-ndf(x)*(result*(a[0]+result*(a[1]+result*(a[2]+result*(a[3]+result*a[4])))));
if (x<=0.) result=1.-result;
}
return result;
}


//function needed to calculate two dimensional cumulative distribution function, see Hull
double fxy(double x, double y, double a, double b, double rho)
{
    double a_s;
    double b_s;
        double result;
    a_s = a / sqrt(2 * (1 - rho * rho));
    b_s = b / sqrt(2 * (1 - rho * rho));
    result = exp(a_s * (2 * x - a_s) + b_s * (2 * y - b_s) + 2 * rho * (x - a_s) * (y - b_s));
return result;
}


// function needed to calculate two dimensional cumulative distribution function, see Hull
// this equals ND2 if a and b and rho are all nonpositive, 
// the generalization for the other cases is ND2 below
double Ntwo(double a, double b, double rho)
{
    static double aij[4]={0.325303,
                          0.4211071,
                          0.1334425,
                          0.006374323};
        static double bij[4]={0.1337764,
                          0.6243247,
                          1.3425378,
                          2.2626645};
    int i;
    int j;
        double result;
    result = 0;
        for(i=0;i<=3;i++) 
                {
                        for(j=0;j<=3;j++)
                        {
                                result+=aij[i] * aij[j] * fxy(bij[i], bij[j], a, b, rho); 
                        }
                }
    result = result * sqrt(1 - rho * rho) / PI;
return result;
}
//calculates cumulative distribution function for a bivariate normal distribution
//see John Hull: Options, Futures and Other Derivatives
double ND2(double a, double b, double rho)
{
    double rho1;
    double rho2;
    double denominator;
        double result;

    if (rho > 0.9999)
        result = nc(min(a, b));
    else if (rho < -0.9999)
        result = max(0, nc(a) - nc(-b));
    else
        {
        if (a * b * rho <= 0) 
                {
            if (a <= 0 && b <= 0 && rho <= 0)
                result = Ntwo(a, b, rho);
            else if (a <= 0 && b * rho >= 0)
                result = nc(a) - Ntwo(a, -b, -rho);
            else if (b <= 0 && rho >= 0)
                result = nc(b) - Ntwo(-a, b, -rho);
            else
                result = nc(a) + nc(b) - 1 + Ntwo(-a, -b, rho);
                }
        else
                {
            denominator = sqrt(a * a - 2 * rho * a * b + b * b);
            rho1 = (rho * a - b) * sgn(a) / denominator;
            rho2 = (rho * b - a) * sgn(b) / denominator;
            result = ND2(a, 0, rho1) + ND2(b, 0, rho2) - (1 - sgn(a) * sgn(b)) / 4;
        }
        if (result < 0) result = 0;
        }
        return result;
}




double cndev(double u)
{
/* returns the inverse of cumulative normal distribution function
Reference> The Full Monte, by Boris Moro, Union Bank of Switzerland
                        RISK 1995(2)*/

static double a[4]={2.50662823884,
                                        -18.61500062529,
                                        41.39119773534,
                                        -25.44106049637};
static double b[4]={-8.47351093090,
                                        23.08336743743,
                                        -21.06224101826,
                                        3.13082909833};
static double c[9]={0.3374754822726147,
                                        0.9761690190917186,
                                        0.1607979714918209,
                                        0.0276438810333863,
                                        0.0038405729373609,
                                        0.0003951896511919,
                                        0.0000321767881768,
                                        0.0000002888167364,
                                        0.0000003960315187};
double x,r;
x=u-0.5;
if (fabs(x)<0.42)
{
r=x*x;
r=x*(((a[3]*r+a[2])*r+a[1])*r+a[0])/
((((b[3]*r+b[2])*r+b[1])*r+b[0])*r+1.0);
return(r);
}
r=u;
if(x>0.0) r=1.0-u;
        r=log(-log(r));
        r=c[0]+r*(c[1]+r*(c[2]+r*(c[3]+r*(c[4]+r*(c[5]+r*(c[6]+
                r*(c[7]+r*c[8])))))));
        if(x<0.0) r=-r;
return(r);
}




// normal.h
// cumulative distribution function for a bivariate normal distribution, 
// see John Hull: Options, Futures and Other Derivatives
double ND2(double a, double b, double rho);
// standard normal density function
double ndf(double t); 
// standard normal cumulative distribution function
double nc(double x); 
double min(double, double);
double max(double, double);

// Spread_pricing.h
// cumulative distribution function for a bivariate normal distribution, 
// see John Hull: Options, Futures and Other Derivatives
double ND2(double a, double b, double rho);
// standard normal density function
double ndf(double t); 
// standard normal cumulative distribution function
double nc(double x); 
double min(double, double);
double max(double, double);
double LFspread(double s1, double s2, double r1, double r2, double rd, double \
                    v1, double v2, double rho, double T, double \
                    a, double b, double K, int phi, int greek);
double Vanilla(double spot, double strike, double vol, double rd, double rf, double tau, int \
                 phi, int greek);
void GaussLegendre
(double lower, double upper, double* x, double* w, int n);

 






























MathFinance logo Footer Pic 1 Footer Pic 2 Footer Pic 3 Footer Pic 4 © MathFinance AG
Privacy Policy  |  Disclaimer  |  AGB / Terms and Conditions