design_pattern_derivative_statistics

Our monte carlo rounte  lacks any indication of the simulation’s convergence. We could make the routine return standard error, or a convergence table, or simply have it return the value for every single path and analyze the results elsewhere. We make the statistics gatherer an class and an input for Monte Carlo simulation function.

class StatisticsMC {
public:
    StatisticsMC(){}
    virtual void DumpOneResult(double result) = 0;
    virtual std::vector<std::vector<double>> GetResultsSoFar() const = 0;
    virtual StatisticsMC* clone() const = 0;
    virtual ~StatisticsMC(){}
};

1. DumpOneResult save a result into the internal data container.

2. GetResultSoFar returns statistics so far, as the simulation continues, let us observe the convergence.

3. Clone method garantees the deep copy.

Suppose we want to see the mean statistics,  we could extend  the base class

class StatisticsMean : public StatisticsMC {
public:
    StatisticsMean() : RunningSum(0.0), PathsDone(0){}
    virtual void DumpOneResult(double result){
        PathsDone++;
        RunningSum += result;
    }
    virtual std::vector<std::vector<double>> GetResultsSoFar() const {
        std::vector<std::vector<double>> Results(1);
        Results[0].resize(1);
        Results[0][0] = RunningSum/PathsDone;

        return Results;
    }
    virtual StatisticsMC* clone() const{ return new StatisticsMean(*this); }
private:
    double RunningSum;
    unsigned long PathsDone;
};

If we use a statistics gatherer and run the simulation it will tell us the relevant statistics for the entire simulation. However, it does not necessarily give us a feel for how well the simulation has converged. One standard method of checking the convergence is to examine the standard error of the simulation; that is measure the sample standard deviation and divide by the square root of the number of paths. We could  use decorator pattern here, which means

1. our class inherits from StasticsMC.

2. our class contains a pointer for StatisticsMC.

3. We need to implement ‘rule of three’， since our class has a pointer.

class ConvergenceTable : public StatisticsMC {
public:
    ConvergenceTable(const StatisticsMC& inner_):StoppingPoint(2), PathsDone(0){
        inner = inner_.clone();
    }
    ConvergenceTable(const ConvergenceTable& original){
        inner = original.inner->clone();
        StoppingPoint = original.StoppingPoint;
        PathsDone = original.PathsDone;
        ResultsSoFar = original.ResultsSoFar;
    }
    ConvergenceTable& operator=(const ConvergenceTable& original){
        if(this != &original){
            delete inner;
            inner = original.inner->clone();
            StoppingPoint = original.StoppingPoint;
            PathsDone = original.PathsDone;
            ResultsSoFar = original.ResultsSoFar;
        }
        return *this;
    }
    virtual ~ConvergenceTable(){
    if (inner != 0)
        delete inner;
    }

    virtual StatisticsMC* clone() const { return new ConvergenceTable(*this); }
    virtual void DumpOneResult(double result);
    virtual std::vector<std::vector<double>> GetResultsSoFar() const;

private:
    StatisticsMC* inner;
    std::vector<std::vector<double>> ResultsSoFar;
    unsigned long StoppingPoint;
    unsigned long PathsDone;
};

void ConvergenceTable::DumpOneResult(double result){
    inner->DumpOneResult(result);
    ++PathsDone;

    if (PathsDone == StoppingPoint){
        StoppingPoint *= 2;
        std::vector<std::vector<double>> thisResult(inner->GetResultsSoFar());
        for(unsigned long i=0; i<thisResult.size(); i++){
            thisResult[i].push_back(PathsDone);
            ResultsSoFar.push_back(thisResult[i]);
        }
    }
}

std::vector<std::vector<double>> ConvergenceTable::GetResultsSoFar() const {
    std::vector<std::vector<double>> tmp(ResultsSoFar);
    if (PathsDone*2 != StoppingPoint){
        std::vector<std::vector<double>> thisResult(inner->GetResultsSoFar());
        for(unsigned long i=0; i<thisResult.size(); i++){
            thisResult[i].push_back(PathsDone);
            tmp.push_back(thisResult[i]);
        }
    }
    return tmp;
}

The Monte carlo simuation function

#ifndef MONTECARLO_HPP_INCLUDED
#define MONTECARLO_HPP_INCLUDED

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <./payoff.hpp>
#include <./vanillaoption.hpp>
#include <./statisticsmc.hpp>

void SimpleMonteCarlo(const VanillaOption& option, double Expiry, double Spot, double Vol,
                      double r, unsigned long NumberOfPaths, StatisticsMC& gather){
    const gsl_rng_type * T;
    gsl_rng * rnd;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    rnd = gsl_rng_alloc (T);

    double variance = Vol*Vol*Expiry;
    double rootVariance = sqrt(variance);
    double itoCorrection = -0.5*variance;
    double movedSpot = Spot*exp(r*Expiry+itoCorrection);
    double discount = exp(-r*Expiry);
    double thisSpot;

    for (unsigned long i=0; i < NumberOfPaths; i++){
        double thisGaussian = gsl_ran_gaussian(rnd, 1.0);
        thisSpot = movedSpot*exp(rootVariance*thisGaussian);
        double thisPayOff = option.OptionPayOff(thisSpot);
        gather.DumpOneResult(thisPayOff*discount);
    }

    gsl_rng_free (rnd);
}
#endif // MONTECARLO_HPP_INCLUDED

We test the payoff, option and statisticsmc classes in main function

#include <iostream>
#include <cmath>
#include <./montecarlo.hpp>
#include <./payoff.hpp>
#include <./statisticsmc.hpp>
#include <./vanillaoption.hpp>

int main(){
    // set parameters
    double spot = 15.0;
    double expiry = 1.0;
    double strike = 15.0;
    double volatility = 0.4;
    double r = 0.01;
    unsigned long NumberPaths = 200000;

    PayOffCall payoffcall(strike);
    VanillaOption option(payoffcall, expiry);
    StatisticsMean gather;
    ConvergenceTable ct(gather);

    // calculate option price and generate statistics
    // SimpleMonteCarlo(option, expiry, spot, volatility, r, NumberPaths, gather);
    SimpleMonteCarlo(option, expiry, spot, volatility, r, NumberPaths, ct);

    // print results
    std::vector<std::vector<double>> results = ct.GetResultsSoFar();
    std::cout << "For Call Option Price: " << std::endl;
    for (unsigned long i=0; i<results.size(); i++){
        for (unsigned long j=0; j<results[i].size(); j++){
            std::cout << results[i][j] << " ";
        }
        std::cout << std::endl;
    }

    return 0;
}

the result is