c++ - Simulating the Geometric Brownian Motion

Question

Background Information:

Consider the Pseudocode:

Question:

What I seem to be having a problem with is the S superscript j part of the Pseudocode of algorithm 1, my professor said that since N = 10,000 and n = 10, I am computing 10 S's. Here is the part of my code that I having trouble with:

double Z[n][N];
    for(int j = 0; j < N; j++){
        for(int i = 0; i < n/2; i++){
            Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
            Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
        }
    }

    /*for(int j = 0; j < N; j++){
        for(int i = 0; i < n; i++){
            cout << Z[j][i] << " ";
        }
        cout << endl;
    }*/

    double S[N][n];
    S[0][0] = S0;
    double t[n+1];
    for(int i = 0; i <= n; i++){
        t[i] = (double)i*T/(n-1);
    }
    for(int j = 1; j <= N; j++){
        for(int i = 1; i <= n; i++){
            S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
        }
    }

To be specific I only require help for computing the S^j part of the code (assuming everything else is correct).

Here is my whole code for completeness:

#include <iostream>
#include <cmath>
#include <math.h>
#include <fstream>
#include <random>

using namespace std;


double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(int n, int N,double u_1, double u_2);
double Box_MullerY(int n, int N,double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);

int main(){
    int n = 10;
    int N = 10000;

    // Calculate BS
    double T = 1.0;
    double K = 50.0;
    double r = 0.1;
    double sigma = 0.3;
    double S0 = 50.0;
    double price = Black_Scholes(T,K,S0,r,sigma);
    //cout << price << endl;

    // Random-shift Halton Sequence
    Random_Shift_Halton(n,N,r,sigma,T, S0);


}

double phi(double x){
    return 0.5 * erfc(-x * M_SQRT1_2);
}

double Black_Scholes(double T, double K, double S0, double r, double sigma){
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
    return C;
}

double Halton_Seq(int index, double base){
    double f = 1.0, r = 0.0;
    while(index > 0){
        f = f/base;
        r = r + f*(fmod(index,base));
        index = index/base;
    }
    return r;
}

double Box_MullerX(int n, int N,double u_1, double u_2){
    double R = -2.0*log(u_1);
    double theta = 2*M_PI*u_2;
    double X = sqrt(R)*cos(theta);

    return X;
}

double Box_MullerY(int n, int N,double u_1, double u_2){
    double R = -2.0*log(u_1);
    double theta = 2*M_PI*u_2;
    double Y = sqrt(R)*sin(theta);

    return Y;
}

double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
    //  Generate the Halton_Sequence
    double hs[N + 1];
    for(int i = 0; i <= N; i++){
        hs[i] = Halton_Seq(i,2.0);
    }
    // Generate two random numbers from U(0,1)
    double u[N+1];
    random_device rd;
    mt19937 e2(rd());
    uniform_real_distribution<double> dist(0, 1);
    for(int i = 0; i <= N; i++){
        u[i] = dist(e2);
    }
    // Add the vector u to each vector hs and take fraction part of the sum
    double shifted_hs[N+1];
    for(int i = 0; i <= N; i++){
        shifted_hs[i] = hs[i] + u[i];
        if(shifted_hs[i] > 1){
            shifted_hs[i] = shifted_hs[i] - 1;
        }
    }

    // Use Geometric Brownian Motion
    /*double Z[N];
    for(int i = 0; i < N; i++){
        if(i % 2 == 0){
            Z[i] = Box_MullerX(n,N,shifted_hs[i+1],shifted_hs[i+2]);
        }else{
            Z[i] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
        }
    }*/
    double Z[n][N];
    for(int j = 0; j < N; j++){
        for(int i = 0; i < n/2; i++){
            Z[2*i][j] = Box_MullerX(n,N,shifted_hs[i],shifted_hs[i+1]);
            Z[2*i+1][j] = Box_MullerY(n,N,shifted_hs[i],shifted_hs[i+1]);
        }
    }

    /*for(int j = 0; j < N; j++){
        for(int i = 0; i < n; i++){
            cout << Z[j][i] << " ";
        }
        cout << endl;
    }*/

    double S[N][n];
    S[0][0] = S0;
    double t[n+1];
    for(int i = 0; i <= n; i++){
        t[i] = (double)i*T/(n-1);
    }
    for(int j = 1; j <= N; j++){
        for(int i = 1; i <= n; i++){
            S[j][i] = S[j-1][i-1]*exp((mu - sigma/2)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j][i]);
        }
    }

    for(int j = 1; j <= N; j++){
        for(int i = 1; i <= n; i++){
            cout << S[j][i] << endl;
        }
    }



    // Use random-shift halton sequence to obtain 40 independent estimates for the price of the option

}

Answer 1

Please explain what are you expect and what are you get. That helps to understand possible variants.

