Simulating stock prices using GBM
- Thread starter earth_imperator
- Start date
Here's a second test/proof method:
Testing the quality of the GBM algorithm used
Example for S0=100, t=1 year, volatility s=0.3 (ie. 30%), drift u=0% :
z=-1 ; m1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 74.08182206817179
z=+1 ; p1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 134.9858807576003
Perform a Monte Carlo simulation with say 10 million runs by using 1 interval (N=1),
and count the generated St that fall within the range of -1SD thru +1SD (above the variables m1SD and p1SD).
Set this in relation to the total of the generated numbers.
It has to give the well known probability for the range from -1SD to +1SD: 68.27% --> cf. the "68–95–99.7 rule".
If it doesn't give on average about 68.27%, than the used GBM algorithm is buggy.
FYI: as already proven, the standard GBM algoritm is indeed buggy as it gives only about 67.7% instead of 68.2%
And, as someone said:
"Technically, z-scoring does not depend on any distributional assumptions, such as normality. It's just a way of describing how far observations are from the mean, no matter what the distribution happens to be. So no harm in z-scoring non-normal variables."
Ie. the use of such z-scoring is distribution agnostic, ie. it works for any kind of distribution, incl. lognormal distribution as in stock price case.
Testing the quality of the GBM algorithm used
Example for S0=100, t=1 year, volatility s=0.3 (ie. 30%), drift u=0% :
z=-1 ; m1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 74.08182206817179
z=+1 ; p1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 134.9858807576003
Perform a Monte Carlo simulation with say 10 million runs by using 1 interval (N=1),
and count the generated St that fall within the range of -1SD thru +1SD (above the variables m1SD and p1SD).
Set this in relation to the total of the generated numbers.
It has to give the well known probability for the range from -1SD to +1SD: 68.27% --> cf. the "68–95–99.7 rule".
If it doesn't give on average about 68.27%, than the used GBM algorithm is buggy.
FYI: as already proven, the standard GBM algoritm is indeed buggy as it gives only about 67.7% instead of 68.2%
And, as someone said:
"Technically, z-scoring does not depend on any distributional assumptions, such as normality. It's just a way of describing how far observations are from the mean, no matter what the distribution happens to be. So no harm in z-scoring non-normal variables."
Ie. the use of such z-scoring is distribution agnostic, ie. it works for any kind of distribution, incl. lognormal distribution as in stock price case.
Last edited:
More insight: look how fucking buggy the standard GBM algorithm really is: it degrades even much worser for higher volatilities! Instead of 68.27% it generates such results:Here's a second test/proof method:
Testing the quality of the GBM algorithm used
Example for S0=100, t=1 year, volatility s=0.3 (ie. 30%), drift u=0% :
z=-1 ; m1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 74.08182206817179
z=+1 ; p1SD = S0 * exp(u * t + s * sqrt(t) * z) gives 134.9858807576003
Perform a Monte Carlo simulation with say 10 million runs by using 1 interval (N=1),
and count the generated St that fall within the range of -1SD thru +1SD (above the variables m1SD and p1SD).
Set this in relation to the total of the generated numbers.
It has to give the well known probability for the range from -1SD to +1SD: 68.27% --> cf. the "68–95–99.7 rule".
If it doesn't give on average about 68.27%, than the used GBM algorithm is buggy.
FYI: as already proven, the standard GBM algoritm is indeed buggy as it gives only about 67.7% instead of 68.2%
And, as someone said:
"Technically, z-scoring does not depend on any distributional assumptions, such as normality. It's just a way of describing how far observations are from the mean, no matter what the distribution happens to be. So no harm in z-scoring non-normal variables."
Ie. the use of such z-scoring is distribution agnostic, ie. it works for any kind of distribution, incl. lognormal distribution as in stock price case.
Code:
Standard GBM algorithm :
S0=100 u=0 t=1 :
s=0.3 : 67.7% (inside -1SD to +1SD)
s=0.5 : 66.8%
s=0.75 : 64.9%
s=1.0 : 62.4%
s=1.25 : 59.4%
s=1.5 : 55.8%And imagine this: this buggy algorithm gets used also in the BSM option pricing model
, which even was awarded the Nobel Prize... 'diots they are!
Last edited:
Has someone a rationale explanation for this related problem ? :
As demonstrated, simulating stock prices with N=1 (ie. just 1 bar in t) should draw in about 68.27% of the cases prices that land within 1SD around the initial price.
Now imagine this: what if say N=3, ie. 3 bars in t: what hit rate should this have? Should it be the same hit rate of 68.27% ?
Anybody know?
Since the outcome is "conic" like a parabola (ie. it expands/grows the more the further t is, cf. pic in OP), then it seems IMO logical that the end result should be the same, ie. 68.27%, but simulations give a higher hit rate... Why?
As demonstrated, simulating stock prices with N=1 (ie. just 1 bar in t) should draw in about 68.27% of the cases prices that land within 1SD around the initial price.
Now imagine this: what if say N=3, ie. 3 bars in t: what hit rate should this have? Should it be the same hit rate of 68.27% ?
Anybody know?
Since the outcome is "conic" like a parabola (ie. it expands/grows the more the further t is, cf. pic in OP), then it seems IMO logical that the end result should be the same, ie. 68.27%, but simulations give a higher hit rate... Why?
Last edited:
I just learned that there are even some people out there, especially in Asia, who really believe that they can predict the future stock price path by simulating it with GBM ! There are even some academic papers doing exactly this! Unbelievable!
Here's one example from an Australian university (a PDF research paper):
https://ro.uow.edu.au/cgi/viewcontent.cgi?article=1705&context=aabfj
"
AABFJ | Volume 10, no. 3, 2016
Simulating Stock Prices Using Geometric Brownian Motion: Evidence from Australian Companies
Krishna Reddy 1 and Vaughan Clinton 2
Abstract
This study uses the geometric Brownian motion (GBM) method to simulate stock price paths, and tests whether the simulated stock prices align with actual stock returns. The sample for this study was based on the large listed Australian companies listed on the S&P/ASX 50 Index. Daily stock price data was obtained from the Thomson One database over the period 1 January 2013 to 31 December 2014. The findings are slightly encouraging as results show that over all time horizons the chances of a stock price simulated using GBM moving in the same direction as real stock prices was a little greater than 50 percent. However, the results
improved slightly when portfolios were formed.
[...]"
What a BS!
Here's one example from an Australian university (a PDF research paper):
https://ro.uow.edu.au/cgi/viewcontent.cgi?article=1705&context=aabfj
"
AABFJ | Volume 10, no. 3, 2016
Simulating Stock Prices Using Geometric Brownian Motion: Evidence from Australian Companies
Krishna Reddy 1 and Vaughan Clinton 2
Abstract
This study uses the geometric Brownian motion (GBM) method to simulate stock price paths, and tests whether the simulated stock prices align with actual stock returns. The sample for this study was based on the large listed Australian companies listed on the S&P/ASX 50 Index. Daily stock price data was obtained from the Thomson One database over the period 1 January 2013 to 31 December 2014. The findings are slightly encouraging as results show that over all time horizons the chances of a stock price simulated using GBM moving in the same direction as real stock prices was a little greater than 50 percent. However, the results
improved slightly when portfolios were formed.
[...]"
What a BS!
This problem already has been solved: it is logical that the hit rate (in %) is higher the higher N is. I even was able to develop a formula for it. Simulations verified its correctness.
For N=2 the value is 76.27%
For N=3 the value is 79.29%
For N=10 the value is 83.33%
For N=100 the value is 84.77%
For N=1000 the value is 84.92%
etc.
Last edited:
Here's C++ test code for testing QuantLib's GBM.
It practically demonstrates and proves that QuantLib's GBM is indeed very buggy:
It practically demonstrates and proves that QuantLib's GBM is indeed very buggy:
Code:
/*
Testing_GBM_of_QuantLib.cpp
Statistically testing QuantLib's GBM implementation via a simple Monte Carlo simulation loop.
This is a minimalistic standalone test solution w/o the Boost test framework.
We want to test/verify whether the GBM in QuantLib passes the "68–95–99.7 rule test" --> https://en.wikipedia.org/wiki/68–95–99.7_rule
It has to give the well known probability distribution of 68.27% for the range from -1SD to +1SD
If it doesn't give on average about 68.27%, than the used GBM algorithm is buggy.
2023-09-08-Fr: init by @earth_imperator
Compile:
g++ -Wall -Wextra -O2 -std=c++11 -o Testing_GBM_of_QuantLib.exe Testing_GBM_of_QuantLib.cpp -lQuantLib
Run:
./Testing_GBM_of_QuantLib.exe
See main() for optional program parameters.
See also:
https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548
https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5851336
https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861025 (timeSteps > 1)
https://www.elitetrader.com/et/threads/any-quantlib-experts-here.376108/page-2#post-5861021
Expected for 1SD:
timeSteps HitRate
1 68.27% see https://en.wikipedia.org/wiki/68–95–99.7_rule
2 76.27% see https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861025
3 79.29%
4 80.79%
5 81.66%
10 83.33%
100 84.77%
1000 84.92%
Results of QuantLib testing:
Since in QuantLib it's practically impossible to set timeSteps=1 (b/c then it would generate always just the same initial value),
then we use timeSteps=2 by default. Then the result must 76.27% (cf. table above). Below we tested many timeSteps:
timeSteps HitRate
1 Not possible to test
2 92.02%
3 89.39%
4 88.12%
5 87.47%
10 86.10%
100 84.66%
1000 85.59% (check)
These results show that QuantLib's GBM is definitely buggy. Q.E.D.
*/
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <string>
#include <ql/quantlib.hpp>
using namespace std;
using namespace QuantLib;
void Testing_GBM_of_QuantLib(const size_t timeSteps, const size_t nRuns, size_t seed,
const double S, const double mu, const double sigma, const double t)
{
while (!seed) seed = time(0);
const size_t nGenAll = nRuns * timeSteps;
const double m1SD = S * exp(mu * t + sigma * sqrt(t) * -1.0); // Sx at -1SD
const double p1SD = S * exp(mu * t + sigma * sqrt(t) * 1.0); // Sx at +1SD
printf("timeSteps=%zu nRuns=%zu seed=%zu S=%lf mu=%lf sigma=%lf t=%lf : nGenAll=%zu : -1SD=%lf +1SD=%lf\n",
timeSteps, nRuns, seed, S, mu, sigma, t, nGenAll, m1SD, p1SD);
// counting all generated values that fall within 1SD around the initial price
size_t cHit = 0;
// Monte Carlo simulation loop
for (size_t i = 0; i < nRuns; ++i)
{
// instantiate Geometric Brownian Motion (GBM) stochastic process
const auto& gbm = boost::shared_ptr<StochasticProcess>(new GeometricBrownianMotionProcess(S, mu, sigma));
// generate sequence of normally distributed random numbers from uniform distribution using Box-Muller transformation
if (i) seed = rand(); // for seed we use rand() for being able to reproduce the whole sequence (for DBG etc)...
MersenneTwisterUniformRng rng(seed);
//
typedef BoxMullerGaussianRng<MersenneTwisterUniformRng> MersenneBoxMuller;
MersenneBoxMuller boxMullerRng(rng);
//
RandomSequenceGenerator<MersenneBoxMuller> gsg(timeSteps, boxMullerRng);
// generate simulated path of stock price
PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> gbmPathGenerator(gbm, t, timeSteps, gsg, false);
const Path& samplePath = gbmPathGenerator.next().value;
// count the hits within 1SD around the initial price
for (size_t j = 0; j < timeSteps; ++j)
{
const auto Sx = samplePath.at(j);
if (Sx >= m1SD && Sx <= p1SD) ++cHit;
}
}
// result:
const double pct = double(cHit) / nGenAll * 100;
printf("cHit=%zu/%zu(%.3lf%%)\n", cHit, nGenAll, pct);
}
int main(int argc, char* argv[])
{
const size_t timeSteps = argc > 1 ? atol(argv[1]) : 2; // must be >= 2 as the first generated Sx equals initial price
size_t nRuns = argc > 2 ? atol(argv[2]) : 0; // 0 means to take nRuns = 1000000 / timeSteps
if (!nRuns) nRuns = 1000000 / timeSteps;
const size_t seed = argc > 3 ? atol(argv[3]) : 0; // initial seed for the RNG. 0 means random initital seed
const double S = argc > 4 ? atof(argv[4]) : 100.0; // initial price
const double mu = argc > 5 ? atof(argv[5]) : 0.0; // drift
const double sigma = argc > 6 ? atof(argv[6]) : 0.30; // stddev, ie. volatility / 100
const double t = argc > 7 ? atof(argv[7]) : 1.0; // time in yrs
Testing_GBM_of_QuantLib(timeSteps, nRuns, seed, S, mu, sigma, t);
return 0;
}
Last edited:
Completing/extending the above table:
Expected HitRate% for 1SD around initial stock price for various GBM timeSteps:
Code:
Expected HitRate% for 1SD around initial stock price for various GBM timeSteps
This table is an extension of https://en.wikipedia.org/wiki/68–95–99.7_rule as there only timeStep=1 is given.
Params: S=100 rPct=0 qPct=0 IV=30(s=0.3) DIY=365.00 DTE=365(t=1 dt=1 dd=365) zFm=-1(SxFm=74.081822) zTo=+1(SxTo=134.985881)
These numbers are exact values, not the result of simulations.
timeSteps Expected_HitRate
1 68.2689% cf. https://en.wikipedia.org/wiki/68–95–99.7_rule
2 76.2695%
3 79.2918%
4 80.7919%
5 81.6648%
6 82.2304%
7 82.6265%
8 82.9197%
9 83.1461%
10 83.3265%
15 83.8653%
20 84.1337%
25 84.2942%
30 84.4010%
35 84.4771%
40 84.5341%
45 84.5785%
50 84.6139%
60 84.6671%
70 84.7050%
80 84.7334%
90 84.7555%
100 84.7732%
500 84.9003%
1000 84.9162%
10000 84.9305%
100000 84.9319%
1000000 84.9320%
10000000 84.9320%
Last edited:
Here*s a new bufixed version (v1.1).
Now QuantLib passes the said test.
Now also timeSteps=1 is made possible too.
Code:
/*
Testing_GBM_of_QuantLib.cpp
v1.1
Statistically testing QuantLib's GBM implementation via a simple Monte Carlo simulation loop.
This is a minimalistic standalone test solution w/o the Boost test framework.
We want to test/verify whether the GBM in QuantLib passes the "68–95–99.7 rule test" --> https://en.wikipedia.org/wiki/68–95–99.7_rule
It has to give the well known probability distribution of 68.27% for the range from -1SD to +1SD
If it doesn't give on average about 68.27%, than the used GBM algorithm is buggy.
History:
2023-09-09-Sa: v1.1 Fixed a bug in accessing the generated path data: now skipping 0th element and still taking timeSteps elements
2023-09-08-Fr: v1.0 Init by @earth_imperator: posted to https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368
Compile:
g++ -Wall -Wextra -O2 -std=c++11 -o Testing_GBM_of_QuantLib.exe Testing_GBM_of_QuantLib.cpp -lQuantLib
Run:
./Testing_GBM_of_QuantLib.exe
See main() for optional program parameters.
See also:
[1] https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548
{2] https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5851336
[3[ https://en.wikipedia.org/wiki/68–95–99.7_rule
[4] https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861530 (timeSteps > 1)
Expected for 1SD:
timeSteps HitRate
1 68.27% cf. [3]
2 76.27% cf. [4]
3 79.29%
4 80.79%
5 81.66%
10 83.33%
100 84.77%
1000 84.92%
Results of QuantLib testing:
timeSteps HitRate
1 68.46%
2 75.59%
3 78.42%
4 79.93%
5 80.94%
10 82.85%
100 84.33%
1000 85.55%
These new test results, after fixing the test program in v1.1, show that
QuantLib's GBM using the BoxMullerGaussianRng<MersenneTwisterUniformRng> is OK.
It seems that Ito's lemma does not get used, b/c otherwise the result would differ much.
*/
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <ql/quantlib.hpp>
using namespace std;
using namespace QuantLib;
void Testing_GBM_of_QuantLib(const size_t timeSteps, const size_t nRuns, size_t seed,
const double S, const double mu, const double sigma, const double t)
{
while (!seed) seed = time(0);
const size_t nGenAll = nRuns * timeSteps;
const double m1SD = S * exp(mu * t + sigma * sqrt(t) * -1.0); // Sx at -1SD
const double p1SD = S * exp(mu * t + sigma * sqrt(t) * 1.0); // Sx at +1SD
printf("timeSteps=%zu nRuns=%zu seed=%zu S=%lf mu=%lf sigma=%lf t=%lf : nGenAll=%zu : -1SD=%lf +1SD=%lf\n",
timeSteps, nRuns, seed, S, mu, sigma, t, nGenAll, m1SD, p1SD);
// counting all generated values that fall within 1SD around the initial price
size_t cHit = 0;
// Monte Carlo simulation loop
for (size_t i = 0; i < nRuns; ++i)
{
// instantiate Geometric Brownian Motion (GBM) stochastic process
const auto& gbm = boost::shared_ptr<StochasticProcess>(new GeometricBrownianMotionProcess(S, mu, sigma));
// generate sequence of normally distributed random numbers from uniform distribution using Box-Muller transformation
if (i) seed = rand(); // for seed we use rand() for being able to reproduce the whole sequence (for DBG etc)...
MersenneTwisterUniformRng rng(seed);
//
typedef BoxMullerGaussianRng<MersenneTwisterUniformRng> MersenneBoxMuller;
MersenneBoxMuller boxMullerRng(rng);
//
RandomSequenceGenerator<MersenneBoxMuller> gsg(timeSteps, boxMullerRng);
// generate simulated path of stock price
PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> gbmPathGenerator(gbm, t, timeSteps, gsg, false);
const Path& samplePath = gbmPathGenerator.next().value;
// count the hits within 1SD around the initial price
for (size_t j = 1; j <= timeSteps; ++j) // BUGFIXED: now skipping 0th element, but still taking timeSteps elements
{
const auto Sx = samplePath.at(j);
if (Sx >= m1SD && Sx <= p1SD) ++cHit;
}
}
// result:
const double pct = double(cHit) / nGenAll * 100;
printf("cHit=%zu/%zu(%.3lf%%)\n", cHit, nGenAll, pct);
}
int main(int argc, char* argv[])
{
const size_t timeSteps = argc > 1 ? atol(argv[1]) : 1; // 1+
size_t nRuns = argc > 2 ? atol(argv[2]) : 0; // 0 means to take nRuns = 1000000 / timeSteps
if (!nRuns) nRuns = 1000000 / timeSteps;
const size_t seed = argc > 3 ? atol(argv[3]) : 0; // initial seed for the RNG. 0 means random initial seed
const double S = argc > 4 ? atof(argv[4]) : 100.0; // initial price
const double mu = argc > 5 ? atof(argv[5]) : 0.0; // drift
const double sigma = argc > 6 ? atof(argv[6]) : 0.30; // stddev, ie. volatility / 100
const double t = argc > 7 ? atof(argv[7]) : 1.0; // time in yrs
Testing_GBM_of_QuantLib(timeSteps, nRuns, seed, S, mu, sigma, t);
return 0;
}
Last edited: