Geometric Brownian Motion; Simulation of Stock Price

Question

I made a GBM function in C++ and I believe I am getting too much of a range of stock prices when I start with an initial price of 100 the output can be from [50,400]. I am not sure what I am doing wrong in my code, I am guessing there is something wrong with the way I seed the random standard normal numbers. Please have a look at the function and let me know if there is anything I shold modify or change.

Here is the function:

std::vector<double> GBM(const int M, const int N, const double T, const double r, const double q, const double sigma, const double S0){
    double dt =  T/N;
    std::vector<double> Z;
    std::vector<double> S;
    S.push_back(S0);
    std::mt19937 e2(time(0));
    std::normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        Z.push_back(dist(e2));
    }
    double drift  = exp(dt*((r - q)-0.5*sigma*sigma));
    double vol = sqrt(sigma*sigma*dt);
    for(int i = 1; i < M; i++){
        S.push_back(S[i-1] * drift * exp(vol*Z[i]));
    }
    return S;
}

Here is the main.cpp file that utilizes the function above:

#include <iostream>
#include "LSM.h"
#include <cmath>
#include <ctime>
#include <Eigen/Core>
#include <Eigen/SVD>
#include <iostream>
#include <vector>
#include <random>

std::vector<double> GBM(const int M, const int N, const double T, const double r, const double q, const double sigma, const double S0);

int main(){
    const double r = 0.04;          // Riskless interest rate
    const double q = 0.0;           // Divident yield
    const double sigma = 0.20;      // Volatility of stock
    const double T = 1;             // Time (expiry)
    const int N = 1000;             // Number of time steps
    const double K = 100.0;         // Strike price
    const double S0 = 100.0;        // Initial stock price
    const int M = 10000;                // Number of paths
    const int R = 2;                // Choice of basis for Laguerre polynomial

    //LSM Option_value(r,q,sigma,T,N,K,S0,M,R);

    std::vector<double> s = GBM(M,N,T,r,q,sigma,S0);
    for(int i = 0; i < M; i++){
        std::cout << s[i] << std::endl;
    }



    return 0;
}

A typical output that one should get starting with an initial stock price of 100 is below:

Welcome to StackOverflow. Please read and follow the posting guidelines in the help documentation. [Minimal, complete, verifiable example](http://stackoverflow.com/help/mcve) applies here. We cannot effectively help you until you post your MCVE code and accurately describe the problem. We should be able to paste your posted code into a text file and reproduce the problem you described. — Prune, Oct 23 '17 at 18:19
@Prune I am not familiar with posting on this site, what is MCVE code? — , Oct 24 '17 at 00:07
Follow the link in my comment. Your update is good; just supply a trivial main program to produce this output, and then make it clear what you *wanted* for output, in contrast to what you *did* get. — Prune, Oct 24 '17 at 00:14
From the comments behind the constants, you want to simulate 10000 paths of an integration from 0 to 1 using 1000 subdivision steps, i.e., a step size of `0.001`. What you are doing is integrating one path over 10000 steps of step size 0.001, that is, from 0 to 10. Please correct your code to reflect your initial intentions. — Lutz Lehmann, Oct 24 '17 at 10:09
I rewrite it according to wiki version https://pastebin.com/6bi8zZCs. May it helps. You need to change your M back to 10000 — vadikrobot, Oct 29 '17 at 22:13

Germán · Accepted Answer · 2017-10-31T08:06:05.630

Function GBM should simulate 1 path every time. So no need to supply M. And the path length is, in your code, defined by N instead of M.

If you implement this change, GBM return the whole simulated path. Then you need to call GBM M times in order to calculate all the simulations.

Also there is no need to store all the random numbers generated.

Based on your sample, something like this:

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

// Random generator initialize (only once).
static std::mt19937 rng(time(0)); 

std::vector<double> GBM(const int N, const double T, const double r,
        const double q, const double sigma, const double S0)
{
    double dt =  T/N;
    std::vector<double> S;
    S.push_back(S0);
    std::normal_distribution<double> dist(0.0, 1.0);
    double drift  = exp(dt*((r - q)-0.5*sigma*sigma));
    double vol = sqrt(sigma*sigma*dt);
    for(int i = 1; i < N; i++){
        double Z = dist(rng);
        S.push_back(S[i-1] * drift * exp(vol*Z));
    }
    return S;
}

int main(){
    const double r = 0.04;          // Riskless interest rate
    const double q = 0.0;           // Divident yield
    const double sigma = 0.20;      // Volatility of stock
    const double T = 1;             // Time (expiry)
    const int N = 1000;             // Number of time steps
    const double S0 = 100.0;        // Initial stock price
    const int M = 100;              // Number of paths

    for (int sindx = 0; sindx < M; sindx++)
    {
        std::vector<double> s = GBM(N,T,r,q,sigma,S0);

        std::cout << "Simulation " << sindx << ": "
            << s[0] << ", " << s[1] << " ... " << s[N-2] << ", " << s[N-1]
            << std::endl;
    }

    return 0;
}

Can you provide working C++ code, then I will happily reward the points — , Oct 30 '17 at 20:28

score 1 · Answer 2 · answered Oct 24 '17 at 10:48

1

From the comments behind the constants, you want to simulate 10000 paths of an integration from 0 to 1 using 1000 subdivision steps, i.e., a step size of 0.001.

What you are doing is integrating one path over 10000 steps of step size 0.001, that is, from 0 to 10.

If you do this correctly, the result should look like a list of

S0 * exp( ((r-q)-0.5*sigma*sigma)*T + sigma*sqrt(T)*Z[i] )

as the value of the GBM at time T only depends on W(T) which is distributed as N(0,T) or sqrt(T)*N(0,1).

answered Oct 24 '17 at 10:48

Lutz Lehmann

25,219
2
22
51

Could you provide a bit more information in regards to simulating 10000 paths of an integration from 0 to 1 using 1000 subdivision steps? I may be mistaken but it seems like I need to use two for loops? – Oct 24 '17 at 17:04
Do you need the intermediate values? Usually for something like american options. If only the price distribution at T=1 is of interest, you can just sample as in the second part. -- Else you need nested loops, one for 1..M and the inner for 1..N or vice versa, and a 2-dimensional array to hold the results. – Lutz Lehmann Oct 24 '17 at 17:13
Well just starting with an initial value of 100, I need to generate M stock paths. So it seems I need a matrix instead of just a vector, I was trying to avoid that in order to make things run faster – Oct 24 '17 at 19:50
I made a bounty if you could provide some workable code I would be happy to give you the 100 points – Oct 25 '17 at 22:15

Geometric Brownian Motion; Simulation of Stock Price

2 Answers2