And i noticed that you are mixing integer and double values in the code. That is not the good practice. Also - Mark Setchell right - check your sigma value. I suppose it should be squared. And try to correctly set the initial values for S.

I tried to fix the errors in your code the way I understood it. The code changed a little because of the older version of the compiler.

#include <iostream>
#include <cmath>
#include <math.h>
#include <random>

using namespace std;

const double M_PI = 3.1415926535897932384626433832795;
const double M_SQRT1_2 = 0.70710678118654752440084436210485;

double phi(long double x);
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double Halton_Seq(int index, double base);
double Box_MullerX(double u_1, double u_2);
double Box_MullerY(double u_1, double u_2);
double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0);


double erfc(double v) {
    double y = 1.0 / ( 1.0 + 0.3275911 * (v >= 0 ? v : -v));
    y = 1 - (((((
        + 1.061405429  * y
        - 1.453152027) * y
        + 1.421413741) * y
        - 0.284496736) * y
        + 0.254829592) * y)
        * exp (-v * v);
    return (v >= 0 ? y : -y);
}

double phi(double x){
    return 0.5 * erfc(-x * M_SQRT1_2);
}

double Black_Scholes(double T, double K, double S0, double r, double sigma){
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
    return C;
}

double Halton_Seq(int index, double base){
    double f = 1.0, r = 0.0;
    double indice = (double)index;
    while(indice >= base){
        f = f/base;
        r = r + f*(fmod(indice,base));
        indice = ceil(indice/base);
    }
    return r;
}

double Box_MullerX(double u_1, double u_2){
    double R = -2.0*log(u_1);
    double theta = 2.0*M_PI*u_2;
    double X = sqrt(R)*cos(theta);

    return X;
}

double Box_MullerY(double u_1, double u_2){
    double R = -2.0*log(u_1);
    double theta = 2.0*M_PI*u_2;
    double Y = sqrt(R)*sin(theta);

    return Y;
}

double Random_Shift_Halton(int n,int N,double mu, double sigma, double T, double S0){
    //  Generate the Halton_Sequence
    double* hs = new double[N+1];
    for(int i = 0; i <= N; i++){
        hs[i] = Halton_Seq(i,2.0);
    }
    // Generate two random numbers from U(0,1)
    double* u = new double[N+1];
    random_device rd;
    mt19937 e2(rd());
    uniform_real_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i <= N; i++){
        u[i] = dist(e2);
    }
    // Add the vector u to each vector hs and take fraction part of the sum
    double* shifted_hs = new double[N+1];
    for(int i = 0; i <= N; i++){
        shifted_hs[i] = hs[i] + u[i];
        if(shifted_hs[i] > 1){
            shifted_hs[i] = shifted_hs[i] - 1.0;
        }
    }

    // Use Geometric Brownian Motion
    double** Z = new double*[N];    
    for(int i = 0; i < N; ++i)
       Z[i] = new double[n];
    for(int j = 0; j < N; j++){
        for(int i = 0; i < n/2; i++){
            Z[j][2*i] = Box_MullerX(shifted_hs[j],shifted_hs[j+1]);
            Z[j][2*i+1] = Box_MullerY(shifted_hs[j],shifted_hs[j+1]);
        }
    }

    double** S = new double*[N+1];  
    for(int i = 0; i <= N; ++i) 
       S[i] = new double[n+1];
    for (int i = 0; i <= N; i++) S[i][0] = S0;
    for (int j = 0; j <= n; j++) S[0][j] = S0;
    double* t = new double[n+1];
    for(int i = 0; i <= n; i++){
        t[i] = (double)i*T/(double)(n-1);
    }
    for(int j = 1; j <= N; j++){
        for(int i = 1; i <= n; i++){
            S[j][i] = S[j-1][i-1]*exp((mu - sigma/2.0)*(t[i] - t[i-1]) + sqrt(sigma)*sqrt(t[i] - t[i-1])*Z[j-1][i-1]);
        }
    }

    // Use random-shift halton sequence to obtain 40 independent estimates for the price of the option
    //....

    delete(hs);
    delete(u);
    delete(shifted_hs);
    for(int i = 0; i < N; ++i) {
       delete(Z[i]);
    }
    for(int i = 0; i <= N; ++i) delete(S[i]);
    delete(S);
    delete(Z);
    delete(t);

    return 0.0;
}

int main()
{
    int n = 10;
    int N = 10000;

    // Calculate BS
    double T = 1.0;
    double K = 50.0;
    double r = 0.1;
    double sigma = 0.3;
    sigma *= sigma;
    double S0 = 50.0;
    double price = Black_Scholes(T,K,S0,r,sigma);
    //cout << price << endl;

    // Random-shift Halton Sequence
    Random_Shift_Halton(n,N,r,sigma,T, S0);
    return 0;
}

Is that what are you looking for?