Menu
▾
▴
quantlib-users
|
From: U.Mutlu <um...@mu...> - 2023-09-09 09:25:31
|
A short standalone test code in C++ with test results posted online at https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368 demonstrates that the GBM method in QuantLib is fatally buggy and has been so since start. Can the experts please check it and comment it. It's suspicious that such a flaw in a very important part of the library (GBM and Monte Carlo) got undetected by the experts and the user community for more that 20+ years. Can the experts verify / confirm / duplicate the findings made there? Or is maybe that test code itself buggy or the test method maybe unscientific? |
|
From: Ioannis R. <qua...@de...> - 2023-09-09 14:44:22
|
If you search within the QuantLib code for BoxMullerGaussianRng, you will see it is used only in the experimental folder. It is therefore not surprising if it doesn't produce the expected results. I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is used extensively in other areas of QuantLib. This type expands to InverseCumulativeRsg< RandomSequenceGenerator< MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good results. Ioannis Rigopoulos, founder of deriscope.com On 9/9/2023 11:25 AM, U.Mutlu wrote: > A short standalone test code in C++ with test results posted online at > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368 > > demonstrates that the GBM method in QuantLib > is fatally buggy and has been so since start. > > Can the experts please check it and comment it. > It's suspicious that such a flaw in a very important part of the > library (GBM and Monte Carlo) > got undetected by the experts and the user community for more that 20+ > years. > > Can the experts verify / confirm / duplicate the findings made there? > Or is maybe that test code itself buggy or the test method maybe > unscientific? > > > > _______________________________________________ > QuantLib-users mailing list > Qua...@li... > https://lists.sourceforge.net/lists/listinfo/quantlib-users -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: Peter C. <pca...@gm...> - 2023-09-09 18:25:05
|
Also notice that the (Ito-) solution of the SDE of
GeometricBrownianMotionProcess
dS(t, S)= \mu S dt + \sigma S dW_t.
see here
https://github.com/OpenSourceRisk/QuantLib/blob/00d2fced875622c1f52d9025b3ab01d04729eee8/ql/processes/geometricbrownianprocess.hpp#L36
is lognormally distributed with parameters
ln S(t) ~ N( \mu * t - 0.5 * \sigma^2 * t, \sigma^2 * t ).
Your testing code seems to assume
ln S(t) ~ N( \mu, \sigma^2 * t).
on the other hand. So maybe this is just a misunderstanding of the
parameters going into the constructor of
GeometricBrownianMotionProcess?
Best
Peter
On Sat, 9 Sept 2023 at 16:45, Ioannis Rigopoulos <qua...@de...> wrote:
>
> If you search within the QuantLib code for BoxMullerGaussianRng, you
> will see it is used only in the experimental folder. It is therefore not
> surprising if it doesn't produce the expected results.
>
> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which
> is used extensively in other areas of QuantLib.
>
> This type expands to InverseCumulativeRsg< RandomSequenceGenerator<
> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me
> good results.
>
> Ioannis Rigopoulos, founder of deriscope.com
>
> On 9/9/2023 11:25 AM, U.Mutlu wrote:
> > A short standalone test code in C++ with test results posted online at
> > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368
> >
> > demonstrates that the GBM method in QuantLib
> > is fatally buggy and has been so since start.
> >
> > Can the experts please check it and comment it.
> > It's suspicious that such a flaw in a very important part of the
> > library (GBM and Monte Carlo)
> > got undetected by the experts and the user community for more that 20+
> > years.
> >
> > Can the experts verify / confirm / duplicate the findings made there?
> > Or is maybe that test code itself buggy or the test method maybe
> > unscientific?
> >
> >
> >
> > _______________________________________________
> > QuantLib-users mailing list
> > Qua...@li...
> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
>
> --
> This email has been checked for viruses by Avast antivirus software.
> www.avast.com
>
>
> _______________________________________________
> QuantLib-users mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
|
|
From: U.Mutlu <um...@mu...> - 2023-09-09 20:32:41
|
I think Ito's lemma is the big culprit here. Ito's lemma is IMO misused in the whole field, and should especially not be used in GBM since GBM uses timeSteps... I'm 100% sure that one really does not need this Ito's lemma stuff, neither in GBM nor in the Black-Scholes-Merton (BSM) option pricing formula. Dunno about fields beyond these 2 cases. Just ask me if you need a proof for the BSM case. The GBM case (using the standard GBM algoritm with Ito's lemma) was proven to be incorrect just some weeks ago, again by me: 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 I wonder why the academia favors and is using this Ito's lemma stuff, as it's really unnecessary, IMO. Peter Caspers wrote on 09/09/23 20:24: > Also notice that the (Ito-) solution of the SDE of > GeometricBrownianMotionProcess > > dS(t, S)= \mu S dt + \sigma S dW_t. > > see here > > https://github.com/OpenSourceRisk/QuantLib/blob/00d2fced875622c1f52d9025b3ab01d04729eee8/ql/processes/geometricbrownianprocess.hpp#L36 > > is lognormally distributed with parameters > > ln S(t) ~ N( \mu * t - 0.5 * \sigma^2 * t, \sigma^2 * t ). > > Your testing code seems to assume > > ln S(t) ~ N( \mu, \sigma^2 * t). > > on the other hand. So maybe this is just a misunderstanding of the > parameters going into the constructor of > GeometricBrownianMotionProcess? > > Best > Peter > > On Sat, 9 Sept 2023 at 16:45, Ioannis Rigopoulos <qua...@de...> wrote: >> >> If you search within the QuantLib code for BoxMullerGaussianRng, you >> will see it is used only in the experimental folder. It is therefore not >> surprising if it doesn't produce the expected results. >> >> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which >> is used extensively in other areas of QuantLib. >> >> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me >> good results. >> >> Ioannis Rigopoulos, founder of deriscope.com >> >> On 9/9/2023 11:25 AM, U.Mutlu wrote: >>> A short standalone test code in C++ with test results posted online at >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368 >>> >>> demonstrates that the GBM method in QuantLib >>> is fatally buggy and has been so since start. >>> >>> Can the experts please check it and comment it. >>> It's suspicious that such a flaw in a very important part of the >>> library (GBM and Monte Carlo) >>> got undetected by the experts and the user community for more that 20+ >>> years. >>> >>> Can the experts verify / confirm / duplicate the findings made there? >>> Or is maybe that test code itself buggy or the test method maybe >>> unscientific? >>> >>> >>> >>> _______________________________________________ >>> QuantLib-users mailing list >>> Qua...@li... >>> https://lists.sourceforge.net/lists/listinfo/quantlib-users >> >> -- >> This email has been checked for viruses by Avast antivirus software. >> www.avast.com >> >> >> _______________________________________________ >> QuantLib-users mailing list >> Qua...@li... >> https://lists.sourceforge.net/lists/listinfo/quantlib-users > > > _______________________________________________ > QuantLib-users mailing list > Qua...@li... > https://lists.sourceforge.net/lists/listinfo/quantlib-users > |
|
From: U.Mutlu <um...@mu...> - 2023-09-10 09:30:07
|
As said in other posting here, after fixing the test program by skipping the first item (the initial price) in the generated sample path, BoxMullerGaussianRng now passes the said test. The bugfixed test code and the new test results can be found here: https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 That important fact that the generated sample path contains not timeSteps elements but 1 + timeSteps elements needs to be documented in the library doc. For example on this API doc page one normally would expect to find such an important information, but it's missing: https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html If you or someone else can change/extend the test program by using the suggested alternative(s) to BoxMullerGaussianRng, I would be happy to hear about it. Thx. Ioannis Rigopoulos wrote on 09/09/23 16:28: > If you search within the QuantLib code for BoxMullerGaussianRng, you will see > it is used only in the experimental folder. It is therefore not surprising if > it doesn't produce the expected results. > > I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is used > extensively in other areas of QuantLib. > > This type expands to InverseCumulativeRsg< RandomSequenceGenerator< > MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good > results. > > Ioannis Rigopoulos, founder of deriscope.com > > On 9/9/2023 11:25 AM, U.Mutlu wrote: >> A short standalone test code in C++ with test results posted online at >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368 >> >> demonstrates that the GBM method in QuantLib >> is fatally buggy and has been so since start. >> >> Can the experts please check it and comment it. >> It's suspicious that such a flaw in a very important part of the library >> (GBM and Monte Carlo) >> got undetected by the experts and the user community for more that 20+ years. >> >> Can the experts verify / confirm / duplicate the findings made there? >> Or is maybe that test code itself buggy or the test method maybe unscientific? |
|
From: Ioannis R. <qua...@de...> - 2023-09-10 11:46:30
|
Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - does not seem to be used in other areas of the standard part of QuantLib. Your numerical results below ... indicate that the simulated values for ln(S(t)) produced with timeSteps = 1 are very likely normally distributed with mean (I use your notation) ln(S(0)) + mu*t and standard deviation sigma * sqrt(t) This result _*is consistent*_ with: a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw _*and*_ b) the fact that in QuantLib the PathGenerator.next().value returns the result of the SDE expression _*without *_applying the ITO correction associated with the fact that dt is not infinitesimal. The b) is also responsible for the lack of convergence in your output towards the theoretical target of 68.27% You would get the correct convergence if you modified your code by using the expressions below: const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * -1.0); // Sx at -1SD const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * 1.0); // Sx at +1SD new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as Peter also pointed out) Again, due to b) one can produce the correct simulated values of a GBM diffused quantity (such as a stock price in the GBM model) by using N time steps, with N very large. Using N = 1 (like in your example), the simulated values will still be lognormally distributed (whence your good result with N = 1), but will be centered at a wrong mean and thus will _*fail *_to represent the correct values expected by the GBM SDE. Ioannis On 9/10/2023 11:29 AM, U.Mutlu wrote: > As said in other posting here, after fixing the test program > by skipping the first item (the initial price) in the > generated sample path, BoxMullerGaussianRng now passes the said test. > > The bugfixed test code and the new test results can be found here: > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 > > > That important fact that the generated sample path contains > not timeSteps elements but 1 + timeSteps elements needs > to be documented in the library doc. > For example on this API doc page one normally would expect > to find such an important information, but it's missing: > https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html > > > If you or someone else can change/extend the test program by using > the suggested alternative(s) to BoxMullerGaussianRng, I would be happy > to hear about it. Thx. > > > Ioannis Rigopoulos wrote on 09/09/23 16:28: >> If you search within the QuantLib code for BoxMullerGaussianRng, you >> will see >> it is used only in the experimental folder. It is therefore not >> surprising if >> it doesn't produce the expected results. >> >> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, >> which is used >> extensively in other areas of QuantLib. >> >> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me >> good >> results. >> >> Ioannis Rigopoulos, founder of deriscope.com >> >> On 9/9/2023 11:25 AM, U.Mutlu wrote: >>> A short standalone test code in C++ with test results posted online at >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861368 >>> >>> >>> demonstrates that the GBM method in QuantLib >>> is fatally buggy and has been so since start. >>> >>> Can the experts please check it and comment it. >>> It's suspicious that such a flaw in a very important part of the >>> library >>> (GBM and Monte Carlo) >>> got undetected by the experts and the user community for more that >>> 20+ years. >>> >>> Can the experts verify / confirm / duplicate the findings made there? >>> Or is maybe that test code itself buggy or the test method maybe >>> unscientific? > -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: U.Mutlu <um...@mu...> - 2023-09-10 13:01:01
|
Ioannis Rigopoulos wrote on 09/10/23 13:06: > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) Thank you for your analysis. I guess you mean "mu * t - 0.5 * sigma * sigma * t", though in our example with t=1 it doesn't make any difference. It gives for timeSteps=1 this result: $ ./Testing_GBM_of_QuantLib.exe 1 timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 +1SD=129.046162 cHit=663067/1000000(66.307%) This is far from the expected value of 68.2689% for timeStep=1. So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) Cf. also the following showing similar results for GBM, which then was compared to the standard lognormal calculation: https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied. Ioannis Rigopoulos wrote on 09/10/23 13:06: > Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - does > not seem to be used in other areas of the standard part of QuantLib. > > Your numerical results below ... > > indicate that the simulated values for ln(S(t)) produced with timeSteps = 1 > are very likely normally distributed with mean (I use your notation) ln(S(0)) > + mu*t and standard deviation sigma * sqrt(t) > > This result _*is consistent*_ with: > > a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw > > _*and*_ > > b) the fact that in QuantLib the PathGenerator.next().value returns the result > of the SDE expression _*without *_applying the ITO correction associated with > the fact that dt is not infinitesimal. > > The b) is also responsible for the lack of convergence in your output towards > the theoretical target of 68.27% > > You would get the correct convergence if you modified your code by using the > expressions below: > > const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * > -1.0); // Sx at -1SD > const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * > 1.0); // Sx at +1SD > > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as > Peter also pointed out) > > Again, due to b) one can produce the correct simulated values of a GBM > diffused quantity (such as a stock price in the GBM model) by using N time > steps, with N very large. Using N = 1 (like in your example), the simulated > values will still be lognormally distributed (whence your good result with N = > 1), but will be centered at a wrong mean and thus will _*fail *_to represent > the correct values expected by the GBM SDE. > > Ioannis > > On 9/10/2023 11:29 AM, U.Mutlu wrote: >> As said in other posting here, after fixing the test program >> by skipping the first item (the initial price) in the >> generated sample path, BoxMullerGaussianRng now passes the said test. >> >> The bugfixed test code and the new test results can be found here: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >> >> >> That important fact that the generated sample path contains >> not timeSteps elements but 1 + timeSteps elements needs >> to be documented in the library doc. >> For example on this API doc page one normally would expect >> to find such an important information, but it's missing: >> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >> >> If you or someone else can change/extend the test program by using >> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >> to hear about it. Thx. >> >> >> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>> If you search within the QuantLib code for BoxMullerGaussianRng, you will see >>> it is used only in the experimental folder. It is therefore not surprising if >>> it doesn't produce the expected results. >>> >>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is used >>> extensively in other areas of QuantLib. >>> >>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good >>> results. >>> >>> Ioannis Rigopoulos, founder of deriscope.com |
|
From: Ioannis R. <qua...@de...> - 2023-09-10 17:30:05
|
Yes, I mean mu * t - 0.5 * sigma * sigma * t Are you sure, you also adjusted the boundaries corresponding to plus and minus 1 SD? I ran an Excel spreadsheet by modifying a spreadsheet I had published several years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk that uses the GeometricBrownianMotionProcess , albeit with the PseudoRandom::rsg_type generator. At any case, I got convergence to 68.27% as I increased the number of time steps and keeping the number of simulations fixed to 100,000, provided I used the correct (Ito-adjusted) means, as below: As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, just like those you use. The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. The simulation results (percentage of values being within the min and max boundaries) are below: As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time steps, which is close to the expected 68.27%. The No Ito-adjusted percentage seems instead to converge to a different value around 67.73% On a different topic, I do not understand your statement "/I've never seen that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where m is the mean of a normally distributed variable. This means that one needs to calculate m in order to do the mentioned computation. If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then ln(S(t)) at any time t will be normally distributed with mean m = ln(S(0)) + mu*t - 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. Do you have any reason to believe that the Ito lemma is somehow not valid so that m = ln(S(0)) + mu*t? Ioannis On 9/10/2023 3:00 PM, U.Mutlu wrote: > Ioannis Rigopoulos wrote on 09/10/23 13:06: > > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) > > Thank you for your analysis. > I guess you mean "mu * t - 0.5 * sigma * sigma * t", > though in our example with t=1 it doesn't make any difference. > It gives for timeSteps=1 this result: > > $ ./Testing_GBM_of_QuantLib.exe 1 > timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 > sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 > +1SD=129.046162 > cHit=663067/1000000(66.307%) > > This is far from the expected value of 68.2689% for timeStep=1. > So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) > > Cf. also the following showing similar results for GBM, > which then was compared to the standard lognormal calculation: > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 > > > And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed > with Ito's lemma applied. > > > Ioannis Rigopoulos wrote on 09/10/23 13:06: >> Thank you for testing the BoxMullerGaussianRng code, which - as I >> wrote - does >> not seem to be used in other areas of the standard part of QuantLib. >> >> Your numerical results below ... >> >> indicate that the simulated values for ln(S(t)) produced with >> timeSteps = 1 >> are very likely normally distributed with mean (I use your notation) >> ln(S(0)) >> + mu*t and standard deviation sigma * sqrt(t) >> >> This result _*is consistent*_ with: >> >> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >> >> _*and*_ >> >> b) the fact that in QuantLib the PathGenerator.next().value returns >> the result >> of the SDE expression _*without *_applying the ITO correction >> associated with >> the fact that dt is not infinitesimal. >> >> The b) is also responsible for the lack of convergence in your output >> towards >> the theoretical target of 68.27% >> >> You would get the correct convergence if you modified your code by >> using the >> expressions below: >> >> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >> sqrt(t) * >> -1.0); // Sx at -1SD >> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >> sqrt(t) * >> 1.0); // Sx at +1SD >> >> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >> (as >> Peter also pointed out) >> >> Again, due to b) one can produce the correct simulated values of a GBM >> diffused quantity (such as a stock price in the GBM model) by using N >> time >> steps, with N very large. Using N = 1 (like in your example), the >> simulated >> values will still be lognormally distributed (whence your good result >> with N = >> 1), but will be centered at a wrong mean and thus will _*fail *_to >> represent >> the correct values expected by the GBM SDE. >> >> Ioannis >> >> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>> As said in other posting here, after fixing the test program >>> by skipping the first item (the initial price) in the >>> generated sample path, BoxMullerGaussianRng now passes the said test. >>> >>> The bugfixed test code and the new test results can be found here: >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>> >>> >>> >>> That important fact that the generated sample path contains >>> not timeSteps elements but 1 + timeSteps elements needs >>> to be documented in the library doc. >>> For example on this API doc page one normally would expect >>> to find such an important information, but it's missing: >>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>> >>> >>> If you or someone else can change/extend the test program by using >>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>> to hear about it. Thx. >>> >>> >>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>> If you search within the QuantLib code for BoxMullerGaussianRng, >>>> you will see >>>> it is used only in the experimental folder. It is therefore not >>>> surprising if >>>> it doesn't produce the expected results. >>>> >>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, >>>> which is used >>>> extensively in other areas of QuantLib. >>>> >>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives >>>> me good >>>> results. >>>> >>>> Ioannis Rigopoulos, founder of deriscope.com > -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: Ioannis R. <qua...@de...> - 2023-09-10 18:17:47
|
I think, I understand now your confusion regarding the usage of Ito lemma. The mean m = ln(S(0)) + mu*t - 0.5*sigma^2*t is the mean of ln(S(t)) at time t. Formally, m := ||ln(S(t))|| = ln(S(0)) + mu*t - 0.5*sigma^2*t , where ||x|| denotes the mean of the random variable x. But, due to Ito again, the above implies: m' := ||S(t)|| = S(0)*exp(mu*t) which means the mean m' of S(t) does not contain the "Ito term" 0.5*sigma^2*t. When investigating the 68.27% rule, we always look in a normally distributed variable, in this case the ln(S(t)), which is centered around m (not around m'). But when we speak in terms of S(t) - rather than of ln(S(t)) - the distribution of S(t) is centered around m', which indeed does not involve the Ito term. I guess, many practitioners prefer to refer to S(t) rather than to ln(S(t)) and therefore also refer to the 1SD double-side interval around m', giving the impression that the Ito term is not present. But if one needs to calculate the precise min and max values of the 1SD double-side interval in the S(t) space, one should use: min of S(t) = m'*exp(mu*t - 0.5*sigma^2*t - sigma*t¹ᐟ²) max of S(t) = m'*exp(mu*t - 0.5*sigma^2*t + sigma*t¹ᐟ²) The above expressions reintroduce the Ito term. Ioannis On 9/10/2023 7:29 PM, Ioannis Rigopoulos wrote: > > Yes, I mean mu * t - 0.5 * sigma * sigma * t > > Are you sure, you also adjusted the boundaries corresponding to plus > and minus 1 SD? > > I ran an Excel spreadsheet by modifying a spreadsheet I had published > several years ago at > https://blog.deriscope.com/index.php/en/excel-value-at-risk that uses > the GeometricBrownianMotionProcess , albeit with the > PseudoRandom::rsg_type generator. > > At any case, I got convergence to 68.27% as I increased the number of > time steps and keeping the number of simulations fixed to 100,000, > provided I used the correct (Ito-adjusted) means, as below: > > As you see above, the "No Ito-adjusted" boundaries are 74.08 and > 134.99, just like those you use. > > The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting > the no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. > > The simulation results (percentage of values being within the min and > max boundaries) are below: > > As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 > time steps, which is close to the expected 68.27%. > > The No Ito-adjusted percentage seems instead to converge to a > different value around 67.73% > > On a different topic, I do not understand your statement "/I've never > seen that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". > > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, > where m is the mean of a normally distributed variable. This means > that one needs to calculate m in order to do the mentioned computation. > > If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then > ln(S(t)) at any time t will be normally distributed with mean m = > ln(S(0)) + mu*t - 0.5*sigma^2*t , where the last term is a consequence > of the Ito lemma. > > Do you have any reason to believe that the Ito lemma is somehow not > valid so that m = ln(S(0)) + mu*t? > > Ioannis > > On 9/10/2023 3:00 PM, U.Mutlu wrote: > >> Ioannis Rigopoulos wrote on 09/10/23 13:06: >> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >> >> Thank you for your analysis. >> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >> though in our example with t=1 it doesn't make any difference. >> It gives for timeSteps=1 this result: >> >> $ ./Testing_GBM_of_QuantLib.exe 1 >> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : >> -1SD=70.822035 +1SD=129.046162 >> cHit=663067/1000000(66.307%) >> >> This is far from the expected value of 68.2689% for timeStep=1. >> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >> >> Cf. also the following showing similar results for GBM, >> which then was compared to the standard lognormal calculation: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >> >> >> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed >> with Ito's lemma applied. >> >> >> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>> Thank you for testing the BoxMullerGaussianRng code, which - as I >>> wrote - does >>> not seem to be used in other areas of the standard part of QuantLib. >>> >>> Your numerical results below ... >>> >>> indicate that the simulated values for ln(S(t)) produced with >>> timeSteps = 1 >>> are very likely normally distributed with mean (I use your notation) >>> ln(S(0)) >>> + mu*t and standard deviation sigma * sqrt(t) >>> >>> This result _*is consistent*_ with: >>> >>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>> >>> _*and*_ >>> >>> b) the fact that in QuantLib the PathGenerator.next().value returns >>> the result >>> of the SDE expression _*without *_applying the ITO correction >>> associated with >>> the fact that dt is not infinitesimal. >>> >>> The b) is also responsible for the lack of convergence in your >>> output towards >>> the theoretical target of 68.27% >>> >>> You would get the correct convergence if you modified your code by >>> using the >>> expressions below: >>> >>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>> sqrt(t) * >>> -1.0); // Sx at -1SD >>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>> sqrt(t) * >>> 1.0); // Sx at +1SD >>> >>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, >>> sigma) (as >>> Peter also pointed out) >>> >>> Again, due to b) one can produce the correct simulated values of a GBM >>> diffused quantity (such as a stock price in the GBM model) by using >>> N time >>> steps, with N very large. Using N = 1 (like in your example), the >>> simulated >>> values will still be lognormally distributed (whence your good >>> result with N = >>> 1), but will be centered at a wrong mean and thus will _*fail *_to >>> represent >>> the correct values expected by the GBM SDE. >>> >>> Ioannis >>> >>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>> As said in other posting here, after fixing the test program >>>> by skipping the first item (the initial price) in the >>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>> >>>> The bugfixed test code and the new test results can be found here: >>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>> >>>> >>>> >>>> That important fact that the generated sample path contains >>>> not timeSteps elements but 1 + timeSteps elements needs >>>> to be documented in the library doc. >>>> For example on this API doc page one normally would expect >>>> to find such an important information, but it's missing: >>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>> >>>> >>>> If you or someone else can change/extend the test program by using >>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>> to hear about it. Thx. >>>> >>>> >>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>> If you search within the QuantLib code for BoxMullerGaussianRng, >>>>> you will see >>>>> it is used only in the experimental folder. It is therefore not >>>>> surprising if >>>>> it doesn't produce the expected results. >>>>> >>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, >>>>> which is used >>>>> extensively in other areas of QuantLib. >>>>> >>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives >>>>> me good >>>>> results. >>>>> >>>>> Ioannis Rigopoulos, founder of deriscope.com >> > > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > Virus-free.www.avast.com > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > > > <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2> > > > _______________________________________________ > QuantLib-users mailing list > Qua...@li... > https://lists.sourceforge.net/lists/listinfo/quantlib-users -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 06:36:45
|
Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you very much. I need some time to study and verify the reults, but it looks very good to solve the confusion and mystery I had experienced when I started this research some weeks ago. I'll try to replicate your results and will report later in detail here. Ioannis Rigopoulos wrote on 09/10/23 19:29: > On a different topic, I do not understand your statement "/I've never seen > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where Sorry for the consfusion, it just means the stock price (Sx after time t, ie. S(t) or S_t) for the lower boundary at -1SD and for the upper boundary at +1SD, ie. your "min of S(t)" and "max of S(t)" in your other posting. In future I'll use a more common term like yours, sorry. Thx Ioannis Rigopoulos wrote on 09/10/23 19:29: > Yes, I mean mu * t - 0.5 * sigma * sigma * t > > Are you sure, you also adjusted the boundaries corresponding to plus and minus > 1 SD? > > I ran an Excel spreadsheet by modifying a spreadsheet I had published several > years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk that > uses the GeometricBrownianMotionProcess , albeit with the > PseudoRandom::rsg_type generator. > > At any case, I got convergence to 68.27% as I increased the number of time > steps and keeping the number of simulations fixed to 100,000, provided I used > the correct (Ito-adjusted) means, as below: > > As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, just > like those you use. > > The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the > no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. > > The simulation results (percentage of values being within the min and max > boundaries) are below: > > As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time steps, > which is close to the expected 68.27%. > > The No Ito-adjusted percentage seems instead to converge to a different value > around 67.73% > > On a different topic, I do not understand your statement "/I've never seen > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". > > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where m is > the mean of a normally distributed variable. This means that one needs to > calculate m in order to do the mentioned computation. > > If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then ln(S(t)) > at any time t will be normally distributed with mean m = ln(S(0)) + mu*t - > 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. > > Do you have any reason to believe that the Ito lemma is somehow not valid so > that m = ln(S(0)) + mu*t? > > Ioannis > > On 9/10/2023 3:00 PM, U.Mutlu wrote: > >> Ioannis Rigopoulos wrote on 09/10/23 13:06: >> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >> >> Thank you for your analysis. >> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >> though in our example with t=1 it doesn't make any difference. >> It gives for timeSteps=1 this result: >> >> $ ./Testing_GBM_of_QuantLib.exe 1 >> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 >> +1SD=129.046162 >> cHit=663067/1000000(66.307%) >> >> This is far from the expected value of 68.2689% for timeStep=1. >> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >> >> Cf. also the following showing similar results for GBM, >> which then was compared to the standard lognormal calculation: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >> >> >> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed with >> Ito's lemma applied. >> >> >> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>> Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - does >>> not seem to be used in other areas of the standard part of QuantLib. >>> >>> Your numerical results below ... >>> >>> indicate that the simulated values for ln(S(t)) produced with timeSteps = 1 >>> are very likely normally distributed with mean (I use your notation) ln(S(0)) >>> + mu*t and standard deviation sigma * sqrt(t) >>> >>> This result _*is consistent*_ with: >>> >>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>> >>> _*and*_ >>> >>> b) the fact that in QuantLib the PathGenerator.next().value returns the result >>> of the SDE expression _*without *_applying the ITO correction associated with >>> the fact that dt is not infinitesimal. >>> >>> The b) is also responsible for the lack of convergence in your output towards >>> the theoretical target of 68.27% >>> >>> You would get the correct convergence if you modified your code by using the >>> expressions below: >>> >>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>> -1.0); // Sx at -1SD >>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>> 1.0); // Sx at +1SD >>> >>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as >>> Peter also pointed out) >>> >>> Again, due to b) one can produce the correct simulated values of a GBM >>> diffused quantity (such as a stock price in the GBM model) by using N time >>> steps, with N very large. Using N = 1 (like in your example), the simulated >>> values will still be lognormally distributed (whence your good result with N = >>> 1), but will be centered at a wrong mean and thus will _*fail *_to represent >>> the correct values expected by the GBM SDE. >>> >>> Ioannis >>> >>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>> As said in other posting here, after fixing the test program >>>> by skipping the first item (the initial price) in the >>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>> >>>> The bugfixed test code and the new test results can be found here: >>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>> >>>> >>>> >>>> That important fact that the generated sample path contains >>>> not timeSteps elements but 1 + timeSteps elements needs >>>> to be documented in the library doc. >>>> For example on this API doc page one normally would expect >>>> to find such an important information, but it's missing: >>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>> >>>> If you or someone else can change/extend the test program by using >>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>> to hear about it. Thx. >>>> >>>> >>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>> If you search within the QuantLib code for BoxMullerGaussianRng, you will >>>>> see >>>>> it is used only in the experimental folder. It is therefore not >>>>> surprising if >>>>> it doesn't produce the expected results. >>>>> >>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is >>>>> used >>>>> extensively in other areas of QuantLib. >>>>> >>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good >>>>> results. >>>>> >>>>> Ioannis Rigopoulos, founder of deriscope.com >> > > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > Virus-free.www.avast.com > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > > > > > > _______________________________________________ > QuantLib-users mailing list > Qua...@li... > https://lists.sourceforge.net/lists/listinfo/quantlib-users > |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 08:20:08
|
Dear Ioannis Rigopoulos, I now analysed your table and must say that you seem to have a much different expectation to, and understanding of that table. The values in your table for timeSteps > 1 can definitely not be correct, b/c why do you think it has to be 68.27 for all timeSteps? This is illogical b/c the more timeSteps the more the hitrate% must be, like in my table. It saturates at about 84.93% Sorry, I have to retract my initial enthusiastic judment of your work below, as it's fundamentally flawed & wrong, and I should have been more careful when I quickly & briefly inspected your posting in the early morning here. Since your results in the table for timeSteps > 1 are wrong, then your formulas must be wrong too, so I had to stop my analysis right there. Please just answer why you think it has to be 68.27% for all timeSteps? U.Mutlu wrote on 09/11/23 08:36: > Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you very much. > I need some time to study and verify the reults, but it looks very good > to solve the confusion and mystery I had experienced when I started this > research some weeks ago. > I'll try to replicate your results and will report later in detail here. > > Ioannis Rigopoulos wrote on 09/10/23 19:29: > > On a different topic, I do not understand your statement "/I've never seen > > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". > > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where > > Sorry for the consfusion, it just means the stock price (Sx after time t, > ie. S(t) or S_t) for the lower boundary at -1SD and for the upper boundary at > +1SD, > ie. your "min of S(t)" and "max of S(t)" in your other posting. > In future I'll use a more common term like yours, sorry. > > Thx > > > Ioannis Rigopoulos wrote on 09/10/23 19:29: >> Yes, I mean mu * t - 0.5 * sigma * sigma * t >> >> Are you sure, you also adjusted the boundaries corresponding to plus and minus >> 1 SD? >> >> I ran an Excel spreadsheet by modifying a spreadsheet I had published several >> years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk that >> uses the GeometricBrownianMotionProcess , albeit with the >> PseudoRandom::rsg_type generator. >> >> At any case, I got convergence to 68.27% as I increased the number of time >> steps and keeping the number of simulations fixed to 100,000, provided I used >> the correct (Ito-adjusted) means, as below: >> >> As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, just >> like those you use. >> >> The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the >> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. >> >> The simulation results (percentage of values being within the min and max >> boundaries) are below: >> >> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time steps, >> which is close to the expected 68.27%. >> >> The No Ito-adjusted percentage seems instead to converge to a different value >> around 67.73% >> >> On a different topic, I do not understand your statement "/I've never seen >> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >> >> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where m is >> the mean of a normally distributed variable. This means that one needs to >> calculate m in order to do the mentioned computation. >> >> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then ln(S(t)) >> at any time t will be normally distributed with mean m = ln(S(0)) + mu*t - >> 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. >> >> Do you have any reason to believe that the Ito lemma is somehow not valid so >> that m = ln(S(0)) + mu*t? >> >> Ioannis >> >> On 9/10/2023 3:00 PM, U.Mutlu wrote: >> >>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>> >>> Thank you for your analysis. >>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>> though in our example with t=1 it doesn't make any difference. >>> It gives for timeSteps=1 this result: >>> >>> $ ./Testing_GBM_of_QuantLib.exe 1 >>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 >>> +1SD=129.046162 >>> cHit=663067/1000000(66.307%) >>> >>> This is far from the expected value of 68.2689% for timeStep=1. >>> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >>> >>> Cf. also the following showing similar results for GBM, >>> which then was compared to the standard lognormal calculation: >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>> >>> >>> >>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed with >>> Ito's lemma applied. >>> >>> >>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>> Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - >>>> does >>>> not seem to be used in other areas of the standard part of QuantLib. >>>> >>>> Your numerical results below ... >>>> >>>> indicate that the simulated values for ln(S(t)) produced with timeSteps = 1 >>>> are very likely normally distributed with mean (I use your notation) ln(S(0)) >>>> + mu*t and standard deviation sigma * sqrt(t) >>>> >>>> This result _*is consistent*_ with: >>>> >>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>> >>>> _*and*_ >>>> >>>> b) the fact that in QuantLib the PathGenerator.next().value returns the >>>> result >>>> of the SDE expression _*without *_applying the ITO correction associated with >>>> the fact that dt is not infinitesimal. >>>> >>>> The b) is also responsible for the lack of convergence in your output towards >>>> the theoretical target of 68.27% >>>> >>>> You would get the correct convergence if you modified your code by using the >>>> expressions below: >>>> >>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>>> -1.0); // Sx at -1SD >>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>>> 1.0); // Sx at +1SD >>>> >>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as >>>> Peter also pointed out) >>>> >>>> Again, due to b) one can produce the correct simulated values of a GBM >>>> diffused quantity (such as a stock price in the GBM model) by using N time >>>> steps, with N very large. Using N = 1 (like in your example), the simulated >>>> values will still be lognormally distributed (whence your good result with >>>> N = >>>> 1), but will be centered at a wrong mean and thus will _*fail *_to represent >>>> the correct values expected by the GBM SDE. >>>> >>>> Ioannis >>>> >>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>> As said in other posting here, after fixing the test program >>>>> by skipping the first item (the initial price) in the >>>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>>> >>>>> The bugfixed test code and the new test results can be found here: >>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>> >>>>> >>>>> >>>>> >>>>> That important fact that the generated sample path contains >>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>> to be documented in the library doc. >>>>> For example on this API doc page one normally would expect >>>>> to find such an important information, but it's missing: >>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>> >>>>> If you or someone else can change/extend the test program by using >>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>>> to hear about it. Thx. >>>>> >>>>> >>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>> If you search within the QuantLib code for BoxMullerGaussianRng, you will >>>>>> see >>>>>> it is used only in the experimental folder. It is therefore not >>>>>> surprising if >>>>>> it doesn't produce the expected results. >>>>>> >>>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is >>>>>> used >>>>>> extensively in other areas of QuantLib. >>>>>> >>>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good >>>>>> results. >>>>>> >>>>>> Ioannis Rigopoulos, founder of deriscope.com >>> >> >> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> >> >> Virus-free.www.avast.com >> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> >> >> >> >> >> >> >> _______________________________________________ >> QuantLib-users mailing list >> Qua...@li... >> https://lists.sourceforge.net/lists/listinfo/quantlib-users >> > > |
|
From: Ioannis R. <qua...@de...> - 2023-09-11 09:48:14
|
Our SDE for the evolution of S is:
dS = μSdt + σSdw
In your code you use:
PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>>
gbmPathGenerator(gbm, t, timeSteps, gsg, false);
const Path& samplePath = gbmPathGenerator.next().value;
The PathGenerator::next().value is designed to return an array of
numbers {S(T₁) , S(T₂) , ... , S(Tᵤ)}
where
u := timeSteps
Δt := t/u
Tᵢ := iΔt for i = 1,...,u
It follows:
Tᵤ = t = the final time when we are interested to know the simulated
values of S (in your test you use t = 1).
QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} incrementally
through a total number of u steps by adding the quantity ΔS := μSΔt +
σSΔw on the previous member of the sequence, with Δw := random number
drawn from N~(0,Δt¹ᐟ²).
The previous member of S(T₁) is assumed the T₀ = 0 and the initial value
S(T₀) = S(0) is assumed deterministic and known.
So, starting with the given S(T₀), QuantLib calculates the S(T₁) using
the formula:
S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw
Then given the S(T₁) as calculated above, it proceeds to calculate S(T₂) as:
S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw
For simplicity I have used the same notation for the Δw in the equations
above, but they are different.
After u steps, the final value S(t) = S(Tᵤ) is calculated as:
S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw
All equations above hold only in an approximate sense when Δt is finite.
They are accurate only in a limiting sense as Δt → 0.
It follows that the value S(t) will be calculated correctly only if Δt →
0, which is equivalent to u → ∞
If you use u = 1 (i.e. one time step), you will still match the 68.27%
rule when the 1SD interval is defined in accordance with a normal
distribution because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to
a normal distribution and if you build a histogram of all simulated
values you will notice that S(T₁) is distributed normally, which is
clearly not consistent with the SDE that predicts lognormal distribution.
The correct procedure is of course to use a very large u and define the
1SD interval in accordance with a lognormal distribution, as I described
in my previous email.
Then you will get convergence to the 68.27% rule.
By the way, this is how my histogram of S(t) (for t = 1) simulated
values looks like when I used 100,000 time steps and 100,000
simulations. It is clearly lognormal in shape. I am sure, if you do the
same when time steps = 1, you will get a normal distribution for S(t),
which would prove that time steps = 1 does not reproduce the
theoretically expected values for S(t).
Ioannis
On 9/11/2023 10:19 AM, U.Mutlu wrote:
> Dear Ioannis Rigopoulos, I now analysed your table and
> must say that you seem to have a much different expectation to,
> and understanding of that table.
>
> The values in your table for timeSteps > 1 can definitely
> not be correct, b/c why do you think it has to be 68.27 for all
> timeSteps?
> This is illogical b/c the more timeSteps the more the hitrate%
> must be, like in my table. It saturates at about 84.93%
>
> Sorry, I have to retract my initial enthusiastic judment
> of your work below, as it's fundamentally flawed & wrong,
> and I should have been more careful when I quickly & briefly
> inspected your posting in the early morning here.
>
> Since your results in the table for timeSteps > 1 are wrong,
> then your formulas must be wrong too, so I had to stop my analysis
> right there.
>
> Please just answer why you think it has to be 68.27% for all timeSteps?
>
>
> U.Mutlu wrote on 09/11/23 08:36:
>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you
>> very much.
>> I need some time to study and verify the reults, but it looks very good
>> to solve the confusion and mystery I had experienced when I started this
>> research some weeks ago.
>> I'll try to replicate your results and will report later in detail here.
>>
>> Ioannis Rigopoulos wrote on 09/10/23 19:29:
>> > On a different topic, I do not understand your statement "/I've
>> never seen
>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/".
>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD,
>> where
>>
>> Sorry for the consfusion, it just means the stock price (Sx after
>> time t,
>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper
>> boundary at
>> +1SD,
>> ie. your "min of S(t)" and "max of S(t)" in your other posting.
>> In future I'll use a more common term like yours, sorry.
>>
>> Thx
>>
>>
>> Ioannis Rigopoulos wrote on 09/10/23 19:29:
>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t
>>>
>>> Are you sure, you also adjusted the boundaries corresponding to plus
>>> and minus
>>> 1 SD?
>>>
>>> I ran an Excel spreadsheet by modifying a spreadsheet I had
>>> published several
>>> years ago at
>>> https://blog.deriscope.com/index.php/en/excel-value-at-risk that
>>> uses the GeometricBrownianMotionProcess , albeit with the
>>> PseudoRandom::rsg_type generator.
>>>
>>> At any case, I got convergence to 68.27% as I increased the number
>>> of time
>>> steps and keeping the number of simulations fixed to 100,000,
>>> provided I used
>>> the correct (Ito-adjusted) means, as below:
>>>
>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and
>>> 134.99, just
>>> like those you use.
>>>
>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by
>>> shifting the
>>> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t.
>>>
>>> The simulation results (percentage of values being within the min
>>> and max
>>> boundaries) are below:
>>>
>>> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000
>>> time steps,
>>> which is close to the expected 68.27%.
>>>
>>> The No Ito-adjusted percentage seems instead to converge to a
>>> different value
>>> around 67.73%
>>>
>>> On a different topic, I do not understand your statement "/I've
>>> never seen
>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/".
>>>
>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD,
>>> where m is
>>> the mean of a normally distributed variable. This means that one
>>> needs to
>>> calculate m in order to do the mentioned computation.
>>>
>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw,
>>> then ln(S(t))
>>> at any time t will be normally distributed with mean m = ln(S(0)) +
>>> mu*t -
>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma.
>>>
>>> Do you have any reason to believe that the Ito lemma is somehow not
>>> valid so
>>> that m = ln(S(0)) + mu*t?
>>>
>>> Ioannis
>>>
>>> On 9/10/2023 3:00 PM, U.Mutlu wrote:
>>>
>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06:
>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma)
>>>>
>>>> Thank you for your analysis.
>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t",
>>>> though in our example with t=1 it doesn't make any difference.
>>>> It gives for timeSteps=1 this result:
>>>>
>>>> $ ./Testing_GBM_of_QuantLib.exe 1
>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000
>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035
>>>> +1SD=129.046162
>>>> cHit=663067/1000000(66.307%)
>>>>
>>>> This is far from the expected value of 68.2689% for timeStep=1.
>>>> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-)
>>>>
>>>> Cf. also the following showing similar results for GBM,
>>>> which then was compared to the standard lognormal calculation:
>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548
>>>>
>>>>
>>>>
>>>>
>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get
>>>> computed with
>>>> Ito's lemma applied.
>>>>
>>>>
>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06:
>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as I
>>>>> wrote -
>>>>> does
>>>>> not seem to be used in other areas of the standard part of QuantLib.
>>>>>
>>>>> Your numerical results below ...
>>>>>
>>>>> indicate that the simulated values for ln(S(t)) produced with
>>>>> timeSteps = 1
>>>>> are very likely normally distributed with mean (I use your
>>>>> notation) ln(S(0))
>>>>> + mu*t and standard deviation sigma * sqrt(t)
>>>>>
>>>>> This result _*is consistent*_ with:
>>>>>
>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw
>>>>>
>>>>> _*and*_
>>>>>
>>>>> b) the fact that in QuantLib the PathGenerator.next().value
>>>>> returns the
>>>>> result
>>>>> of the SDE expression _*without *_applying the ITO correction
>>>>> associated with
>>>>> the fact that dt is not infinitesimal.
>>>>>
>>>>> The b) is also responsible for the lack of convergence in your
>>>>> output towards
>>>>> the theoretical target of 68.27%
>>>>>
>>>>> You would get the correct convergence if you modified your code by
>>>>> using the
>>>>> expressions below:
>>>>>
>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma *
>>>>> sqrt(t) *
>>>>> -1.0); // Sx at -1SD
>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma *
>>>>> sqrt(t) *
>>>>> 1.0); // Sx at +1SD
>>>>>
>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t,
>>>>> sigma) (as
>>>>> Peter also pointed out)
>>>>>
>>>>> Again, due to b) one can produce the correct simulated values of a
>>>>> GBM
>>>>> diffused quantity (such as a stock price in the GBM model) by
>>>>> using N time
>>>>> steps, with N very large. Using N = 1 (like in your example), the
>>>>> simulated
>>>>> values will still be lognormally distributed (whence your good
>>>>> result with
>>>>> N =
>>>>> 1), but will be centered at a wrong mean and thus will _*fail *_to
>>>>> represent
>>>>> the correct values expected by the GBM SDE.
>>>>>
>>>>> Ioannis
>>>>>
>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote:
>>>>>> As said in other posting here, after fixing the test program
>>>>>> by skipping the first item (the initial price) in the
>>>>>> generated sample path, BoxMullerGaussianRng now passes the said
>>>>>> test.
>>>>>>
>>>>>> The bugfixed test code and the new test results can be found here:
>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> That important fact that the generated sample path contains
>>>>>> not timeSteps elements but 1 + timeSteps elements needs
>>>>>> to be documented in the library doc.
>>>>>> For example on this API doc page one normally would expect
>>>>>> to find such an important information, but it's missing:
>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html
>>>>>>
>>>>>>
>>>>>> If you or someone else can change/extend the test program by using
>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be
>>>>>> happy
>>>>>> to hear about it. Thx.
>>>>>>
>>>>>>
>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28:
>>>>>>> If you search within the QuantLib code for BoxMullerGaussianRng,
>>>>>>> you will
>>>>>>> see
>>>>>>> it is used only in the experimental folder. It is therefore not
>>>>>>> surprising if
>>>>>>> it doesn't produce the expected results.
>>>>>>>
>>>>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type,
>>>>>>> which is
>>>>>>> used
>>>>>>> extensively in other areas of QuantLib.
>>>>>>>
>>>>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator<
>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and
>>>>>>> gives me good
>>>>>>> results.
>>>>>>>
>>>>>>> Ioannis Rigopoulos, founder of deriscope.com
>>>>
>>>
>>> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
>>>
>>>
>>> Virus-free.www.avast.com
>>> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> QuantLib-users mailing list
>>> Qua...@li...
>>> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>>>
>>
>>
>
--
This email has been checked for viruses by Avast antivirus software.
www.avast.com |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 11:34:00
|
Ioannis Rigopoulos wrote on 09/11/23 11:47: > By the way, this is how my histogram of S(t) (for t = 1) simulated values > looks like when I used 100,000 time steps and 100,000 simulations. It is > clearly lognormal in shape. Yes, it's lognormally. My values are lognormlly as well, You just restate the obvious fact. > I am sure, if you do the same when time steps = 1, > you will get a normal distribution for S(t), which would prove > that time steps = 1 does not reproduce the theoretically expected > values for S(t). Not true, simply from the fact that the set boundary is lognormally: -1SD: 74.08182206817179 +1SD: 134.9858807576003 Ie. it's not symmetric (Gaussian, normally distributed) around the initial mean (100), but is lognormally (and therefore unsymmetric). FYI: in my other test program (that does not use QuantLib) I get in such simulations exact results almost identical to the expected theoretical results! So, your math and arguing is wrong, IMO. You are too much fixated on timeSteps=1. Maybe you have a much different understanding of what timeSteps really means in practice. It's _not_ something like the number of simulations that you simply can change and expect to get the same result; no, it has a clearly defined context and gives different results for different timeSteps. And you haven't answered the question why you think it has to be 68.27% for all timeSteps (which IMO is wrong). For those interested, here are more details about the above said: https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 And, I today did a brief summary of these ongoing discussions on GBM here: https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919 Ioannis Rigopoulos wrote on 09/11/23 11:47: > Our SDE for the evolution of S is: > > dS = μSdt + σSdw > > In your code you use: > > > PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> > gbmPathGenerator(gbm, t, timeSteps, gsg, false); > const Path& samplePath = gbmPathGenerator.next().value; > > The PathGenerator::next().value is designed to return an array of numbers > {S(T₁) , S(T₂) , ... , S(Tᵤ)} > > where > > u := timeSteps > > Δt := t/u > > Tᵢ := iΔt for i = 1,...,u > > It follows: > > Tᵤ = t = the final time when we are interested to know the simulated values of > S (in your test you use t = 1). > > QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} incrementally > through a total number of u steps by adding the quantity ΔS := μSΔt + σSΔw on > the previous member of the sequence, with Δw := random number drawn from > N~(0,Δt¹ᐟ²). > > The previous member of S(T₁) is assumed the T₀ = 0 and the initial value S(T₀) > = S(0) is assumed deterministic and known. > > So, starting with the given S(T₀), QuantLib calculates the S(T₁) using the > formula: > > S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw > > Then given the S(T₁) as calculated above, it proceeds to calculate S(T₂) as: > > S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw > > For simplicity I have used the same notation for the Δw in the equations > above, but they are different. > > After u steps, the final value S(t) = S(Tᵤ) is calculated as: > > S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw > > All equations above hold only in an approximate sense when Δt is finite. > > They are accurate only in a limiting sense as Δt → 0. > > It follows that the value S(t) will be calculated correctly only if Δt → 0, > which is equivalent to u → ∞ > > If you use u = 1 (i.e. one time step), you will still match the 68.27% rule > when the 1SD interval is defined in accordance with a normal distribution > because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal > distribution and if you build a histogram of all simulated values you will > notice that S(T₁) is distributed normally, which is clearly not consistent > with the SDE that predicts lognormal distribution. > > The correct procedure is of course to use a very large u and define the 1SD > interval in accordance with a lognormal distribution, as I described in my > previous email. > > Then you will get convergence to the 68.27% rule. > > By the way, this is how my histogram of S(t) (for t = 1) simulated values > looks like when I used 100,000 time steps and 100,000 simulations. It is > clearly lognormal in shape. I am sure, if you do the same when time steps = 1, > you will get a normal distribution for S(t), which would prove that time steps > = 1 does not reproduce the theoretically expected values for S(t). > > Ioannis > > On 9/11/2023 10:19 AM, U.Mutlu wrote: >> Dear Ioannis Rigopoulos, I now analysed your table and >> must say that you seem to have a much different expectation to, >> and understanding of that table. >> >> The values in your table for timeSteps > 1 can definitely >> not be correct, b/c why do you think it has to be 68.27 for all timeSteps? >> This is illogical b/c the more timeSteps the more the hitrate% >> must be, like in my table. It saturates at about 84.93% >> >> Sorry, I have to retract my initial enthusiastic judment >> of your work below, as it's fundamentally flawed & wrong, >> and I should have been more careful when I quickly & briefly >> inspected your posting in the early morning here. >> >> Since your results in the table for timeSteps > 1 are wrong, >> then your formulas must be wrong too, so I had to stop my analysis right there. >> >> Please just answer why you think it has to be 68.27% for all timeSteps? >> >> >> U.Mutlu wrote on 09/11/23 08:36: >>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you very much. >>> I need some time to study and verify the reults, but it looks very good >>> to solve the confusion and mystery I had experienced when I started this >>> research some weeks ago. >>> I'll try to replicate your results and will report later in detail here. >>> >>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>> > On a different topic, I do not understand your statement "/I've never seen >>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where >>> >>> Sorry for the consfusion, it just means the stock price (Sx after time t, >>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper boundary at >>> +1SD, >>> ie. your "min of S(t)" and "max of S(t)" in your other posting. >>> In future I'll use a more common term like yours, sorry. >>> >>> Thx >>> >>> >>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t >>>> >>>> Are you sure, you also adjusted the boundaries corresponding to plus and >>>> minus >>>> 1 SD? >>>> >>>> I ran an Excel spreadsheet by modifying a spreadsheet I had published several >>>> years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk that >>>> uses the GeometricBrownianMotionProcess , albeit with the >>>> PseudoRandom::rsg_type generator. >>>> >>>> At any case, I got convergence to 68.27% as I increased the number of time >>>> steps and keeping the number of simulations fixed to 100,000, provided I used >>>> the correct (Ito-adjusted) means, as below: >>>> >>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, just >>>> like those you use. >>>> >>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the >>>> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. >>>> >>>> The simulation results (percentage of values being within the min and max >>>> boundaries) are below: >>>> >>>> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time >>>> steps, >>>> which is close to the expected 68.27%. >>>> >>>> The No Ito-adjusted percentage seems instead to converge to a different value >>>> around 67.73% >>>> >>>> On a different topic, I do not understand your statement "/I've never seen >>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>> >>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where m is >>>> the mean of a normally distributed variable. This means that one needs to >>>> calculate m in order to do the mentioned computation. >>>> >>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then >>>> ln(S(t)) >>>> at any time t will be normally distributed with mean m = ln(S(0)) + mu*t - >>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. >>>> >>>> Do you have any reason to believe that the Ito lemma is somehow not valid so >>>> that m = ln(S(0)) + mu*t? >>>> >>>> Ioannis >>>> >>>> On 9/10/2023 3:00 PM, U.Mutlu wrote: >>>> >>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>>>> >>>>> Thank you for your analysis. >>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>>>> though in our example with t=1 it doesn't make any difference. >>>>> It gives for timeSteps=1 this result: >>>>> >>>>> $ ./Testing_GBM_of_QuantLib.exe 1 >>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 >>>>> +1SD=129.046162 >>>>> cHit=663067/1000000(66.307%) >>>>> >>>>> This is far from the expected value of 68.2689% for timeStep=1. >>>>> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >>>>> >>>>> Cf. also the following showing similar results for GBM, >>>>> which then was compared to the standard lognormal calculation: >>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>>> >>>>> >>>>> >>>>> >>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed with >>>>> Ito's lemma applied. >>>>> >>>>> >>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - >>>>>> does >>>>>> not seem to be used in other areas of the standard part of QuantLib. >>>>>> >>>>>> Your numerical results below ... >>>>>> >>>>>> indicate that the simulated values for ln(S(t)) produced with timeSteps = 1 >>>>>> are very likely normally distributed with mean (I use your notation) >>>>>> ln(S(0)) >>>>>> + mu*t and standard deviation sigma * sqrt(t) >>>>>> >>>>>> This result _*is consistent*_ with: >>>>>> >>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>>>> >>>>>> _*and*_ >>>>>> >>>>>> b) the fact that in QuantLib the PathGenerator.next().value returns the >>>>>> result >>>>>> of the SDE expression _*without *_applying the ITO correction associated >>>>>> with >>>>>> the fact that dt is not infinitesimal. >>>>>> >>>>>> The b) is also responsible for the lack of convergence in your output >>>>>> towards >>>>>> the theoretical target of 68.27% >>>>>> >>>>>> You would get the correct convergence if you modified your code by using >>>>>> the >>>>>> expressions below: >>>>>> >>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>>>>> -1.0); // Sx at -1SD >>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * sqrt(t) * >>>>>> 1.0); // Sx at +1SD >>>>>> >>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as >>>>>> Peter also pointed out) >>>>>> >>>>>> Again, due to b) one can produce the correct simulated values of a GBM >>>>>> diffused quantity (such as a stock price in the GBM model) by using N time >>>>>> steps, with N very large. Using N = 1 (like in your example), the simulated >>>>>> values will still be lognormally distributed (whence your good result with >>>>>> N = >>>>>> 1), but will be centered at a wrong mean and thus will _*fail *_to >>>>>> represent >>>>>> the correct values expected by the GBM SDE. >>>>>> >>>>>> Ioannis >>>>>> >>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>>>> As said in other posting here, after fixing the test program >>>>>>> by skipping the first item (the initial price) in the >>>>>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>>>>> >>>>>>> The bugfixed test code and the new test results can be found here: >>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> That important fact that the generated sample path contains >>>>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>>>> to be documented in the library doc. >>>>>>> For example on this API doc page one normally would expect >>>>>>> to find such an important information, but it's missing: >>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>>>> >>>>>>> If you or someone else can change/extend the test program by using >>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>>>>> to hear about it. Thx. >>>>>>> >>>>>>> >>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>>>> If you search within the QuantLib code for BoxMullerGaussianRng, you will >>>>>>>> see >>>>>>>> it is used only in the experimental folder. It is therefore not >>>>>>>> surprising if >>>>>>>> it doesn't produce the expected results. >>>>>>>> >>>>>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, which is >>>>>>>> used >>>>>>>> extensively in other areas of QuantLib. >>>>>>>> >>>>>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me good >>>>>>>> results. >>>>>>>> >>>>>>>> Ioannis Rigopoulos, founder of deriscope.com |
|
From: Ioannis R. <qua...@de...> - 2023-09-11 11:43:41
|
Could you please send us a link to a file that contains the simulated values produced by your program when you use one time step? Since you use 1,000,000 simulations, the file should contain 1,000,000 numbers. On 9/11/2023 1:33 PM, U.Mutlu wrote: > Ioannis Rigopoulos wrote on 09/11/23 11:47: > > By the way, this is how my histogram of S(t) (for t = 1) simulated > values > > looks like when I used 100,000 time steps and 100,000 simulations. > It is > > clearly lognormal in shape. > > Yes, it's lognormally. My values are lognormlly as well, > You just restate the obvious fact. > > > I am sure, if you do the same when time steps = 1, > > you will get a normal distribution for S(t), which would prove > > that time steps = 1 does not reproduce the theoretically expected > > values for S(t). > > Not true, simply from the fact that the set boundary is lognormally: > -1SD: 74.08182206817179 > +1SD: 134.9858807576003 > Ie. it's not symmetric (Gaussian, normally distributed) around the > initial mean (100), but is lognormally (and therefore unsymmetric). > > FYI: in my other test program (that does not use QuantLib) > I get in such simulations exact results almost identical > to the expected theoretical results! So, your math and arguing is > wrong, IMO. > > You are too much fixated on timeSteps=1. Maybe you have > a much different understanding of what timeSteps really > means in practice. It's _not_ something like the number > of simulations that you simply can change and expect > to get the same result; no, it has a clearly defined context > and gives different results for different timeSteps. > > And you haven't answered the question why you think it has > to be 68.27% for all timeSteps (which IMO is wrong). > > For those interested, here are more details about the above said: > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 > > > And, I today did a brief summary of these ongoing discussions on GBM > here: > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919 > > > > Ioannis Rigopoulos wrote on 09/11/23 11:47: >> Our SDE for the evolution of S is: >> >> dS = μSdt + σSdw >> >> In your code you use: >> >> >> PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> >> gbmPathGenerator(gbm, t, timeSteps, gsg, false); >> const Path& samplePath = gbmPathGenerator.next().value; >> >> The PathGenerator::next().value is designed to return an array of >> numbers >> {S(T₁) , S(T₂) , ... , S(Tᵤ)} >> >> where >> >> u := timeSteps >> >> Δt := t/u >> >> Tᵢ := iΔt for i = 1,...,u >> >> It follows: >> >> Tᵤ = t = the final time when we are interested to know the simulated >> values of >> S (in your test you use t = 1). >> >> QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} incrementally >> through a total number of u steps by adding the quantity ΔS := μSΔt + >> σSΔw on >> the previous member of the sequence, with Δw := random number drawn from >> N~(0,Δt¹ᐟ²). >> >> The previous member of S(T₁) is assumed the T₀ = 0 and the initial >> value S(T₀) >> = S(0) is assumed deterministic and known. >> >> So, starting with the given S(T₀), QuantLib calculates the S(T₁) >> using the >> formula: >> >> S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw >> >> Then given the S(T₁) as calculated above, it proceeds to calculate >> S(T₂) as: >> >> S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw >> >> For simplicity I have used the same notation for the Δw in the equations >> above, but they are different. >> >> After u steps, the final value S(t) = S(Tᵤ) is calculated as: >> >> S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw >> >> All equations above hold only in an approximate sense when Δt is finite. >> >> They are accurate only in a limiting sense as Δt → 0. >> >> It follows that the value S(t) will be calculated correctly only if >> Δt → 0, >> which is equivalent to u → ∞ >> >> If you use u = 1 (i.e. one time step), you will still match the >> 68.27% rule >> when the 1SD interval is defined in accordance with a normal >> distribution >> because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal >> distribution and if you build a histogram of all simulated values you >> will >> notice that S(T₁) is distributed normally, which is clearly not >> consistent >> with the SDE that predicts lognormal distribution. >> >> The correct procedure is of course to use a very large u and define >> the 1SD >> interval in accordance with a lognormal distribution, as I described >> in my >> previous email. >> >> Then you will get convergence to the 68.27% rule. >> >> By the way, this is how my histogram of S(t) (for t = 1) simulated >> values >> looks like when I used 100,000 time steps and 100,000 simulations. It is >> clearly lognormal in shape. I am sure, if you do the same when time >> steps = 1, >> you will get a normal distribution for S(t), which would prove that >> time steps >> = 1 does not reproduce the theoretically expected values for S(t). >> >> Ioannis >> >> On 9/11/2023 10:19 AM, U.Mutlu wrote: >>> Dear Ioannis Rigopoulos, I now analysed your table and >>> must say that you seem to have a much different expectation to, >>> and understanding of that table. >>> >>> The values in your table for timeSteps > 1 can definitely >>> not be correct, b/c why do you think it has to be 68.27 for all >>> timeSteps? >>> This is illogical b/c the more timeSteps the more the hitrate% >>> must be, like in my table. It saturates at about 84.93% >>> >>> Sorry, I have to retract my initial enthusiastic judment >>> of your work below, as it's fundamentally flawed & wrong, >>> and I should have been more careful when I quickly & briefly >>> inspected your posting in the early morning here. >>> >>> Since your results in the table for timeSteps > 1 are wrong, >>> then your formulas must be wrong too, so I had to stop my analysis >>> right there. >>> >>> Please just answer why you think it has to be 68.27% for all timeSteps? >>> >>> >>> U.Mutlu wrote on 09/11/23 08:36: >>>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you >>>> very much. >>>> I need some time to study and verify the reults, but it looks very >>>> good >>>> to solve the confusion and mystery I had experienced when I started >>>> this >>>> research some weeks ago. >>>> I'll try to replicate your results and will report later in detail >>>> here. >>>> >>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>> > On a different topic, I do not understand your statement "/I've >>>> never seen >>>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m >>>> +1SD, where >>>> >>>> Sorry for the consfusion, it just means the stock price (Sx after >>>> time t, >>>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper >>>> boundary at >>>> +1SD, >>>> ie. your "min of S(t)" and "max of S(t)" in your other posting. >>>> In future I'll use a more common term like yours, sorry. >>>> >>>> Thx >>>> >>>> >>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t >>>>> >>>>> Are you sure, you also adjusted the boundaries corresponding to >>>>> plus and >>>>> minus >>>>> 1 SD? >>>>> >>>>> I ran an Excel spreadsheet by modifying a spreadsheet I had >>>>> published several >>>>> years ago at >>>>> https://blog.deriscope.com/index.php/en/excel-value-at-risk that >>>>> uses the GeometricBrownianMotionProcess , albeit with the >>>>> PseudoRandom::rsg_type generator. >>>>> >>>>> At any case, I got convergence to 68.27% as I increased the number >>>>> of time >>>>> steps and keeping the number of simulations fixed to 100,000, >>>>> provided I used >>>>> the correct (Ito-adjusted) means, as below: >>>>> >>>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and >>>>> 134.99, just >>>>> like those you use. >>>>> >>>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by >>>>> shifting the >>>>> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. >>>>> >>>>> The simulation results (percentage of values being within the min >>>>> and max >>>>> boundaries) are below: >>>>> >>>>> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 >>>>> time >>>>> steps, >>>>> which is close to the expected 68.27%. >>>>> >>>>> The No Ito-adjusted percentage seems instead to converge to a >>>>> different value >>>>> around 67.73% >>>>> >>>>> On a different topic, I do not understand your statement "/I've >>>>> never seen >>>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>> >>>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, >>>>> where m is >>>>> the mean of a normally distributed variable. This means that one >>>>> needs to >>>>> calculate m in order to do the mentioned computation. >>>>> >>>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then >>>>> ln(S(t)) >>>>> at any time t will be normally distributed with mean m = ln(S(0)) >>>>> + mu*t - >>>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito >>>>> lemma. >>>>> >>>>> Do you have any reason to believe that the Ito lemma is somehow >>>>> not valid so >>>>> that m = ln(S(0)) + mu*t? >>>>> >>>>> Ioannis >>>>> >>>>> On 9/10/2023 3:00 PM, U.Mutlu wrote: >>>>> >>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>>>>> >>>>>> Thank you for your analysis. >>>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>>>>> though in our example with t=1 it doesn't make any difference. >>>>>> It gives for timeSteps=1 this result: >>>>>> >>>>>> $ ./Testing_GBM_of_QuantLib.exe 1 >>>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : >>>>>> -1SD=70.822035 >>>>>> +1SD=129.046162 >>>>>> cHit=663067/1000000(66.307%) >>>>>> >>>>>> This is far from the expected value of 68.2689% for timeStep=1. >>>>>> So, there must be a bug somewhere. IMO it's our suspect friend >>>>>> Ito :-) >>>>>> >>>>>> Cf. also the following showing similar results for GBM, >>>>>> which then was compared to the standard lognormal calculation: >>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get >>>>>> computed with >>>>>> Ito's lemma applied. >>>>>> >>>>>> >>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as >>>>>>> I wrote - >>>>>>> does >>>>>>> not seem to be used in other areas of the standard part of >>>>>>> QuantLib. >>>>>>> >>>>>>> Your numerical results below ... >>>>>>> >>>>>>> indicate that the simulated values for ln(S(t)) produced with >>>>>>> timeSteps = 1 >>>>>>> are very likely normally distributed with mean (I use your >>>>>>> notation) >>>>>>> ln(S(0)) >>>>>>> + mu*t and standard deviation sigma * sqrt(t) >>>>>>> >>>>>>> This result _*is consistent*_ with: >>>>>>> >>>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>>>>> >>>>>>> _*and*_ >>>>>>> >>>>>>> b) the fact that in QuantLib the PathGenerator.next().value >>>>>>> returns the >>>>>>> result >>>>>>> of the SDE expression _*without *_applying the ITO correction >>>>>>> associated >>>>>>> with >>>>>>> the fact that dt is not infinitesimal. >>>>>>> >>>>>>> The b) is also responsible for the lack of convergence in your >>>>>>> output >>>>>>> towards >>>>>>> the theoretical target of 68.27% >>>>>>> >>>>>>> You would get the correct convergence if you modified your code >>>>>>> by using >>>>>>> the >>>>>>> expressions below: >>>>>>> >>>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma >>>>>>> * sqrt(t) * >>>>>>> -1.0); // Sx at -1SD >>>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma >>>>>>> * sqrt(t) * >>>>>>> 1.0); // Sx at +1SD >>>>>>> >>>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, >>>>>>> sigma) (as >>>>>>> Peter also pointed out) >>>>>>> >>>>>>> Again, due to b) one can produce the correct simulated values of >>>>>>> a GBM >>>>>>> diffused quantity (such as a stock price in the GBM model) by >>>>>>> using N time >>>>>>> steps, with N very large. Using N = 1 (like in your example), >>>>>>> the simulated >>>>>>> values will still be lognormally distributed (whence your good >>>>>>> result with >>>>>>> N = >>>>>>> 1), but will be centered at a wrong mean and thus will _*fail *_to >>>>>>> represent >>>>>>> the correct values expected by the GBM SDE. >>>>>>> >>>>>>> Ioannis >>>>>>> >>>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>>>>> As said in other posting here, after fixing the test program >>>>>>>> by skipping the first item (the initial price) in the >>>>>>>> generated sample path, BoxMullerGaussianRng now passes the said >>>>>>>> test. >>>>>>>> >>>>>>>> The bugfixed test code and the new test results can be found here: >>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> That important fact that the generated sample path contains >>>>>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>>>>> to be documented in the library doc. >>>>>>>> For example on this API doc page one normally would expect >>>>>>>> to find such an important information, but it's missing: >>>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>>>>> >>>>>>>> >>>>>>>> If you or someone else can change/extend the test program by using >>>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would >>>>>>>> be happy >>>>>>>> to hear about it. Thx. >>>>>>>> >>>>>>>> >>>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>>>>> If you search within the QuantLib code for >>>>>>>>> BoxMullerGaussianRng, you will >>>>>>>>> see >>>>>>>>> it is used only in the experimental folder. It is therefore not >>>>>>>>> surprising if >>>>>>>>> it doesn't produce the expected results. >>>>>>>>> >>>>>>>>> I use myself the MultiPathGenerator with >>>>>>>>> PseudoRandom::rsg_type, which is >>>>>>>>> used >>>>>>>>> extensively in other areas of QuantLib. >>>>>>>>> >>>>>>>>> This type expands to InverseCumulativeRsg< >>>>>>>>> RandomSequenceGenerator< >>>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and >>>>>>>>> gives me good >>>>>>>>> results. >>>>>>>>> >>>>>>>>> Ioannis Rigopoulos, founder of deriscope.com > > -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 14:19:17
|
Ioannis Rigopoulos wrote on 09/11/23 13:43: > Could you please send us a link to a file that contains the simulated values > produced by your program when you use one time step? > > Since you use 1,000,000 simulations, the file should contain 1,000,000 numbers. On the following web page you can find download links to even 2 such GBM data sets I just have created and uploaded to a server: https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861992 " 2 GBM generated test files in CSV format (with header line; see there), with 1,000,000 generated stock prices (S0=100, IV=30, t=1y, timeSteps=1, r=0, q=0), each about 29 MB uncompressed (about 12 MB zip-compressed). This GBM generator does not use Ito's lemma. ... " > On 9/11/2023 1:33 PM, U.Mutlu wrote: >> Ioannis Rigopoulos wrote on 09/11/23 11:47: >> > By the way, this is how my histogram of S(t) (for t = 1) simulated values >> > looks like when I used 100,000 time steps and 100,000 simulations. It is >> > clearly lognormal in shape. >> >> Yes, it's lognormally. My values are lognormlly as well, >> You just restate the obvious fact. >> >> > I am sure, if you do the same when time steps = 1, >> > you will get a normal distribution for S(t), which would prove >> > that time steps = 1 does not reproduce the theoretically expected >> > values for S(t). >> >> Not true, simply from the fact that the set boundary is lognormally: >> -1SD: 74.08182206817179 >> +1SD: 134.9858807576003 >> Ie. it's not symmetric (Gaussian, normally distributed) around the >> initial mean (100), but is lognormally (and therefore unsymmetric). >> >> FYI: in my other test program (that does not use QuantLib) >> I get in such simulations exact results almost identical >> to the expected theoretical results! So, your math and arguing is wrong, IMO. >> >> You are too much fixated on timeSteps=1. Maybe you have >> a much different understanding of what timeSteps really >> means in practice. It's _not_ something like the number >> of simulations that you simply can change and expect >> to get the same result; no, it has a clearly defined context >> and gives different results for different timeSteps. >> >> And you haven't answered the question why you think it has >> to be 68.27% for all timeSteps (which IMO is wrong). >> >> For those interested, here are more details about the above said: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >> >> >> And, I today did a brief summary of these ongoing discussions on GBM here: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919 >> >> >> >> Ioannis Rigopoulos wrote on 09/11/23 11:47: >>> Our SDE for the evolution of S is: >>> >>> dS = μSdt + σSdw >>> >>> In your code you use: >>> >>> >>> PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> >>> gbmPathGenerator(gbm, t, timeSteps, gsg, false); >>> const Path& samplePath = gbmPathGenerator.next().value; >>> >>> The PathGenerator::next().value is designed to return an array of numbers >>> {S(T₁) , S(T₂) , ... , S(Tᵤ)} >>> >>> where >>> >>> u := timeSteps >>> >>> Δt := t/u >>> >>> Tᵢ := iΔt for i = 1,...,u >>> >>> It follows: >>> >>> Tᵤ = t = the final time when we are interested to know the simulated values of >>> S (in your test you use t = 1). >>> >>> QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} incrementally >>> through a total number of u steps by adding the quantity ΔS := μSΔt + σSΔw on >>> the previous member of the sequence, with Δw := random number drawn from >>> N~(0,Δt¹ᐟ²). >>> >>> The previous member of S(T₁) is assumed the T₀ = 0 and the initial value S(T₀) >>> = S(0) is assumed deterministic and known. >>> >>> So, starting with the given S(T₀), QuantLib calculates the S(T₁) using the >>> formula: >>> >>> S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw >>> >>> Then given the S(T₁) as calculated above, it proceeds to calculate S(T₂) as: >>> >>> S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw >>> >>> For simplicity I have used the same notation for the Δw in the equations >>> above, but they are different. >>> >>> After u steps, the final value S(t) = S(Tᵤ) is calculated as: >>> >>> S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw >>> >>> All equations above hold only in an approximate sense when Δt is finite. >>> >>> They are accurate only in a limiting sense as Δt → 0. >>> >>> It follows that the value S(t) will be calculated correctly only if Δt → 0, >>> which is equivalent to u → ∞ >>> >>> If you use u = 1 (i.e. one time step), you will still match the 68.27% rule >>> when the 1SD interval is defined in accordance with a normal distribution >>> because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal >>> distribution and if you build a histogram of all simulated values you will >>> notice that S(T₁) is distributed normally, which is clearly not consistent >>> with the SDE that predicts lognormal distribution. >>> >>> The correct procedure is of course to use a very large u and define the 1SD >>> interval in accordance with a lognormal distribution, as I described in my >>> previous email. >>> >>> Then you will get convergence to the 68.27% rule. >>> >>> By the way, this is how my histogram of S(t) (for t = 1) simulated values >>> looks like when I used 100,000 time steps and 100,000 simulations. It is >>> clearly lognormal in shape. I am sure, if you do the same when time steps = 1, >>> you will get a normal distribution for S(t), which would prove that time steps >>> = 1 does not reproduce the theoretically expected values for S(t). >>> >>> Ioannis >>> >>> On 9/11/2023 10:19 AM, U.Mutlu wrote: >>>> Dear Ioannis Rigopoulos, I now analysed your table and >>>> must say that you seem to have a much different expectation to, >>>> and understanding of that table. >>>> >>>> The values in your table for timeSteps > 1 can definitely >>>> not be correct, b/c why do you think it has to be 68.27 for all timeSteps? >>>> This is illogical b/c the more timeSteps the more the hitrate% >>>> must be, like in my table. It saturates at about 84.93% >>>> >>>> Sorry, I have to retract my initial enthusiastic judment >>>> of your work below, as it's fundamentally flawed & wrong, >>>> and I should have been more careful when I quickly & briefly >>>> inspected your posting in the early morning here. >>>> >>>> Since your results in the table for timeSteps > 1 are wrong, >>>> then your formulas must be wrong too, so I had to stop my analysis right >>>> there. >>>> >>>> Please just answer why you think it has to be 68.27% for all timeSteps? >>>> >>>> >>>> U.Mutlu wrote on 09/11/23 08:36: >>>>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you very much. >>>>> I need some time to study and verify the reults, but it looks very good >>>>> to solve the confusion and mystery I had experienced when I started this >>>>> research some weeks ago. >>>>> I'll try to replicate your results and will report later in detail here. >>>>> >>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>> > On a different topic, I do not understand your statement "/I've never >>>>> seen >>>>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where >>>>> >>>>> Sorry for the consfusion, it just means the stock price (Sx after time t, >>>>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper >>>>> boundary at >>>>> +1SD, >>>>> ie. your "min of S(t)" and "max of S(t)" in your other posting. >>>>> In future I'll use a more common term like yours, sorry. >>>>> >>>>> Thx >>>>> >>>>> >>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t >>>>>> >>>>>> Are you sure, you also adjusted the boundaries corresponding to plus and >>>>>> minus >>>>>> 1 SD? >>>>>> >>>>>> I ran an Excel spreadsheet by modifying a spreadsheet I had published >>>>>> several >>>>>> years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk >>>>>> that >>>>>> uses the GeometricBrownianMotionProcess , albeit with the >>>>>> PseudoRandom::rsg_type generator. >>>>>> >>>>>> At any case, I got convergence to 68.27% as I increased the number of time >>>>>> steps and keeping the number of simulations fixed to 100,000, provided I >>>>>> used >>>>>> the correct (Ito-adjusted) means, as below: >>>>>> >>>>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, >>>>>> just >>>>>> like those you use. >>>>>> >>>>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the >>>>>> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. >>>>>> >>>>>> The simulation results (percentage of values being within the min and max >>>>>> boundaries) are below: >>>>>> >>>>>> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time >>>>>> steps, >>>>>> which is close to the expected 68.27%. >>>>>> >>>>>> The No Ito-adjusted percentage seems instead to converge to a different >>>>>> value >>>>>> around 67.73% >>>>>> >>>>>> On a different topic, I do not understand your statement "/I've never seen >>>>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>>> >>>>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where >>>>>> m is >>>>>> the mean of a normally distributed variable. This means that one needs to >>>>>> calculate m in order to do the mentioned computation. >>>>>> >>>>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then >>>>>> ln(S(t)) >>>>>> at any time t will be normally distributed with mean m = ln(S(0)) + mu*t - >>>>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. >>>>>> >>>>>> Do you have any reason to believe that the Ito lemma is somehow not >>>>>> valid so >>>>>> that m = ln(S(0)) + mu*t? >>>>>> >>>>>> Ioannis >>>>>> >>>>>> On 9/10/2023 3:00 PM, U.Mutlu wrote: >>>>>> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>>>>>> >>>>>>> Thank you for your analysis. >>>>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>>>>>> though in our example with t=1 it doesn't make any difference. >>>>>>> It gives for timeSteps=1 this result: >>>>>>> >>>>>>> $ ./Testing_GBM_of_QuantLib.exe 1 >>>>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>>>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 >>>>>>> +1SD=129.046162 >>>>>>> cHit=663067/1000000(66.307%) >>>>>>> >>>>>>> This is far from the expected value of 68.2689% for timeStep=1. >>>>>>> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >>>>>>> >>>>>>> Cf. also the following showing similar results for GBM, >>>>>>> which then was compared to the standard lognormal calculation: >>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed with >>>>>>> Ito's lemma applied. >>>>>>> >>>>>>> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as I wrote - >>>>>>>> does >>>>>>>> not seem to be used in other areas of the standard part of QuantLib. >>>>>>>> >>>>>>>> Your numerical results below ... >>>>>>>> >>>>>>>> indicate that the simulated values for ln(S(t)) produced with >>>>>>>> timeSteps = 1 >>>>>>>> are very likely normally distributed with mean (I use your notation) >>>>>>>> ln(S(0)) >>>>>>>> + mu*t and standard deviation sigma * sqrt(t) >>>>>>>> >>>>>>>> This result _*is consistent*_ with: >>>>>>>> >>>>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>>>>>> >>>>>>>> _*and*_ >>>>>>>> >>>>>>>> b) the fact that in QuantLib the PathGenerator.next().value returns the >>>>>>>> result >>>>>>>> of the SDE expression _*without *_applying the ITO correction associated >>>>>>>> with >>>>>>>> the fact that dt is not infinitesimal. >>>>>>>> >>>>>>>> The b) is also responsible for the lack of convergence in your output >>>>>>>> towards >>>>>>>> the theoretical target of 68.27% >>>>>>>> >>>>>>>> You would get the correct convergence if you modified your code by using >>>>>>>> the >>>>>>>> expressions below: >>>>>>>> >>>>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>> sqrt(t) * >>>>>>>> -1.0); // Sx at -1SD >>>>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>> sqrt(t) * >>>>>>>> 1.0); // Sx at +1SD >>>>>>>> >>>>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) (as >>>>>>>> Peter also pointed out) >>>>>>>> >>>>>>>> Again, due to b) one can produce the correct simulated values of a GBM >>>>>>>> diffused quantity (such as a stock price in the GBM model) by using N >>>>>>>> time >>>>>>>> steps, with N very large. Using N = 1 (like in your example), the >>>>>>>> simulated >>>>>>>> values will still be lognormally distributed (whence your good result >>>>>>>> with >>>>>>>> N = >>>>>>>> 1), but will be centered at a wrong mean and thus will _*fail *_to >>>>>>>> represent >>>>>>>> the correct values expected by the GBM SDE. >>>>>>>> >>>>>>>> Ioannis >>>>>>>> >>>>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>>>>>> As said in other posting here, after fixing the test program >>>>>>>>> by skipping the first item (the initial price) in the >>>>>>>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>>>>>>> >>>>>>>>> The bugfixed test code and the new test results can be found here: >>>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> That important fact that the generated sample path contains >>>>>>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>>>>>> to be documented in the library doc. >>>>>>>>> For example on this API doc page one normally would expect >>>>>>>>> to find such an important information, but it's missing: >>>>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>>>>>> >>>>>>>>> >>>>>>>>> If you or someone else can change/extend the test program by using >>>>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>>>>>>> to hear about it. Thx. >>>>>>>>> >>>>>>>>> >>>>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>>>>>> If you search within the QuantLib code for BoxMullerGaussianRng, you >>>>>>>>>> will >>>>>>>>>> see >>>>>>>>>> it is used only in the experimental folder. It is therefore not >>>>>>>>>> surprising if >>>>>>>>>> it doesn't produce the expected results. >>>>>>>>>> >>>>>>>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, >>>>>>>>>> which is >>>>>>>>>> used >>>>>>>>>> extensively in other areas of QuantLib. >>>>>>>>>> >>>>>>>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me >>>>>>>>>> good >>>>>>>>>> results. >>>>>>>>>> >>>>>>>>>> Ioannis Rigopoulos, founder of deriscope.com |
|
From: Ioannis R. <qua...@de...> - 2023-09-11 15:29:12
|
Thank you U.Mutlu. The histograms indicate that the numbers represent indeed lognormal distributions. The average of all numbers in file GBM_ts1_A amounts to *104.5935229* The average of all numbers in file GBM_ts1_B amounts to *104.6047793* Both data sets exhibit the correct log standard deviation, very close to *0.3*. *Since you have assumed a drift of 0, the SDE represents a martingale and therefore the average stock price at t=1 owes to be very close to 100.* *This proves that the simulated numbers cannot represent the GBM diffusion of a stock price starting at 100 and having drift = 0.* Furthermore, the quantity S(0)exp(0.5*sigma^2*t) equals *104.602786*, which is extremely close to the above averages and indicates that the Ito term 0.5*sigma^2*t is embedded in your results, albeit in a wrong fashion, since - as mentioned above - it should not affect the average of S(t). It seems that you have set up your code in such a way that ln(S) is driftless and therefore S acquires a positive Ito drift of 0.5*sigma^2. It also seems that the number of time steps are much greater than 1, otherwise the distribution would appear as normal. At any case, since I cannot run your own code and verify the results, I will finish my feedback with this email. Ioannis On 9/11/2023 4:19 PM, U.Mutlu wrote: > Ioannis Rigopoulos wrote on 09/11/23 13:43: >> Could you please send us a link to a file that contains the simulated >> values >> produced by your program when you use one time step? >> >> Since you use 1,000,000 simulations, the file should contain >> 1,000,000 numbers. > > On the following web page you can find download links to even 2 such > GBM data sets I just have created and uploaded to a server: > https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861992 > > > " > 2 GBM generated test files in CSV format (with header line; see there), > with 1,000,000 generated stock prices (S0=100, IV=30, t=1y, > timeSteps=1, r=0, q=0), > each about 29 MB uncompressed (about 12 MB zip-compressed). > > This GBM generator does not use Ito's lemma. > > ... > " > > >> On 9/11/2023 1:33 PM, U.Mutlu wrote: >>> Ioannis Rigopoulos wrote on 09/11/23 11:47: >>> > By the way, this is how my histogram of S(t) (for t = 1) simulated >>> values >>> > looks like when I used 100,000 time steps and 100,000 simulations. >>> It is >>> > clearly lognormal in shape. >>> >>> Yes, it's lognormally. My values are lognormlly as well, >>> You just restate the obvious fact. >>> >>> > I am sure, if you do the same when time steps = 1, >>> > you will get a normal distribution for S(t), which would prove >>> > that time steps = 1 does not reproduce the theoretically expected >>> > values for S(t). >>> >>> Not true, simply from the fact that the set boundary is lognormally: >>> -1SD: 74.08182206817179 >>> +1SD: 134.9858807576003 >>> Ie. it's not symmetric (Gaussian, normally distributed) around the >>> initial mean (100), but is lognormally (and therefore unsymmetric). >>> >>> FYI: in my other test program (that does not use QuantLib) >>> I get in such simulations exact results almost identical >>> to the expected theoretical results! So, your math and arguing is >>> wrong, IMO. >>> >>> You are too much fixated on timeSteps=1. Maybe you have >>> a much different understanding of what timeSteps really >>> means in practice. It's _not_ something like the number >>> of simulations that you simply can change and expect >>> to get the same result; no, it has a clearly defined context >>> and gives different results for different timeSteps. >>> >>> And you haven't answered the question why you think it has >>> to be 68.27% for all timeSteps (which IMO is wrong). >>> >>> For those interested, here are more details about the above said: >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>> >>> >>> >>> And, I today did a brief summary of these ongoing discussions on GBM >>> here: >>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919 >>> >>> >>> >>> >>> Ioannis Rigopoulos wrote on 09/11/23 11:47: >>>> Our SDE for the evolution of S is: >>>> >>>> dS = μSdt + σSdw >>>> >>>> In your code you use: >>>> >>>> >>>> PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> >>>> gbmPathGenerator(gbm, t, timeSteps, gsg, false); >>>> const Path& samplePath = gbmPathGenerator.next().value; >>>> >>>> The PathGenerator::next().value is designed to return an array of >>>> numbers >>>> {S(T₁) , S(T₂) , ... , S(Tᵤ)} >>>> >>>> where >>>> >>>> u := timeSteps >>>> >>>> Δt := t/u >>>> >>>> Tᵢ := iΔt for i = 1,...,u >>>> >>>> It follows: >>>> >>>> Tᵤ = t = the final time when we are interested to know the >>>> simulated values of >>>> S (in your test you use t = 1). >>>> >>>> QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} >>>> incrementally >>>> through a total number of u steps by adding the quantity ΔS := μSΔt >>>> + σSΔw on >>>> the previous member of the sequence, with Δw := random number drawn >>>> from >>>> N~(0,Δt¹ᐟ²). >>>> >>>> The previous member of S(T₁) is assumed the T₀ = 0 and the initial >>>> value S(T₀) >>>> = S(0) is assumed deterministic and known. >>>> >>>> So, starting with the given S(T₀), QuantLib calculates the S(T₁) >>>> using the >>>> formula: >>>> >>>> S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw >>>> >>>> Then given the S(T₁) as calculated above, it proceeds to calculate >>>> S(T₂) as: >>>> >>>> S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw >>>> >>>> For simplicity I have used the same notation for the Δw in the >>>> equations >>>> above, but they are different. >>>> >>>> After u steps, the final value S(t) = S(Tᵤ) is calculated as: >>>> >>>> S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw >>>> >>>> All equations above hold only in an approximate sense when Δt is >>>> finite. >>>> >>>> They are accurate only in a limiting sense as Δt → 0. >>>> >>>> It follows that the value S(t) will be calculated correctly only if >>>> Δt → 0, >>>> which is equivalent to u → ∞ >>>> >>>> If you use u = 1 (i.e. one time step), you will still match the >>>> 68.27% rule >>>> when the 1SD interval is defined in accordance with a normal >>>> distribution >>>> because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal >>>> distribution and if you build a histogram of all simulated values >>>> you will >>>> notice that S(T₁) is distributed normally, which is clearly not >>>> consistent >>>> with the SDE that predicts lognormal distribution. >>>> >>>> The correct procedure is of course to use a very large u and define >>>> the 1SD >>>> interval in accordance with a lognormal distribution, as I >>>> described in my >>>> previous email. >>>> >>>> Then you will get convergence to the 68.27% rule. >>>> >>>> By the way, this is how my histogram of S(t) (for t = 1) simulated >>>> values >>>> looks like when I used 100,000 time steps and 100,000 simulations. >>>> It is >>>> clearly lognormal in shape. I am sure, if you do the same when time >>>> steps = 1, >>>> you will get a normal distribution for S(t), which would prove that >>>> time steps >>>> = 1 does not reproduce the theoretically expected values for S(t). >>>> >>>> Ioannis >>>> >>>> On 9/11/2023 10:19 AM, U.Mutlu wrote: >>>>> Dear Ioannis Rigopoulos, I now analysed your table and >>>>> must say that you seem to have a much different expectation to, >>>>> and understanding of that table. >>>>> >>>>> The values in your table for timeSteps > 1 can definitely >>>>> not be correct, b/c why do you think it has to be 68.27 for all >>>>> timeSteps? >>>>> This is illogical b/c the more timeSteps the more the hitrate% >>>>> must be, like in my table. It saturates at about 84.93% >>>>> >>>>> Sorry, I have to retract my initial enthusiastic judment >>>>> of your work below, as it's fundamentally flawed & wrong, >>>>> and I should have been more careful when I quickly & briefly >>>>> inspected your posting in the early morning here. >>>>> >>>>> Since your results in the table for timeSteps > 1 are wrong, >>>>> then your formulas must be wrong too, so I had to stop my analysis >>>>> right >>>>> there. >>>>> >>>>> Please just answer why you think it has to be 68.27% for all >>>>> timeSteps? >>>>> >>>>> >>>>> U.Mutlu wrote on 09/11/23 08:36: >>>>>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you >>>>>> very much. >>>>>> I need some time to study and verify the reults, but it looks >>>>>> very good >>>>>> to solve the confusion and mystery I had experienced when I >>>>>> started this >>>>>> research some weeks ago. >>>>>> I'll try to replicate your results and will report later in >>>>>> detail here. >>>>>> >>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>>> > On a different topic, I do not understand your statement >>>>>> "/I've never >>>>>> seen >>>>>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m >>>>>> +1SD, where >>>>>> >>>>>> Sorry for the consfusion, it just means the stock price (Sx after >>>>>> time t, >>>>>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper >>>>>> boundary at >>>>>> +1SD, >>>>>> ie. your "min of S(t)" and "max of S(t)" in your other posting. >>>>>> In future I'll use a more common term like yours, sorry. >>>>>> >>>>>> Thx >>>>>> >>>>>> >>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t >>>>>>> >>>>>>> Are you sure, you also adjusted the boundaries corresponding to >>>>>>> plus and >>>>>>> minus >>>>>>> 1 SD? >>>>>>> >>>>>>> I ran an Excel spreadsheet by modifying a spreadsheet I had >>>>>>> published >>>>>>> several >>>>>>> years ago at >>>>>>> https://blog.deriscope.com/index.php/en/excel-value-at-risk >>>>>>> that >>>>>>> uses the GeometricBrownianMotionProcess , albeit with the >>>>>>> PseudoRandom::rsg_type generator. >>>>>>> >>>>>>> At any case, I got convergence to 68.27% as I increased the >>>>>>> number of time >>>>>>> steps and keeping the number of simulations fixed to 100,000, >>>>>>> provided I >>>>>>> used >>>>>>> the correct (Ito-adjusted) means, as below: >>>>>>> >>>>>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and >>>>>>> 134.99, >>>>>>> just >>>>>>> like those you use. >>>>>>> >>>>>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by >>>>>>> shifting the >>>>>>> no-adjusted boundaries - in log terms - to the left by >>>>>>> 0.5*sigma^2*t. >>>>>>> >>>>>>> The simulation results (percentage of values being within the >>>>>>> min and max >>>>>>> boundaries) are below: >>>>>>> >>>>>>> As you see, the Ito-adjusted percentage reaches 68.30% for >>>>>>> 100,000 time >>>>>>> steps, >>>>>>> which is close to the expected 68.27%. >>>>>>> >>>>>>> The No Ito-adjusted percentage seems instead to converge to a >>>>>>> different >>>>>>> value >>>>>>> around 67.73% >>>>>>> >>>>>>> On a different topic, I do not understand your statement "/I've >>>>>>> never seen >>>>>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>>>> >>>>>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m >>>>>>> +1SD, where >>>>>>> m is >>>>>>> the mean of a normally distributed variable. This means that one >>>>>>> needs to >>>>>>> calculate m in order to do the mentioned computation. >>>>>>> >>>>>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, >>>>>>> then >>>>>>> ln(S(t)) >>>>>>> at any time t will be normally distributed with mean m = >>>>>>> ln(S(0)) + mu*t - >>>>>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito >>>>>>> lemma. >>>>>>> >>>>>>> Do you have any reason to believe that the Ito lemma is somehow not >>>>>>> valid so >>>>>>> that m = ln(S(0)) + mu*t? >>>>>>> >>>>>>> Ioannis >>>>>>> >>>>>>> On 9/10/2023 3:00 PM, U.Mutlu wrote: >>>>>>> >>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>>>>>>> >>>>>>>> Thank you for your analysis. >>>>>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>>>>>>> though in our example with t=1 it doesn't make any difference. >>>>>>>> It gives for timeSteps=1 this result: >>>>>>>> >>>>>>>> $ ./Testing_GBM_of_QuantLib.exe 1 >>>>>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>>>>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : >>>>>>>> -1SD=70.822035 >>>>>>>> +1SD=129.046162 >>>>>>>> cHit=663067/1000000(66.307%) >>>>>>>> >>>>>>>> This is far from the expected value of 68.2689% for timeStep=1. >>>>>>>> So, there must be a bug somewhere. IMO it's our suspect friend >>>>>>>> Ito :-) >>>>>>>> >>>>>>>> Cf. also the following showing similar results for GBM, >>>>>>>> which then was compared to the standard lognormal calculation: >>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get >>>>>>>> computed with >>>>>>>> Ito's lemma applied. >>>>>>>> >>>>>>>> >>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>>>> Thank you for testing the BoxMullerGaussianRng code, which - >>>>>>>>> as I wrote - >>>>>>>>> does >>>>>>>>> not seem to be used in other areas of the standard part of >>>>>>>>> QuantLib. >>>>>>>>> >>>>>>>>> Your numerical results below ... >>>>>>>>> >>>>>>>>> indicate that the simulated values for ln(S(t)) produced with >>>>>>>>> timeSteps = 1 >>>>>>>>> are very likely normally distributed with mean (I use your >>>>>>>>> notation) >>>>>>>>> ln(S(0)) >>>>>>>>> + mu*t and standard deviation sigma * sqrt(t) >>>>>>>>> >>>>>>>>> This result _*is consistent*_ with: >>>>>>>>> >>>>>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>>>>>>> >>>>>>>>> _*and*_ >>>>>>>>> >>>>>>>>> b) the fact that in QuantLib the PathGenerator.next().value >>>>>>>>> returns the >>>>>>>>> result >>>>>>>>> of the SDE expression _*without *_applying the ITO correction >>>>>>>>> associated >>>>>>>>> with >>>>>>>>> the fact that dt is not infinitesimal. >>>>>>>>> >>>>>>>>> The b) is also responsible for the lack of convergence in your >>>>>>>>> output >>>>>>>>> towards >>>>>>>>> the theoretical target of 68.27% >>>>>>>>> >>>>>>>>> You would get the correct convergence if you modified your >>>>>>>>> code by using >>>>>>>>> the >>>>>>>>> expressions below: >>>>>>>>> >>>>>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>>> sqrt(t) * >>>>>>>>> -1.0); // Sx at -1SD >>>>>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>>> sqrt(t) * >>>>>>>>> 1.0); // Sx at +1SD >>>>>>>>> >>>>>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, >>>>>>>>> sigma) (as >>>>>>>>> Peter also pointed out) >>>>>>>>> >>>>>>>>> Again, due to b) one can produce the correct simulated values >>>>>>>>> of a GBM >>>>>>>>> diffused quantity (such as a stock price in the GBM model) by >>>>>>>>> using N >>>>>>>>> time >>>>>>>>> steps, with N very large. Using N = 1 (like in your example), the >>>>>>>>> simulated >>>>>>>>> values will still be lognormally distributed (whence your good >>>>>>>>> result >>>>>>>>> with >>>>>>>>> N = >>>>>>>>> 1), but will be centered at a wrong mean and thus will _*fail >>>>>>>>> *_to >>>>>>>>> represent >>>>>>>>> the correct values expected by the GBM SDE. >>>>>>>>> >>>>>>>>> Ioannis >>>>>>>>> >>>>>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>>>>>>> As said in other posting here, after fixing the test program >>>>>>>>>> by skipping the first item (the initial price) in the >>>>>>>>>> generated sample path, BoxMullerGaussianRng now passes the >>>>>>>>>> said test. >>>>>>>>>> >>>>>>>>>> The bugfixed test code and the new test results can be found >>>>>>>>>> here: >>>>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> That important fact that the generated sample path contains >>>>>>>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>>>>>>> to be documented in the library doc. >>>>>>>>>> For example on this API doc page one normally would expect >>>>>>>>>> to find such an important information, but it's missing: >>>>>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> If you or someone else can change/extend the test program by >>>>>>>>>> using >>>>>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would >>>>>>>>>> be happy >>>>>>>>>> to hear about it. Thx. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>>>>>>> If you search within the QuantLib code for >>>>>>>>>>> BoxMullerGaussianRng, you >>>>>>>>>>> will >>>>>>>>>>> see >>>>>>>>>>> it is used only in the experimental folder. It is therefore not >>>>>>>>>>> surprising if >>>>>>>>>>> it doesn't produce the expected results. >>>>>>>>>>> >>>>>>>>>>> I use myself the MultiPathGenerator with >>>>>>>>>>> PseudoRandom::rsg_type, >>>>>>>>>>> which is >>>>>>>>>>> used >>>>>>>>>>> extensively in other areas of QuantLib. >>>>>>>>>>> >>>>>>>>>>> This type expands to InverseCumulativeRsg< >>>>>>>>>>> RandomSequenceGenerator< >>>>>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and >>>>>>>>>>> gives me >>>>>>>>>>> good >>>>>>>>>>> results. >>>>>>>>>>> >>>>>>>>>>> Ioannis Rigopoulos, founder of deriscope.com > > -- This email has been checked for viruses by Avast antivirus software. www.avast.com |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 15:57:05
|
Ioannis Rigopoulos wrote on 09/11/23 17:28: > Thank you U.Mutlu. > > The histograms indicate that the numbers represent indeed lognormal distributions. > > The average of all numbers in file GBM_ts1_A amounts to *104.5935229* > > The average of all numbers in file GBM_ts1_B amounts to *104.6047793* > > Both data sets exhibit the correct log standard deviation, very close to *0.3*. So far so good, and these are the most simplest things, but with the rest of your analysis I don't agree. What do other experts say on the following analysis of Ioannis? : > *Since you have assumed a drift of 0, the SDE represents a martingale and > therefore the average stock price at t=1 owes to be very close to 100.* > > *This proves that the simulated numbers cannot represent the GBM diffusion of > a stock price starting at 100 and having drift = 0.* > > Furthermore, the quantity S(0)exp(0.5*sigma^2*t) equals *104.602786*, which is > extremely close to the above averages and indicates that the Ito term > 0.5*sigma^2*t is embedded in your results, albeit in a wrong fashion, since - > as mentioned above - it should not affect the average of S(t). > > It seems that you have set up your code in such a way that ln(S) is driftless > and therefore S acquires a positive Ito drift of 0.5*sigma^2. > > It also seems that the number of time steps are much greater than 1, otherwise > the distribution would appear as normal. > > At any case, since I cannot run your own code and verify the results, I will > finish my feedback with this email. > > Ioannis Oh oh, that's too much incorrect assumptions & interpretations by you. Nevertheless, thank you for attempting to analyse it. Maybe some other experts too will analyse the data and post their findings here. I claim that Ito's lemma is the culprit and that it's unnecessary both in GBM and BSM (Black-Scholes-Merton option pricing formula). > On 9/11/2023 4:19 PM, U.Mutlu wrote: >> Ioannis Rigopoulos wrote on 09/11/23 13:43: >>> Could you please send us a link to a file that contains the simulated values >>> produced by your program when you use one time step? >>> >>> Since you use 1,000,000 simulations, the file should contain 1,000,000 >>> numbers. >> >> On the following web page you can find download links to even 2 such >> GBM data sets I just have created and uploaded to a server: >> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861992 >> >> >> " >> 2 GBM generated test files in CSV format (with header line; see there), >> with 1,000,000 generated stock prices (S0=100, IV=30, t=1y, timeSteps=1, >> r=0, q=0), >> each about 29 MB uncompressed (about 12 MB zip-compressed). >> >> This GBM generator does not use Ito's lemma. >> >> ... >> " >> >> >>> On 9/11/2023 1:33 PM, U.Mutlu wrote: >>>> Ioannis Rigopoulos wrote on 09/11/23 11:47: >>>> > By the way, this is how my histogram of S(t) (for t = 1) simulated values >>>> > looks like when I used 100,000 time steps and 100,000 simulations. It is >>>> > clearly lognormal in shape. >>>> >>>> Yes, it's lognormally. My values are lognormlly as well, >>>> You just restate the obvious fact. >>>> >>>> > I am sure, if you do the same when time steps = 1, >>>> > you will get a normal distribution for S(t), which would prove >>>> > that time steps = 1 does not reproduce the theoretically expected >>>> > values for S(t). >>>> >>>> Not true, simply from the fact that the set boundary is lognormally: >>>> -1SD: 74.08182206817179 >>>> +1SD: 134.9858807576003 >>>> Ie. it's not symmetric (Gaussian, normally distributed) around the >>>> initial mean (100), but is lognormally (and therefore unsymmetric). >>>> >>>> FYI: in my other test program (that does not use QuantLib) >>>> I get in such simulations exact results almost identical >>>> to the expected theoretical results! So, your math and arguing is wrong, IMO. >>>> >>>> You are too much fixated on timeSteps=1. Maybe you have >>>> a much different understanding of what timeSteps really >>>> means in practice. It's _not_ something like the number >>>> of simulations that you simply can change and expect >>>> to get the same result; no, it has a clearly defined context >>>> and gives different results for different timeSteps. >>>> >>>> And you haven't answered the question why you think it has >>>> to be 68.27% for all timeSteps (which IMO is wrong). >>>> >>>> For those interested, here are more details about the above said: >>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>> >>>> >>>> >>>> And, I today did a brief summary of these ongoing discussions on GBM here: >>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919 >>>> >>>> >>>> >>>> >>>> Ioannis Rigopoulos wrote on 09/11/23 11:47: >>>>> Our SDE for the evolution of S is: >>>>> >>>>> dS = μSdt + σSdw >>>>> >>>>> In your code you use: >>>>> >>>>> >>>>> PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>> >>>>> gbmPathGenerator(gbm, t, timeSteps, gsg, false); >>>>> const Path& samplePath = gbmPathGenerator.next().value; >>>>> >>>>> The PathGenerator::next().value is designed to return an array of numbers >>>>> {S(T₁) , S(T₂) , ... , S(Tᵤ)} >>>>> >>>>> where >>>>> >>>>> u := timeSteps >>>>> >>>>> Δt := t/u >>>>> >>>>> Tᵢ := iΔt for i = 1,...,u >>>>> >>>>> It follows: >>>>> >>>>> Tᵤ = t = the final time when we are interested to know the simulated >>>>> values of >>>>> S (in your test you use t = 1). >>>>> >>>>> QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)} incrementally >>>>> through a total number of u steps by adding the quantity ΔS := μSΔt + >>>>> σSΔw on >>>>> the previous member of the sequence, with Δw := random number drawn from >>>>> N~(0,Δt¹ᐟ²). >>>>> >>>>> The previous member of S(T₁) is assumed the T₀ = 0 and the initial value >>>>> S(T₀) >>>>> = S(0) is assumed deterministic and known. >>>>> >>>>> So, starting with the given S(T₀), QuantLib calculates the S(T₁) using the >>>>> formula: >>>>> >>>>> S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw >>>>> >>>>> Then given the S(T₁) as calculated above, it proceeds to calculate S(T₂) as: >>>>> >>>>> S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw >>>>> >>>>> For simplicity I have used the same notation for the Δw in the equations >>>>> above, but they are different. >>>>> >>>>> After u steps, the final value S(t) = S(Tᵤ) is calculated as: >>>>> >>>>> S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw >>>>> >>>>> All equations above hold only in an approximate sense when Δt is finite. >>>>> >>>>> They are accurate only in a limiting sense as Δt → 0. >>>>> >>>>> It follows that the value S(t) will be calculated correctly only if Δt → 0, >>>>> which is equivalent to u → ∞ >>>>> >>>>> If you use u = 1 (i.e. one time step), you will still match the 68.27% rule >>>>> when the 1SD interval is defined in accordance with a normal distribution >>>>> because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal >>>>> distribution and if you build a histogram of all simulated values you will >>>>> notice that S(T₁) is distributed normally, which is clearly not consistent >>>>> with the SDE that predicts lognormal distribution. >>>>> >>>>> The correct procedure is of course to use a very large u and define the 1SD >>>>> interval in accordance with a lognormal distribution, as I described in my >>>>> previous email. >>>>> >>>>> Then you will get convergence to the 68.27% rule. >>>>> >>>>> By the way, this is how my histogram of S(t) (for t = 1) simulated values >>>>> looks like when I used 100,000 time steps and 100,000 simulations. It is >>>>> clearly lognormal in shape. I am sure, if you do the same when time steps >>>>> = 1, >>>>> you will get a normal distribution for S(t), which would prove that time >>>>> steps >>>>> = 1 does not reproduce the theoretically expected values for S(t). >>>>> >>>>> Ioannis >>>>> >>>>> On 9/11/2023 10:19 AM, U.Mutlu wrote: >>>>>> Dear Ioannis Rigopoulos, I now analysed your table and >>>>>> must say that you seem to have a much different expectation to, >>>>>> and understanding of that table. >>>>>> >>>>>> The values in your table for timeSteps > 1 can definitely >>>>>> not be correct, b/c why do you think it has to be 68.27 for all timeSteps? >>>>>> This is illogical b/c the more timeSteps the more the hitrate% >>>>>> must be, like in my table. It saturates at about 84.93% >>>>>> >>>>>> Sorry, I have to retract my initial enthusiastic judment >>>>>> of your work below, as it's fundamentally flawed & wrong, >>>>>> and I should have been more careful when I quickly & briefly >>>>>> inspected your posting in the early morning here. >>>>>> >>>>>> Since your results in the table for timeSteps > 1 are wrong, >>>>>> then your formulas must be wrong too, so I had to stop my analysis right >>>>>> there. >>>>>> >>>>>> Please just answer why you think it has to be 68.27% for all timeSteps? >>>>>> >>>>>> >>>>>> U.Mutlu wrote on 09/11/23 08:36: >>>>>>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you very >>>>>>> much. >>>>>>> I need some time to study and verify the reults, but it looks very good >>>>>>> to solve the confusion and mystery I had experienced when I started this >>>>>>> research some weeks ago. >>>>>>> I'll try to replicate your results and will report later in detail here. >>>>>>> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>>>> > On a different topic, I do not understand your statement "/I've never >>>>>>> seen >>>>>>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>>>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, >>>>>>> where >>>>>>> >>>>>>> Sorry for the consfusion, it just means the stock price (Sx after time t, >>>>>>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper >>>>>>> boundary at >>>>>>> +1SD, >>>>>>> ie. your "min of S(t)" and "max of S(t)" in your other posting. >>>>>>> In future I'll use a more common term like yours, sorry. >>>>>>> >>>>>>> Thx >>>>>>> >>>>>>> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29: >>>>>>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t >>>>>>>> >>>>>>>> Are you sure, you also adjusted the boundaries corresponding to plus and >>>>>>>> minus >>>>>>>> 1 SD? >>>>>>>> >>>>>>>> I ran an Excel spreadsheet by modifying a spreadsheet I had published >>>>>>>> several >>>>>>>> years ago at https://blog.deriscope.com/index.php/en/excel-value-at-risk >>>>>>>> that >>>>>>>> uses the GeometricBrownianMotionProcess , albeit with the >>>>>>>> PseudoRandom::rsg_type generator. >>>>>>>> >>>>>>>> At any case, I got convergence to 68.27% as I increased the number of >>>>>>>> time >>>>>>>> steps and keeping the number of simulations fixed to 100,000, provided I >>>>>>>> used >>>>>>>> the correct (Ito-adjusted) means, as below: >>>>>>>> >>>>>>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and 134.99, >>>>>>>> just >>>>>>>> like those you use. >>>>>>>> >>>>>>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by shifting the >>>>>>>> no-adjusted boundaries - in log terms - to the left by 0.5*sigma^2*t. >>>>>>>> >>>>>>>> The simulation results (percentage of values being within the min and max >>>>>>>> boundaries) are below: >>>>>>>> >>>>>>>> As you see, the Ito-adjusted percentage reaches 68.30% for 100,000 time >>>>>>>> steps, >>>>>>>> which is close to the expected 68.27%. >>>>>>>> >>>>>>>> The No Ito-adjusted percentage seems instead to converge to a different >>>>>>>> value >>>>>>>> around 67.73% >>>>>>>> >>>>>>>> On a different topic, I do not understand your statement "/I've never >>>>>>>> seen >>>>>>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/". >>>>>>>> >>>>>>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m +1SD, where >>>>>>>> m is >>>>>>>> the mean of a normally distributed variable. This means that one needs to >>>>>>>> calculate m in order to do the mentioned computation. >>>>>>>> >>>>>>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw, then >>>>>>>> ln(S(t)) >>>>>>>> at any time t will be normally distributed with mean m = ln(S(0)) + >>>>>>>> mu*t - >>>>>>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito lemma. >>>>>>>> >>>>>>>> Do you have any reason to believe that the Ito lemma is somehow not >>>>>>>> valid so >>>>>>>> that m = ln(S(0)) + mu*t? >>>>>>>> >>>>>>>> Ioannis >>>>>>>> >>>>>>>> On 9/10/2023 3:00 PM, U.Mutlu wrote: >>>>>>>> >>>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, sigma) >>>>>>>>> >>>>>>>>> Thank you for your analysis. >>>>>>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t", >>>>>>>>> though in our example with t=1 it doesn't make any difference. >>>>>>>>> It gives for timeSteps=1 this result: >>>>>>>>> >>>>>>>>> $ ./Testing_GBM_of_QuantLib.exe 1 >>>>>>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000 mu=0.000000 >>>>>>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 : -1SD=70.822035 >>>>>>>>> +1SD=129.046162 >>>>>>>>> cHit=663067/1000000(66.307%) >>>>>>>>> >>>>>>>>> This is far from the expected value of 68.2689% for timeStep=1. >>>>>>>>> So, there must be a bug somewhere. IMO it's our suspect friend Ito :-) >>>>>>>>> >>>>>>>>> Cf. also the following showing similar results for GBM, >>>>>>>>> which then was compared to the standard lognormal calculation: >>>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548 >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get computed >>>>>>>>> with >>>>>>>>> Ito's lemma applied. >>>>>>>>> >>>>>>>>> >>>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06: >>>>>>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as I >>>>>>>>>> wrote - >>>>>>>>>> does >>>>>>>>>> not seem to be used in other areas of the standard part of QuantLib. >>>>>>>>>> >>>>>>>>>> Your numerical results below ... >>>>>>>>>> >>>>>>>>>> indicate that the simulated values for ln(S(t)) produced with >>>>>>>>>> timeSteps = 1 >>>>>>>>>> are very likely normally distributed with mean (I use your notation) >>>>>>>>>> ln(S(0)) >>>>>>>>>> + mu*t and standard deviation sigma * sqrt(t) >>>>>>>>>> >>>>>>>>>> This result _*is consistent*_ with: >>>>>>>>>> >>>>>>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw >>>>>>>>>> >>>>>>>>>> _*and*_ >>>>>>>>>> >>>>>>>>>> b) the fact that in QuantLib the PathGenerator.next().value returns the >>>>>>>>>> result >>>>>>>>>> of the SDE expression _*without *_applying the ITO correction >>>>>>>>>> associated >>>>>>>>>> with >>>>>>>>>> the fact that dt is not infinitesimal. >>>>>>>>>> >>>>>>>>>> The b) is also responsible for the lack of convergence in your output >>>>>>>>>> towards >>>>>>>>>> the theoretical target of 68.27% >>>>>>>>>> >>>>>>>>>> You would get the correct convergence if you modified your code by >>>>>>>>>> using >>>>>>>>>> the >>>>>>>>>> expressions below: >>>>>>>>>> >>>>>>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>>>> sqrt(t) * >>>>>>>>>> -1.0); // Sx at -1SD >>>>>>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma * >>>>>>>>>> sqrt(t) * >>>>>>>>>> 1.0); // Sx at +1SD >>>>>>>>>> >>>>>>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t, >>>>>>>>>> sigma) (as >>>>>>>>>> Peter also pointed out) >>>>>>>>>> >>>>>>>>>> Again, due to b) one can produce the correct simulated values of a GBM >>>>>>>>>> diffused quantity (such as a stock price in the GBM model) by using N >>>>>>>>>> time >>>>>>>>>> steps, with N very large. Using N = 1 (like in your example), the >>>>>>>>>> simulated >>>>>>>>>> values will still be lognormally distributed (whence your good result >>>>>>>>>> with >>>>>>>>>> N = >>>>>>>>>> 1), but will be centered at a wrong mean and thus will _*fail *_to >>>>>>>>>> represent >>>>>>>>>> the correct values expected by the GBM SDE. >>>>>>>>>> >>>>>>>>>> Ioannis >>>>>>>>>> >>>>>>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote: >>>>>>>>>>> As said in other posting here, after fixing the test program >>>>>>>>>>> by skipping the first item (the initial price) in the >>>>>>>>>>> generated sample path, BoxMullerGaussianRng now passes the said test. >>>>>>>>>>> >>>>>>>>>>> The bugfixed test code and the new test results can be found here: >>>>>>>>>>> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666 >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> That important fact that the generated sample path contains >>>>>>>>>>> not timeSteps elements but 1 + timeSteps elements needs >>>>>>>>>>> to be documented in the library doc. >>>>>>>>>>> For example on this API doc page one normally would expect >>>>>>>>>>> to find such an important information, but it's missing: >>>>>>>>>>> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> If you or someone else can change/extend the test program by using >>>>>>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would be happy >>>>>>>>>>> to hear about it. Thx. >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28: >>>>>>>>>>>> If you search within the QuantLib code for BoxMullerGaussianRng, you >>>>>>>>>>>> will >>>>>>>>>>>> see >>>>>>>>>>>> it is used only in the experimental folder. It is therefore not >>>>>>>>>>>> surprising if >>>>>>>>>>>> it doesn't produce the expected results. >>>>>>>>>>>> >>>>>>>>>>>> I use myself the MultiPathGenerator with PseudoRandom::rsg_type, >>>>>>>>>>>> which is >>>>>>>>>>>> used >>>>>>>>>>>> extensively in other areas of QuantLib. >>>>>>>>>>>> >>>>>>>>>>>> This type expands to InverseCumulativeRsg< RandomSequenceGenerator< >>>>>>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and gives me >>>>>>>>>>>> good >>>>>>>>>>>> results. >>>>>>>>>>>> >>>>>>>>>>>> Ioannis Rigopoulos, founder of deriscope.com >> >> > > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > Virus-free.www.avast.com > <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> > > > > > > _______________________________________________ > QuantLib-users mailing list > Qua...@li... > https://lists.sourceforge.net/lists/listinfo/quantlib-users > |
|
From: Luigi B. <lui...@gm...> - 2023-09-11 17:01:32
|
You talk as if Ito's lemma was something that one chooses to apply because
one likes it. It's not. Just as there is a rule for performing a change
of variable in an integral or a derivative, there is also a rule for
performing a change of variable in a stochastic differential equation, and
that's Ito's lemma. If you change variables from x to log(x) or whatever
else, you have to use it. If you don't change variables, you don't have to
use it.
This said: in your sample implementation, you generate paths from 0 to t
using a number of timesteps. For the sake of example, let's say t=1 year
and you use 12 timesteps, so one step per month. You calculate the price
plus or minus 1SD, where 1SD is the standard deviation of the expected
distribution at t=1. Then you do the following to calculate percentage
hits:
for (size_t j = 1; j < timeSteps; ++j)
{
const auto Sx = samplePath.at(j);
if (Sx >= m1SD && Sx <= p1SD) ++cHit;
}
that is, for every point in the path you check whether it's between S - 1SD
and S + 1SD. But that makes no sense. You can only check the point at the
end of the path, because it's the only point that belongs to the expected
distribution at t=1 year. In the case of 12 timesteps, the point at j=1
belongs to the expected distribution at t=1 month, which has a different,
much smaller SD. The point at j=2 belongs to the expected distribution at
t=2 months, which has a different SD. Of course your percentages increase
with the timesteps, because you're sampling points from all times between 0
and 1 (which belong to much narrower distributions) and compare them with
the SD at t=1 year. If you compare the points at each timestep with their
correct respective distributions, you'll find about the same percentage
(the usual 68%) falling within the respective SDs.
Hope this helps clear this up—I won't be able to give much more time to
this, if at all.
Regards,
Luigi
On Mon, Sep 11, 2023 at 6:00 PM U.Mutlu <um...@mu...> wrote:
> Ioannis Rigopoulos wrote on 09/11/23 17:28:
> > Thank you U.Mutlu.
> >
> > The histograms indicate that the numbers represent indeed lognormal
> distributions.
> >
> > The average of all numbers in file GBM_ts1_A amounts to *104.5935229*
> >
> > The average of all numbers in file GBM_ts1_B amounts to *104.6047793*
> >
> > Both data sets exhibit the correct log standard deviation, very close to
> *0.3*.
>
> So far so good, and these are the most simplest things,
> but with the rest of your analysis I don't agree.
>
> What do other experts say on the following analysis of Ioannis? :
>
> > *Since you have assumed a drift of 0, the SDE represents a martingale and
> > therefore the average stock price at t=1 owes to be very close to 100.*
> >
> > *This proves that the simulated numbers cannot represent the GBM
> diffusion of
> > a stock price starting at 100 and having drift = 0.*
> >
> > Furthermore, the quantity S(0)exp(0.5*sigma^2*t) equals *104.602786*,
> which is
> > extremely close to the above averages and indicates that the Ito term
> > 0.5*sigma^2*t is embedded in your results, albeit in a wrong fashion,
> since -
> > as mentioned above - it should not affect the average of S(t).
> >
> > It seems that you have set up your code in such a way that ln(S) is
> driftless
> > and therefore S acquires a positive Ito drift of 0.5*sigma^2.
> >
> > It also seems that the number of time steps are much greater than 1,
> otherwise
> > the distribution would appear as normal.
> >
> > At any case, since I cannot run your own code and verify the results, I
> will
> > finish my feedback with this email.
> >
> > Ioannis
>
> Oh oh, that's too much incorrect assumptions & interpretations by you.
> Nevertheless, thank you for attempting to analyse it.
>
> Maybe some other experts too will analyse the data and post their findings
> here.
>
> I claim that Ito's lemma is the culprit and that it's unnecessary
> both in GBM and BSM (Black-Scholes-Merton option pricing formula).
>
>
> > On 9/11/2023 4:19 PM, U.Mutlu wrote:
> >> Ioannis Rigopoulos wrote on 09/11/23 13:43:
> >>> Could you please send us a link to a file that contains the simulated
> values
> >>> produced by your program when you use one time step?
> >>>
> >>> Since you use 1,000,000 simulations, the file should contain 1,000,000
> >>> numbers.
> >>
> >> On the following web page you can find download links to even 2 such
> >> GBM data sets I just have created and uploaded to a server:
> >>
> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861992
> >>
> >>
> >> "
> >> 2 GBM generated test files in CSV format (with header line; see there),
> >> with 1,000,000 generated stock prices (S0=100, IV=30, t=1y, timeSteps=1,
> >> r=0, q=0),
> >> each about 29 MB uncompressed (about 12 MB zip-compressed).
> >>
> >> This GBM generator does not use Ito's lemma.
> >>
> >> ...
> >> "
> >>
> >>
> >>> On 9/11/2023 1:33 PM, U.Mutlu wrote:
> >>>> Ioannis Rigopoulos wrote on 09/11/23 11:47:
> >>>> > By the way, this is how my histogram of S(t) (for t = 1) simulated
> values
> >>>> > looks like when I used 100,000 time steps and 100,000 simulations.
> It is
> >>>> > clearly lognormal in shape.
> >>>>
> >>>> Yes, it's lognormally. My values are lognormlly as well,
> >>>> You just restate the obvious fact.
> >>>>
> >>>> > I am sure, if you do the same when time steps = 1,
> >>>> > you will get a normal distribution for S(t), which would prove
> >>>> > that time steps = 1 does not reproduce the theoretically expected
> >>>> > values for S(t).
> >>>>
> >>>> Not true, simply from the fact that the set boundary is lognormally:
> >>>> -1SD: 74.08182206817179
> >>>> +1SD: 134.9858807576003
> >>>> Ie. it's not symmetric (Gaussian, normally distributed) around the
> >>>> initial mean (100), but is lognormally (and therefore unsymmetric).
> >>>>
> >>>> FYI: in my other test program (that does not use QuantLib)
> >>>> I get in such simulations exact results almost identical
> >>>> to the expected theoretical results! So, your math and arguing is
> wrong, IMO.
> >>>>
> >>>> You are too much fixated on timeSteps=1. Maybe you have
> >>>> a much different understanding of what timeSteps really
> >>>> means in practice. It's _not_ something like the number
> >>>> of simulations that you simply can change and expect
> >>>> to get the same result; no, it has a clearly defined context
> >>>> and gives different results for different timeSteps.
> >>>>
> >>>> And you haven't answered the question why you think it has
> >>>> to be 68.27% for all timeSteps (which IMO is wrong).
> >>>>
> >>>> For those interested, here are more details about the above said:
> >>>>
> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548
> >>>>
> >>>>
> >>>>
> >>>> And, I today did a brief summary of these ongoing discussions on GBM
> here:
> >>>>
> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-5#post-5861919
> >>>>
> >>>>
> >>>>
> >>>>
> >>>> Ioannis Rigopoulos wrote on 09/11/23 11:47:
> >>>>> Our SDE for the evolution of S is:
> >>>>>
> >>>>> dS = μSdt + σSdw
> >>>>>
> >>>>> In your code you use:
> >>>>>
> >>>>>
> >>>>> PathGenerator<RandomSequenceGenerator<MersenneBoxMuller>>
> >>>>> gbmPathGenerator(gbm, t, timeSteps, gsg, false);
> >>>>> const Path& samplePath = gbmPathGenerator.next().value;
> >>>>>
> >>>>> The PathGenerator::next().value is designed to return an array of
> numbers
> >>>>> {S(T₁) , S(T₂) , ... , S(Tᵤ)}
> >>>>>
> >>>>> where
> >>>>>
> >>>>> u := timeSteps
> >>>>>
> >>>>> Δt := t/u
> >>>>>
> >>>>> Tᵢ := iΔt for i = 1,...,u
> >>>>>
> >>>>> It follows:
> >>>>>
> >>>>> Tᵤ = t = the final time when we are interested to know the simulated
> >>>>> values of
> >>>>> S (in your test you use t = 1).
> >>>>>
> >>>>> QuantLib builds the sequence {S(T₁) , S(T₂) , ... , S(Tᵤ)}
> incrementally
> >>>>> through a total number of u steps by adding the quantity ΔS := μSΔt +
> >>>>> σSΔw on
> >>>>> the previous member of the sequence, with Δw := random number drawn
> from
> >>>>> N~(0,Δt¹ᐟ²).
> >>>>>
> >>>>> The previous member of S(T₁) is assumed the T₀ = 0 and the initial
> value
> >>>>> S(T₀)
> >>>>> = S(0) is assumed deterministic and known.
> >>>>>
> >>>>> So, starting with the given S(T₀), QuantLib calculates the S(T₁)
> using the
> >>>>> formula:
> >>>>>
> >>>>> S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw
> >>>>>
> >>>>> Then given the S(T₁) as calculated above, it proceeds to calculate
> S(T₂) as:
> >>>>>
> >>>>> S(T₂) := S(T₁) + μS(T₁)Δt + σS(T₁)Δw
> >>>>>
> >>>>> For simplicity I have used the same notation for the Δw in the
> equations
> >>>>> above, but they are different.
> >>>>>
> >>>>> After u steps, the final value S(t) = S(Tᵤ) is calculated as:
> >>>>>
> >>>>> S(t) = S(Tᵤ) := S(Tᵤ₋₁) + μS(Tᵤ₋₁)Δt + σS(Tᵤ₋₁)Δw
> >>>>>
> >>>>> All equations above hold only in an approximate sense when Δt is
> finite.
> >>>>>
> >>>>> They are accurate only in a limiting sense as Δt → 0.
> >>>>>
> >>>>> It follows that the value S(t) will be calculated correctly only if
> Δt → 0,
> >>>>> which is equivalent to u → ∞
> >>>>>
> >>>>> If you use u = 1 (i.e. one time step), you will still match the
> 68.27% rule
> >>>>> when the 1SD interval is defined in accordance with a normal
> distribution
> >>>>> because S(T₁) := S(T₀) + μS(T₀)Δt + σS(T₀)Δw corresponds to a normal
> >>>>> distribution and if you build a histogram of all simulated values
> you will
> >>>>> notice that S(T₁) is distributed normally, which is clearly not
> consistent
> >>>>> with the SDE that predicts lognormal distribution.
> >>>>>
> >>>>> The correct procedure is of course to use a very large u and define
> the 1SD
> >>>>> interval in accordance with a lognormal distribution, as I described
> in my
> >>>>> previous email.
> >>>>>
> >>>>> Then you will get convergence to the 68.27% rule.
> >>>>>
> >>>>> By the way, this is how my histogram of S(t) (for t = 1) simulated
> values
> >>>>> looks like when I used 100,000 time steps and 100,000 simulations.
> It is
> >>>>> clearly lognormal in shape. I am sure, if you do the same when time
> steps
> >>>>> = 1,
> >>>>> you will get a normal distribution for S(t), which would prove that
> time
> >>>>> steps
> >>>>> = 1 does not reproduce the theoretically expected values for S(t).
> >>>>>
> >>>>> Ioannis
> >>>>>
> >>>>> On 9/11/2023 10:19 AM, U.Mutlu wrote:
> >>>>>> Dear Ioannis Rigopoulos, I now analysed your table and
> >>>>>> must say that you seem to have a much different expectation to,
> >>>>>> and understanding of that table.
> >>>>>>
> >>>>>> The values in your table for timeSteps > 1 can definitely
> >>>>>> not be correct, b/c why do you think it has to be 68.27 for all
> timeSteps?
> >>>>>> This is illogical b/c the more timeSteps the more the hitrate%
> >>>>>> must be, like in my table. It saturates at about 84.93%
> >>>>>>
> >>>>>> Sorry, I have to retract my initial enthusiastic judment
> >>>>>> of your work below, as it's fundamentally flawed & wrong,
> >>>>>> and I should have been more careful when I quickly & briefly
> >>>>>> inspected your posting in the early morning here.
> >>>>>>
> >>>>>> Since your results in the table for timeSteps > 1 are wrong,
> >>>>>> then your formulas must be wrong too, so I had to stop my analysis
> right
> >>>>>> there.
> >>>>>>
> >>>>>> Please just answer why you think it has to be 68.27% for all
> timeSteps?
> >>>>>>
> >>>>>>
> >>>>>> U.Mutlu wrote on 09/11/23 08:36:
> >>>>>>> Mr. Ioannis Rigopoulos, that's a brilliant work by you. Thank you
> very
> >>>>>>> much.
> >>>>>>> I need some time to study and verify the reults, but it looks very
> good
> >>>>>>> to solve the confusion and mystery I had experienced when I
> started this
> >>>>>>> research some weeks ago.
> >>>>>>> I'll try to replicate your results and will report later in detail
> here.
> >>>>>>>
> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29:
> >>>>>>> > On a different topic, I do not understand your statement "/I've
> never
> >>>>>>> seen
> >>>>>>> > that Sx@-1SD and Sx@+1SD get computed with Ito's lemma
> applied/".
> >>>>>>> > What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m
> +1SD,
> >>>>>>> where
> >>>>>>>
> >>>>>>> Sorry for the consfusion, it just means the stock price (Sx after
> time t,
> >>>>>>> ie. S(t) or S_t) for the lower boundary at -1SD and for the upper
> >>>>>>> boundary at
> >>>>>>> +1SD,
> >>>>>>> ie. your "min of S(t)" and "max of S(t)" in your other posting.
> >>>>>>> In future I'll use a more common term like yours, sorry.
> >>>>>>>
> >>>>>>> Thx
> >>>>>>>
> >>>>>>>
> >>>>>>> Ioannis Rigopoulos wrote on 09/10/23 19:29:
> >>>>>>>> Yes, I mean mu * t - 0.5 * sigma * sigma * t
> >>>>>>>>
> >>>>>>>> Are you sure, you also adjusted the boundaries corresponding to
> plus and
> >>>>>>>> minus
> >>>>>>>> 1 SD?
> >>>>>>>>
> >>>>>>>> I ran an Excel spreadsheet by modifying a spreadsheet I had
> published
> >>>>>>>> several
> >>>>>>>> years ago at
> https://blog.deriscope.com/index.php/en/excel-value-at-risk
> >>>>>>>> that
> >>>>>>>> uses the GeometricBrownianMotionProcess , albeit with the
> >>>>>>>> PseudoRandom::rsg_type generator.
> >>>>>>>>
> >>>>>>>> At any case, I got convergence to 68.27% as I increased the
> number of
> >>>>>>>> time
> >>>>>>>> steps and keeping the number of simulations fixed to 100,000,
> provided I
> >>>>>>>> used
> >>>>>>>> the correct (Ito-adjusted) means, as below:
> >>>>>>>>
> >>>>>>>> As you see above, the "No Ito-adjusted" boundaries are 74.08 and
> 134.99,
> >>>>>>>> just
> >>>>>>>> like those you use.
> >>>>>>>>
> >>>>>>>> The Ito-adjusted boundaries are 70.82 and 129.05 produced by
> shifting the
> >>>>>>>> no-adjusted boundaries - in log terms - to the left by
> 0.5*sigma^2*t.
> >>>>>>>>
> >>>>>>>> The simulation results (percentage of values being within the min
> and max
> >>>>>>>> boundaries) are below:
> >>>>>>>>
> >>>>>>>> As you see, the Ito-adjusted percentage reaches 68.30% for
> 100,000 time
> >>>>>>>> steps,
> >>>>>>>> which is close to the expected 68.27%.
> >>>>>>>>
> >>>>>>>> The No Ito-adjusted percentage seems instead to converge to a
> different
> >>>>>>>> value
> >>>>>>>> around 67.73%
> >>>>>>>>
> >>>>>>>> On a different topic, I do not understand your statement "/I've
> never
> >>>>>>>> seen
> >>>>>>>> that Sx@-1SD and Sx@+1SD get computed with Ito's lemma applied/".
> >>>>>>>>
> >>>>>>>> What do the Sx@-1SD and Sx@+1SD mean? I suppose m -1SD and m
> +1SD, where
> >>>>>>>> m is
> >>>>>>>> the mean of a normally distributed variable. This means that one
> needs to
> >>>>>>>> calculate m in order to do the mentioned computation.
> >>>>>>>>
> >>>>>>>> If your starting point is an SDE like dS = S*mu*dt + S*sigma*dw,
> then
> >>>>>>>> ln(S(t))
> >>>>>>>> at any time t will be normally distributed with mean m = ln(S(0))
> +
> >>>>>>>> mu*t -
> >>>>>>>> 0.5*sigma^2*t , where the last term is a consequence of the Ito
> lemma.
> >>>>>>>>
> >>>>>>>> Do you have any reason to believe that the Ito lemma is somehow
> not
> >>>>>>>> valid so
> >>>>>>>> that m = ln(S(0)) + mu*t?
> >>>>>>>>
> >>>>>>>> Ioannis
> >>>>>>>>
> >>>>>>>> On 9/10/2023 3:00 PM, U.Mutlu wrote:
> >>>>>>>>
> >>>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06:
> >>>>>>>>> > new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t,
> sigma)
> >>>>>>>>>
> >>>>>>>>> Thank you for your analysis.
> >>>>>>>>> I guess you mean "mu * t - 0.5 * sigma * sigma * t",
> >>>>>>>>> though in our example with t=1 it doesn't make any difference.
> >>>>>>>>> It gives for timeSteps=1 this result:
> >>>>>>>>>
> >>>>>>>>> $ ./Testing_GBM_of_QuantLib.exe 1
> >>>>>>>>> timeSteps=1 nRuns=1000000 seed=1694350548 S=100.000000
> mu=0.000000
> >>>>>>>>> sigma=0.300000 t=1.000000 : nGenAll=1000000 USE_ITO=1 :
> -1SD=70.822035
> >>>>>>>>> +1SD=129.046162
> >>>>>>>>> cHit=663067/1000000(66.307%)
> >>>>>>>>>
> >>>>>>>>> This is far from the expected value of 68.2689% for timeStep=1.
> >>>>>>>>> So, there must be a bug somewhere. IMO it's our suspect friend
> Ito :-)
> >>>>>>>>>
> >>>>>>>>> Cf. also the following showing similar results for GBM,
> >>>>>>>>> which then was compared to the standard lognormal calculation:
> >>>>>>>>>
> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-3#post-5849548
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> And, honestly, I've never seen that Sx@-1SD and Sx@+1SD get
> computed
> >>>>>>>>> with
> >>>>>>>>> Ito's lemma applied.
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> Ioannis Rigopoulos wrote on 09/10/23 13:06:
> >>>>>>>>>> Thank you for testing the BoxMullerGaussianRng code, which - as
> I
> >>>>>>>>>> wrote -
> >>>>>>>>>> does
> >>>>>>>>>> not seem to be used in other areas of the standard part of
> QuantLib.
> >>>>>>>>>>
> >>>>>>>>>> Your numerical results below ...
> >>>>>>>>>>
> >>>>>>>>>> indicate that the simulated values for ln(S(t)) produced with
> >>>>>>>>>> timeSteps = 1
> >>>>>>>>>> are very likely normally distributed with mean (I use your
> notation)
> >>>>>>>>>> ln(S(0))
> >>>>>>>>>> + mu*t and standard deviation sigma * sqrt(t)
> >>>>>>>>>>
> >>>>>>>>>> This result _*is consistent*_ with:
> >>>>>>>>>>
> >>>>>>>>>> a) the SDE: dS(t) = S(t)*mu*dt + S(t)*sigma*dw
> >>>>>>>>>>
> >>>>>>>>>> _*and*_
> >>>>>>>>>>
> >>>>>>>>>> b) the fact that in QuantLib the PathGenerator.next().value
> returns the
> >>>>>>>>>> result
> >>>>>>>>>> of the SDE expression _*without *_applying the ITO correction
> >>>>>>>>>> associated
> >>>>>>>>>> with
> >>>>>>>>>> the fact that dt is not infinitesimal.
> >>>>>>>>>>
> >>>>>>>>>> The b) is also responsible for the lack of convergence in your
> output
> >>>>>>>>>> towards
> >>>>>>>>>> the theoretical target of 68.27%
> >>>>>>>>>>
> >>>>>>>>>> You would get the correct convergence if you modified your code
> by
> >>>>>>>>>> using
> >>>>>>>>>> the
> >>>>>>>>>> expressions below:
> >>>>>>>>>>
> >>>>>>>>>> const double m1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma
> *
> >>>>>>>>>> sqrt(t) *
> >>>>>>>>>> -1.0); // Sx at -1SD
> >>>>>>>>>> const double p1SD = S * exp(mu * t - 0.5*sigma^2*t + sigma
> *
> >>>>>>>>>> sqrt(t) *
> >>>>>>>>>> 1.0); // Sx at +1SD
> >>>>>>>>>>
> >>>>>>>>>> new GeometricBrownianMotionProcess(S, mu - 0.5*sigma^2*t,
> >>>>>>>>>> sigma) (as
> >>>>>>>>>> Peter also pointed out)
> >>>>>>>>>>
> >>>>>>>>>> Again, due to b) one can produce the correct simulated values
> of a GBM
> >>>>>>>>>> diffused quantity (such as a stock price in the GBM model) by
> using N
> >>>>>>>>>> time
> >>>>>>>>>> steps, with N very large. Using N = 1 (like in your example),
> the
> >>>>>>>>>> simulated
> >>>>>>>>>> values will still be lognormally distributed (whence your good
> result
> >>>>>>>>>> with
> >>>>>>>>>> N =
> >>>>>>>>>> 1), but will be centered at a wrong mean and thus will _*fail
> *_to
> >>>>>>>>>> represent
> >>>>>>>>>> the correct values expected by the GBM SDE.
> >>>>>>>>>>
> >>>>>>>>>> Ioannis
> >>>>>>>>>>
> >>>>>>>>>> On 9/10/2023 11:29 AM, U.Mutlu wrote:
> >>>>>>>>>>> As said in other posting here, after fixing the test program
> >>>>>>>>>>> by skipping the first item (the initial price) in the
> >>>>>>>>>>> generated sample path, BoxMullerGaussianRng now passes the
> said test.
> >>>>>>>>>>>
> >>>>>>>>>>> The bugfixed test code and the new test results can be found
> here:
> >>>>>>>>>>>
> https://www.elitetrader.com/et/threads/simulating-stock-prices-using-gbm.375533/page-4#post-5861666
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>> That important fact that the generated sample path contains
> >>>>>>>>>>> not timeSteps elements but 1 + timeSteps elements needs
> >>>>>>>>>>> to be documented in the library doc.
> >>>>>>>>>>> For example on this API doc page one normally would expect
> >>>>>>>>>>> to find such an important information, but it's missing:
> >>>>>>>>>>>
> https://www.quantlib.org/reference/class_quant_lib_1_1_path_generator.html
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>> If you or someone else can change/extend the test program by
> using
> >>>>>>>>>>> the suggested alternative(s) to BoxMullerGaussianRng, I would
> be happy
> >>>>>>>>>>> to hear about it. Thx.
> >>>>>>>>>>>
> >>>>>>>>>>>
> >>>>>>>>>>> Ioannis Rigopoulos wrote on 09/09/23 16:28:
> >>>>>>>>>>>> If you search within the QuantLib code for
> BoxMullerGaussianRng, you
> >>>>>>>>>>>> will
> >>>>>>>>>>>> see
> >>>>>>>>>>>> it is used only in the experimental folder. It is therefore
> not
> >>>>>>>>>>>> surprising if
> >>>>>>>>>>>> it doesn't produce the expected results.
> >>>>>>>>>>>>
> >>>>>>>>>>>> I use myself the MultiPathGenerator with
> PseudoRandom::rsg_type,
> >>>>>>>>>>>> which is
> >>>>>>>>>>>> used
> >>>>>>>>>>>> extensively in other areas of QuantLib.
> >>>>>>>>>>>>
> >>>>>>>>>>>> This type expands to InverseCumulativeRsg<
> RandomSequenceGenerator<
> >>>>>>>>>>>> MersenneTwisterUniformRng > , InverseCumulativeNormal > and
> gives me
> >>>>>>>>>>>> good
> >>>>>>>>>>>> results.
> >>>>>>>>>>>>
> >>>>>>>>>>>> Ioannis Rigopoulos, founder of deriscope.com
> >>
> >>
> >
> > <
> https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient
> >
> > Virus-free.www.avast.com
> > <
> https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient
> >
> >
> >
> >
> >
> >
> > _______________________________________________
> > QuantLib-users mailing list
> > Qua...@li...
> > https://lists.sourceforge.net/lists/listinfo/quantlib-users
> >
>
>
>
>
> _______________________________________________
> QuantLib-users mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-users
>
|
|
From: U.Mutlu <um...@mu...> - 2023-09-11 17:41:57
|
Luigi, your below posting still has not made to my mail client, maybe the mailman takes longer to send it out to the list members, or again your mail server could not deliver to my mail server. I just saw your posting only in the archive. Luigi B. <lui...@gm...> - 2023-09-11 17:01:32 > for (size_t j = 1; j < timeSteps; ++j) Attention! The above loop is incorrect! It must be for (size_t j = 1; j <= timeSteps; ++j) That's very important! I think many QuantLib users even use this buggy one: for (size_t j = 0; j < timeSteps; ++j) but that one is totally wrong! I'll reply later to the rest of your posting. |
|
From: U.Mutlu <um...@mu...> - 2023-09-11 18:29:17
|
U.Mutlu wrote on 09/11/23 19:41: > Luigi, your below posting still has not made to my mail client, > maybe the mailman takes longer to send it out to the list members, > or again your mail server could not deliver to my mail server. > I just saw your posting only in the archive. > > Luigi B. <lui...@gm...> - 2023-09-11 17:01:32 > > for (size_t j = 1; j < timeSteps; ++j) > > Attention! The above loop is incorrect! > It must be > for (size_t j = 1; j <= timeSteps; ++j) > That's very important! > > I think many QuantLib users even use this buggy one: > for (size_t j = 0; j < timeSteps; ++j) > but that one is totally wrong! > > > I'll reply later to the rest of your posting. I created the table below that shows the expected values for the various timeSteps. It's correctness was also confirmed by simulations. If that table does not answer your questions, let me know and I'll try to clarify. " 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. These numbers are exact values, not the result of simulations. Params used: 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) timeSteps Expected_HitRate 1 68.2689% cf. the above wiki link 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% " Regarding Ito's lemma you write:. > You talk as if Ito's lemma was something that one chooses to apply because > one likes it. It's not. Just as there is a rule for performing a change > of variable in an integral or a derivative, there is also a rule for > performing a change of variable in a stochastic differential equation, and > that's Ito's lemma. If you change variables from x to log(x) or whatever > else, you have to use it. If you don't change variables, you don't have to > use it. As you say, I indeed don't need to change any variables, so then I don't need this Ito's lemma stuff. I'm a practitioner from the field, not a theoretician. I just know this, and can even prove it: Ito's lemma is definitely not needed, neither in GBM nor in Black-Scholes-Merton (BSM). At the moment I can't say anything about the possible use of Ito's lemma in any other field beyond these two, if any. |
|
From: Luigi B. <lui...@gm...> - 2023-09-12 12:19:01
|
I had seen your table, but please, read again the part of my email where I explain that you're calculating it by taking points from different expected distributions. Luigi On Mon, Sep 11, 2023 at 8:29 PM U.Mutlu <um...@mu...> wrote: > U.Mutlu wrote on 09/11/23 19:41: > > Luigi, your below posting still has not made to my mail client, > > maybe the mailman takes longer to send it out to the list members, > > or again your mail server could not deliver to my mail server. > > I just saw your posting only in the archive. > > > > Luigi B. <lui...@gm...> - 2023-09-11 17:01:32 > > > for (size_t j = 1; j < timeSteps; ++j) > > > > Attention! The above loop is incorrect! > > It must be > > for (size_t j = 1; j <= timeSteps; ++j) > > That's very important! > > > > I think many QuantLib users even use this buggy one: > > for (size_t j = 0; j < timeSteps; ++j) > > but that one is totally wrong! > > > > > > I'll reply later to the rest of your posting. > > I created the table below that shows the expected values > for the various timeSteps. > It's correctness was also confirmed by simulations. > If that table does not answer your questions, let me know and I'll try to > clarify. > > " > 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. > These numbers are exact values, not the result of simulations. > Params used: 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) > timeSteps Expected_HitRate > 1 68.2689% cf. the above wiki link > 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% > " > > Regarding Ito's lemma you write:. > > You talk as if Ito's lemma was something that one chooses to apply > because > > one likes it. It's not. Just as there is a rule for performing a > change > > of variable in an integral or a derivative, there is also a rule for > > performing a change of variable in a stochastic differential equation, > and > > that's Ito's lemma. If you change variables from x to log(x) or > whatever > > else, you have to use it. If you don't change variables, you don't > have to > > use it. > > As you say, I indeed don't need to change any variables, > so then I don't need this Ito's lemma stuff. > I'm a practitioner from the field, not a theoretician. > I just know this, and can even prove it: Ito's lemma > is definitely not needed, neither in GBM nor in Black-Scholes-Merton (BSM). > At the moment I can't say anything about the possible > use of Ito's lemma in any other field beyond these two, if any. > > > > |
|
From: U.Mutlu <um...@mu...> - 2023-09-13 00:56:11
|
The said table is intended for this practical situation:
You know (or will know) the stock price only at the end
of the timesteps, and are interested to know from these
datapoints the expectancy% at end for the pre-calculated
range from -1SD to +1SD for t of the initial stock price.
It's the same procedure like for timeSteps=1 (which is easier to understand).
Normally you are interested in the outcome after t (ie. timeSteps=1 giving the
usual 68.27%).
But then you would also like to know the result if you
watch more datapoints than just the one at the end,
ie. in your example timeSteps=12 in t=1y.
Your observation below is correct I think (except
the sentence "But that makes no sense"), but that
all is already counted for in the table, IMO.
For timeSteps=12 in your example the expectancy is 83.5962% (ie. p=0.835962).
Ie. 68.27% now suddenly becomes 83.5982% only b/c we watch the stock price
more times (1 vs 12 times, at end of equal time intervals).
Do you get a different value for the above said situation?
Or asked differently: how should the table ideally be made up instead?
I'll also create a cumulative version for 0 to t_i.
Maybe you mean that one.
This method can be important for bots who watch the
stock price at such fixed intervals, like weeks or months,
and need to know the probability... and/or for (auto-)trading
systems based on such stochastics.
Luigi Ballabio wrote on 09/11/23 19:01:
> in your sample implementation, you generate paths from 0 to t using
> a number of timesteps. For the sake of example, let's say t=1 year and you
> use 12 timesteps, so one step per month. You calculate the price plus or
> minus 1SD, where 1SD is the standard deviation of the expected distribution at
> t=1. Then you do the following to calculate percentage hits:
>
> for (size_t j = 1; j <= timeSteps; ++j)
> {
> const auto Sx = samplePath.at(j);
> if (Sx >= m1SD && Sx <= p1SD) ++cHit;
> }
>
> that is, for every point in the path you check whether it's between S - 1SD
> and S + 1SD. But that makes no sense. You can only check the point at the
> end of the path, because it's the only point that belongs to the expected
> distribution at t=1 year. In the case of 12 timesteps, the point at j=1
> belongs to the expected distribution at t=1 month, which has a different, much
> smaller SD. The point at j=2 belongs to the expected distribution at t=2
> months, which has a different SD. Of course your percentages increase with
> the timesteps, because you're sampling points from all times between 0 and 1
> (which belong to much narrower distributions) and compare them with the SD at
> t=1 year. If you compare the points at each timestep with their correct
> respective distributions, you'll find about the same percentage (the usual
> 68%) falling within the respective SDs.
>
> Hope this helps clear this up—I won't be able to give much more time to this,
> if at all.
|
|
From: Luigi B. <lui...@gm...> - 2023-09-13 09:36:18
|
I just mean that you're taking points at the t=1/12, 2/12 etc and seeing
whether they're inside the SD of the expected distribution at t=1. That
looks to me like comparing apples to oranges. If you have some use for
that, though, good.
Luigi
On Wed, Sep 13, 2023 at 2:55 AM U.Mutlu <um...@mu...> wrote:
> The said table is intended for this practical situation:
> You know (or will know) the stock price only at the end
> of the timesteps, and are interested to know from these
> datapoints the expectancy% at end for the pre-calculated
> range from -1SD to +1SD for t of the initial stock price.
>
> It's the same procedure like for timeSteps=1 (which is easier to
> understand).
>
> Normally you are interested in the outcome after t (ie. timeSteps=1 giving
> the
> usual 68.27%).
> But then you would also like to know the result if you
> watch more datapoints than just the one at the end,
> ie. in your example timeSteps=12 in t=1y.
>
> Your observation below is correct I think (except
> the sentence "But that makes no sense"), but that
> all is already counted for in the table, IMO.
>
> For timeSteps=12 in your example the expectancy is 83.5962% (ie.
> p=0.835962).
> Ie. 68.27% now suddenly becomes 83.5982% only b/c we watch the stock price
> more times (1 vs 12 times, at end of equal time intervals).
>
> Do you get a different value for the above said situation?
> Or asked differently: how should the table ideally be made up instead?
>
> I'll also create a cumulative version for 0 to t_i.
> Maybe you mean that one.
>
> This method can be important for bots who watch the
> stock price at such fixed intervals, like weeks or months,
> and need to know the probability... and/or for (auto-)trading
> systems based on such stochastics.
>
>
> Luigi Ballabio wrote on 09/11/23 19:01:
> > in your sample implementation, you generate paths from 0 to t using
> > a number of timesteps. For the sake of example, let's say t=1 year and
> you
> > use 12 timesteps, so one step per month. You calculate the price plus or
> > minus 1SD, where 1SD is the standard deviation of the expected
> distribution at
> > t=1. Then you do the following to calculate percentage hits:
> >
> > for (size_t j = 1; j <= timeSteps; ++j)
> > {
> > const auto Sx = samplePath.at(j);
> > if (Sx >= m1SD && Sx <= p1SD) ++cHit;
> > }
> >
> > that is, for every point in the path you check whether it's between S -
> 1SD
> > and S + 1SD. But that makes no sense. You can only check the point at
> the
> > end of the path, because it's the only point that belongs to the expected
> > distribution at t=1 year. In the case of 12 timesteps, the point at j=1
> > belongs to the expected distribution at t=1 month, which has a
> different, much
> > smaller SD. The point at j=2 belongs to the expected distribution at t=2
> > months, which has a different SD. Of course your percentages increase
> with
> > the timesteps, because you're sampling points from all times between 0
> and 1
> > (which belong to much narrower distributions) and compare them with the
> SD at
> > t=1 year. If you compare the points at each timestep with their correct
> > respective distributions, you'll find about the same percentage (the
> usual
> > 68%) falling within the respective SDs.
> >
> > Hope this helps clear this up—I won't be able to give much more time to
> this,
> > if at all.
>
>
>
>
>
|
