design_pattern_derivative_randomnumber

We have implemented payoff, option, statistics class. In our calculation for option price, we have used gsl library to generate random numbers, gsl is a high quality numerical library, however,  to practise oop in c++, we now try to insulate the random number stream used by a particular simulation from the rest of the program, using c++ design pattern we learned.

There are two major advantages on using a class structure. One is that we can decorate it. For example, suppose we wish to use anti-thetic sampling. We could write a decorator class that does anti-thetic sampling. This can then be combined with any random number generator we have written, and plugged into the Monte Carlo simulator, with no changes to the simulator class. The other advantage is that we might decide not to use pseudo-random (i.e. random) numbers but low-discrepancy numbers instead. Low-discrepancy numbers (sometimes called quasi-random numbers) are sequences of numbers designed to do a good job of filling out space. If we are using a random number class, we could replace this class with a generator for low-discrepancy numbers without changing the interior of our code.

Design considerations. First of all, since we may have many random number generators and we want to be able to add new ones later on without recoding, we use an abstract base class to specify an interface.

class RandomBase {
public:
    RandomBase(unsigned long dimension_):dimension(dimension_){}
    unsigned long getDimension() { return dimension; }
    void setDimension(unsigned long dimension_) { dimension = dimension_; }

    virtual RandomBase* clone() const = 0;
    virtual void getUniforms(std::vector<double>& v) = 0;
    virtual void setSeed(unsigned long Seed) = 0;
    virtual void reSet() = 0;

    virtual void getGaussians(std::vector<double>& v){
        double tmp;
        getUniforms(v);
        for(unsigned long i=0; i<dimension; i++){
            tmp = v[i];
            //std::cout << tmp << std::endl;
            v[i] = InverseCumulativeNormal(tmp);
        }
    }
    virtual ~RandomBase(){}

private:
    unsigned long dimension;
};

We have several member function for this abstract class,

1. set the dimensionality,

2. obtain an array of uniforms of size dimensionality from the generator.

3. states the dimensionality.

4. obtain standard Gaussian draws.

5. reset the generator to its initial state, and to set the seed of the generator.

6. clone method, for deep copy.

In the book, adapter pattern is applied. We now have a random number generator, which cannot fit into our RandomBase framework directly,

#ifndef PARKMILLER_HPP_INCLUDED
#define PARKMILLER_HPP_INCLUDED

const long a = 16807;
const long m = 2147483647;
const long q = 127773;
const long r = 2836;

class ParkMiller {
public:
    ParkMiller(long seed_ = 1):seed(seed_){
        if (seed == 0)
            seed = 1;
    }
    void setSeed(long seed_){
        seed = seed_;
        if (seed == 0)
            seed = 1;
    }
    long getInteger(){
        long k;
        k = seed/q;
        seed = a*(seed-k*q) - r*k;
        if (seed < 0)
            seed += m;
        return seed;
    }

    static unsigned long max(){ return m-1; }
    static unsigned long min(){ return 1; };

private:
    long seed;
};
#endif // PARKMILLER_HPP_INCLUDED

Because ParkMiller class is not an abstract class,  to use adapter pattern, w e then define a new class, which contains an object of ParkMiller class,

class RandomParkMiller : public RandomBase {
public:
    RandomParkMiller(unsigned long dimension_, unsigned long seed_=1):RandomBase(dimension_), inner(seed_), iniSeed(seed_){
        reciprocal = 1 / (1.0 + inner.max());
    }
    virtual RandomBase* clone() const {
        return new RandomParkMiller(*this);
    }
    virtual void getUniforms(std::vector<double>& v){
        v.resize(getDimension());
        for (unsigned long i=0; i < getDimension(); i++){
            double tmp = inner.getInteger()*reciprocal;
            v[i] = tmp;
            //std::cout << v[i] << std::endl;
        }
    }

    virtual void setSeed(unsigned long Seed_){
        iniSeed = Seed_;
        inner.setSeed(Seed_);
    }
    virtual void reSet(){
        inner.setSeed(iniSeed);
    }

private:
    ParkMiller inner;
    unsigned long iniSeed;
    double reciprocal;
};

Then, we make ParkMiller fit into our random number generator framework. In Monte Carlo simulation, one important technique is variance reduction. Below, we use anti-thetic sampling (decorator pattern) to reduce variance. Decorator pattern is similar with adapter pattern above, however, we do not want anti-thetic class to know which kind of random number generator it will decorate, therefore, unlike the example above, we cannot put an object as its data member, we need a pointer. Just as option class has a pointer of payoff type. This time, we use a wrapper class

template< class T>
class Wrapper {
public:
    Wrapper() { DataPtr =0;}

    Wrapper(const T& inner){ DataPtr = inner.clone(); }

    ~Wrapper(){
        if (DataPtr !=0) delete DataPtr;
    }

    Wrapper(const Wrapper<T>& original){
        if (original.DataPtr !=0)
            DataPtr = original.DataPtr->clone();
        else
            DataPtr=0;
    }

    Wrapper& operator=(const Wrapper<T>& original){
        if (this != &original){
            if (DataPtr!=0)
                delete DataPtr;

            DataPtr = (original.DataPtr !=0) ? original.DataPtr->clone() : 0;
        }

        return *this;
    }


    T& operator*()
    {
        return *DataPtr;
    }

    const T& operator*() const
    {
        return *DataPtr;
    }

    const T* const operator->() const{
        return DataPtr;
    }

    T* operator->(){
        return DataPtr;
    }

private:
    T* DataPtr;
};

Based on the wrapper class, we define

class AntiThetic : public RandomBase {
public:
    AntiThetic(const Wrapper<RandomBase>& inner_):RandomBase(*inner_), inner(inner_){
        inner->reSet();
        OddEven = true;
        nextv.resize(getDimension());
    }

    virtual RandomBase* clone() const { return new AntiThetic(*this); }

    virtual void getUniforms(std::vector<double>& v){
        if(OddEven){
            inner->getUniforms(v);
            nextv.resize(v.size());
            for(unsigned long i=0; i<v.size(); i++){
                nextv[i] = 1.0 - v[i];
            }
            OddEven = false;
        }else{
            v = nextv;
            OddEven = true;
        }
    }

    virtual void setSeed(unsigned long Seed_) {
        inner->setSeed(Seed_);
        OddEven = true;
    }

    virtual void reSet(){
        inner->reSet();
        OddEven = true;
    }

private:
    Wrapper<RandomBase> inner;
    bool OddEven;
    std::vector<double> nextv;
};

We use random class in Monte Carlo simulation 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>
#include <./random.hpp>

void SimpleMonteCarlo(const VanillaOption& option, double Expiry, double Spot, double Vol, double r,
                      unsigned long NumberOfPaths, StatisticsMC& gather, RandomBase& generator){

    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;
    double thisGaussian;
    double thisPayOff;

    generator.setDimension(1);
    std::vector<double> rnorm(1);

    for (unsigned long i=0; i < NumberOfPaths; i++){
        generator.getGaussians(rnorm);
        thisGaussian = rnorm[0];
        //std::cout << thisGaussian << std::endl;
        thisSpot = movedSpot*exp(rootVariance*thisGaussian);
        thisPayOff = option.OptionPayOff(thisSpot);
        gather.DumpOneResult(thisPayOff*discount);
    }
}
#endif // MONTECARLO_HPP_INCLUDED

The 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 = 20000;

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

    RandomParkMiller generator(1);
    AntiThetic atgenerator(generator);

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

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

Result is