Files

main.cpp

#include <iostream>
#include <cmath>

using namespace std;


int recursivelySumBetween1AndN( int n ) {
    if (n==1) {
        return 1;
    }
    return n + recursivelySumBetween1AndN( n-1);
}


const double root2Pi = sqrt( 2.0 * 3.141592653589793 );
double normcdf( double x ) {
    if (x<0) {
        return 1-normcdf(-x);
    }
    double k = 1/(1 + 0.2316419*x);
    double poly = k*(0.319381530 + k*(-0.356563782 + k*(1.781477937
                    + k*(-1.821255978 + 1.330274429*k))));
    double approx = 1.0 - 1.0/root2Pi * exp(-0.5*x*x) * poly;
    return approx;
}


double hornerFunction( double x, double a0, double a1) {
    return a0 + x*a1;
}

double hornerFunction( double x, double a0, double a1, double a2) {
    return a0 + x*hornerFunction( x, a1, a2);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3) {
    return a0 + x*hornerFunction( x, a1, a2, a3);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6, double a7) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6, a7);
}

double hornerFunction( double x, double a0, double a1, double a2, double a3, double a4,
                       double a5, double a6, double a7, double a8) {
    return a0 + x*hornerFunction( x, a1, a2, a3, a4, a5, a6, a7, a8);
}

/**
 *  Arguably this is a little easier to read than the original normcdf
 *  function as it makes the use of horner's method obvious.
 */
double normcdfHorner( double x ) {
    if (x<=0) {
        return 1-normcdf(-x);
    }
    double k = 1/(1 + 0.2316419*x);
    double poly = hornerFunction(k,
                                 0.0, 0.319381530, -0.356563782,
                                 1.781477937,-1.821255978,1.330274429);
    double approx = 1.0 - 1.0/root2Pi * exp(-0.5*x*x) * poly;
    return approx;
}

/**
 *   Exercise 7
 */
const double a0 = 2.50662823884;
const double a1 = -18.61500062529;
const double a2 = 41.39119773534;
const double a3 = -25.44106049637;
const double b1 = -8.47351093090;
const double b2 = 23.08336743743;
const double b3 = -21.06224101826;
const double b4 = 3.13082909833;
const double c0 = 0.3374754822726147;
const double c1 = 0.9761690190917186;
const double c2 = 0.1607979714918209;
const double c3 = 0.0276438810333863;
const double c4 = 0.0038405729373609;
const double c5 = 0.0003951896511919;
const double c6 = 0.0000321767881768;
const double c7 = 0.0000002888167364;
const double c8 = 0.0000003960315187;

double norminv( double x ) {
    // We use Moro's algorithm
    double y = x - 0.5;
    if (y<0.42 && y>-0.42) {
        double r = y*y;
        return y*hornerFunction(r,a0,a1,a2,a3)/hornerFunction(r,1.0,b1,b2,b3,b4);
    } else {
        double r;
        if (y<0.0) {
            r = x;
        } else {
            r = 1.0 - x;
        }
        double s = log( -log( r ));
        double t = hornerFunction(s,c0,c1,c2,c3,c4,c5,c6,c7,c8);
        if (x>0.5) {
            return t;
        } else {
            return -t;
        }
    }
}



double blackScholesCallPrice( double strike, double timeToMaturity,
                              double spot, double volatility, double riskFreeRate ) {
    double numerator = log( spot/strike )
                + (riskFreeRate + volatility*volatility*0.5) * timeToMaturity;
    double denominator = volatility * sqrt( timeToMaturity );
    double d1 = numerator/denominator;
    double d2 = d1 - denominator;
    double t1 = normcdf( d1 ) * spot;
    double t2 = normcdf( d2 ) * strike * exp( - riskFreeRate*timeToMaturity);
    return t1 - t2;
}





int main() {

    // Test exercise 1
    cout << "Exercise 1\n";
    cout << "The sum of the numbers between 1 and 100 is ";
    cout << recursivelySumBetween1AndN(100);
    cout << "\n";


    // Test normcdf
    cout << "\nExercise 4\n";
    cout << "normcdf(1.96)=";
    cout << normcdf(1.96)<<"\n";
    cout << "normcdf(-1.96)=";
    cout << normcdf(-1.96)<<"\n";

    // Test normcdfHorner
    cout << "\nExercise 6\n";
    cout << "normcdfHorner(1.96)=";
    cout << normcdfHorner(1.96)<<"\n";
    cout << "normcdfHorner(-1.96)=";
    cout << normcdfHorner(-1.96)<<"\n";

    // Test normimv
    cout << "\nExercise 7\n";
    cout << "This should print out approximately 0, 0.1, 0.2, ...., 1.0\n";
    cout << normcdf( norminv( 0.0 )) << "\n";
    cout << normcdf( norminv( 0.1 )) << "\n";
    cout << normcdf( norminv( 0.2 )) << "\n";
    cout << normcdf( norminv( 0.3 )) << "\n";
    cout << normcdf( norminv( 0.4 )) << "\n";
    cout << normcdf( norminv( 0.5 )) << "\n";
    cout << normcdf( norminv( 0.6 )) << "\n";
    cout << normcdf( norminv( 0.7 )) << "\n";
    cout << normcdf( norminv( 0.8 )) << "\n";
    cout << normcdf( norminv( 0.9 )) << "\n";
    cout << normcdf( norminv( 1.0 )) << "\n";
    cout << "This should print out approximately 0.01, 0.02,0.98,0.99\n";
    cout << normcdf( norminv( 0.01 )) << "\n";
    cout << normcdf( norminv( 0.02)) << "\n";
    cout << normcdf( norminv( 0.98)) << "\n";
    cout << normcdf( norminv( 0.99 )) << "\n";

    // Test option price
    double strike = 100.0;
    double spot = 110.0;
    double vol = 0.1;
    double riskFreeRate = 0.03;
    double timeToMaturity = 0.5;
    cout << "\nExercise8\n";
    cout << "The call price calculated is ";
    cout << blackScholesCallPrice( strike, timeToMaturity, spot, vol, riskFreeRate );
    cout << "\nExpected the value 11.677\n";
    // I simply used an online calculator to get this value

    return 0;
}