Monte Carlo Asian Option Pricing OOP Framework

Box-Muller  Algorithm to Generate Normal Random Numbers

double gaussian_box_muller(){
    double x = 0.0;
    double y = 0.0;
    double euclid_sq = 0.0;

    do{
        x = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
        y = 2.0 * rand() / static_cast<double>(RAND_MAX)-1;
        euclid_sq = x*x + y*y;
    } while (euclid_sq >= 1.0);

    return x*sqrt(-2*log(euclid_sq)/euclid_sq);
}

But, since C++11, we do not need to write basic distribution random number generator, head file <random> provides us with certain distributions.

Basic geometric Brownian motion simulation, note that we need to pass seed to the random number generator.

/* geometric brownian motion */
void spot_price_bm(std::vector<double>& spot_prices,
                   const double& r,
                   const double& v,
                   const double& T,
                   const int& seed){

    double dt =  T/static_cast<double>(spot_prices.size());
    double drift = exp((r-0.5*v*v)*dt);
    double vol = sqrt(v*v*dt);

    std::default_random_engine generator(seed);
    std::normal_distribution<double> normal(0.0, 1.0);
    for (unsigned i=1; i<spot_prices.size(); i++){
        double rvn = normal(generator);
        spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvn);
    }
}

Jump-diffusion process simulation.

/* jump diffusion process, with normal distribution */
void spot_price_bmj(std::vector<double>& spot_prices,
                   const double& r,
                   const double& v,
                   const double& T,
                   const double& lambda,
                   const double& muj,
                   const double& sigmaj,
                   const int& seed){
    double dt =  T/static_cast<double>(spot_prices.size());
    double drift = exp((r-0.5*v*v-lambda*(muj+0.5*sigmaj*sigmaj-1))*dt);
    double vol = sqrt(v*v*dt);

    std::default_random_engine generator(seed);
    std::normal_distribution<double> rvnbm(0.0, 1.0);
    std::normal_distribution<double> rvnjump(muj, sigmaj);
    std::poisson_distribution<int> rvp(lambda);

    for (unsigned i=1; i<spot_prices.size(); i++){
        double rvbm = rvnbm(generator);
        int rnp = rvp(generator);
        double sumj = 0;

        if (rnp==0){
            spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvbm);
        } else {
            for (int i=0; i<rnp; i++){
                sumj += rvnjump(generator);
            }
            spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvbm)*exp(sumj);
        }
    }
}

The payoff classes

#ifndef PAYOFF_HPP_INCLUDED
#define PAYOFF_HPP_INCLUDED

class PayOff{
public:
    PayOff(){};
    virtual ~PayOff(){};
    virtual double operator()(const double& S) const = 0;
};

class PayOffCall : public PayOff{
private:
    double K;
public:
    PayOffCall(const double& K_):K(K_){};
    virtual ~PayOffCall() {};
    virtual double operator()(const double& S) const { return std::max(S-K, 0.0); }
};

class PayOffPut : public PayOff{
private:
    double K;
public:
    PayOffPut(const double& K_):K(K_){};
    virtual ~PayOffPut() {};
    virtual double operator()(const double& S) const { return std::max(K-S, 0.0); }
};

#endif // PAYOFF_HPP_INCLUDED

The option class, option class has a pointer to payoff classes

#ifndef ASIANOPTION_HPP_INCLUDED
#define ASIANOPTION_HPP_INCLUDED

# include "payoff.hpp"
# include "path_generator.hpp"

class AsianOption{
protected:
    PayOff* payoff;
public:
    AsianOption(PayOff* payoff_):payoff(payoff_){};
    virtual ~AsianOption() {};
    virtual double pay_off_price(const std::vector<double>& spot_prices) const = 0;
};

class AsianOptionArithmetic : public AsianOption {
public:
    AsianOptionArithmetic(PayOff* payoff_):AsianOption(payoff_){};
    virtual ~AsianOptionArithmetic() {};

    virtual double pay_off_price(const std::vector<double>& spot_prices) const {
        double sum = 0.0, ave = 0.0;
        for (unsigned i=0; i<spot_prices.size(); i++){
            sum += spot_prices[i];
        }
        ave = sum / static_cast<double>(spot_prices.size());
        //std::cout << ave << '\t' << (*payoff)(ave) << std::endl;
        return (*payoff)(ave);
    }
};

class AsianOptionGeometric : public AsianOption {
public:
    AsianOptionGeometric(PayOff* payoff_):AsianOption(payoff_){};
    virtual ~AsianOptionGeometric() {};

    virtual double pay_off_price(const std::vector<double>& spot_prices) const {
        double sum = 0.0, ave = 0.0;
        for (unsigned i=0; i<spot_prices.size(); i++){
            sum += log(spot_prices[i]);
        }
        ave = exp(sum / static_cast<double>(spot_prices.size()));
        return (*payoff)(ave);
    }
};

#endif // ASIANOPTION_HPP_INCLUDED

Main funciton

#include <iostream>
#include <random>
#include <vector>

#include "payoff.hpp"
#include "AsianOption.hpp"
#include "path_generator.hpp"

using namespace std;

int main(){
    // set parameters for simulation
    int npaths = 10000;
    int nsteps = 100;

    // set parameters for model
    double S = 30.0;
    double K = 30.0;
    double r = 0.01;
    double v = 0.35;
    double T = 1.0;
    double df =  exp(-r*T);

    std::vector<double> spot_price(nsteps, S);
    PayOff* call_payoff = new PayOffCall(K);
    AsianOptionArithmetic asiana(call_payoff);

    double payoffsum = 0;
    for(int i=0; i<npaths; i++){
        spot_price_bm(spot_price, r, v, T, i);
        payoffsum += asiana.pay_off_price(spot_price);
    }

    double optionprice = payoffsum*df/npaths;
    cout << "Option Price: " << optionprice << endl;

    delete call_payoff;

    return 0;
}

path_generator contains a set of functions for stochastic process simulation

#ifndef PATH_GENERATOR_HPP_INCLUDED
#define PATH_GENERATOR_HPP_INCLUDED

#include <vector>
#include <cmath>
#include <ctime>
#include <random>

/* geometric brownian motion */
void spot_price_bm(std::vector<double>& spot_prices,
                   const double& r,
                   const double& v,
                   const double& T,
                   const int& seed){

    double dt =  T/static_cast<double>(spot_prices.size());
    double drift = exp((r-0.5*v*v)*dt);
    double vol = sqrt(v*v*dt);

    std::default_random_engine generator(seed);
    std::normal_distribution<double> normal(0.0, 1.0);
    for (unsigned i=1; i<spot_prices.size(); i++){
        double rvn = normal(generator);
        spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvn);
    }
}

/* jump diffusion process, with normal distribution */
void spot_price_bmj(std::vector<double>& spot_prices,
                   const double& r,
                   const double& v,
                   const double& T,
                   const double& lambda,
                   const double& muj,
                   const double& sigmaj,
                   const int& seed){
    double dt =  T/static_cast<double>(spot_prices.size());
    double drift = exp((r-0.5*v*v-lambda*(muj+0.5*sigmaj*sigmaj-1))*dt);
    double vol = sqrt(v*v*dt);

    std::default_random_engine generator(seed);
    std::normal_distribution<double> rvnbm(0.0, 1.0);
    std::normal_distribution<double> rvnjump(muj, sigmaj);
    std::poisson_distribution<int> rvp(lambda);

    for (unsigned i=1; i<spot_prices.size(); i++){
        double rvbm = rvnbm(generator);
        int rnp = rvp(generator);
        double sumj = 0;

        if (rnp==0){
            spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvbm);
        } else {
            for (int i=0; i<rnp; i++){
                sumj += rvnjump(generator);
            }
            spot_prices[i] = spot_prices[i-1]*drift*exp(vol*rvbm)*exp(sumj);
        }
    }
}

#endif // PATH_GENERATOR_HPP_INCLUDED