The LogNormal isn't a solution to GBM. It's close--close enough that people say it a lot--but it's not an exact solution. Apply Fokker-Plank to the GBM SDE to get a PDE, and plug in the LogNormal, and you'll find it doesn't actually solve the PDE.
What do you conclude from that fact? For N=1 both lognormal and GBM _should_ (ie. has to, must) give the same result, where N is the number of intervals (aka "bars") in t. As demonstrated, GBM fails this basic test.
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.
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: 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!
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?
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!
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.
Here's C++ test code for testing QuantLib's GBM. 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; }
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%
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; }