For the friends who would like to know how the further nested integrals can be calculated here is a small program for third order nesting integral.[$]\int _0^t s \text{$\sigma $dz}(s)\int _0^s \text{$\sigma $dz}(u)\int _0^u \text{$\sigma $dz}(v)[$] function [] = SDE_Compare_Zintegrals_Infiniti_public03( )s1 = RandStream('mt19937ar');RandStream.setDefaultStream(s1);savedState = s1.State; z0=0;sigma0=.710/2;%sigma0=0.710/4%sigma0=0.710*4T_index=1000;%T_index=2000; %T_index=20; dt=.005/4.0;T=(T_index)*dt;paths=1000000; zz(1:paths)=z0; I2(1:paths)=0.0;%integral zdz0 over time by traditional short stepped monte carlo.I1(1:paths)=0.0;zzT(1:paths)=0.0;%One step integral zdz calculated by my method for nn=1:T_index %Loop over time.t=(nn)*dt; t2=(nn-1)*dt; Random1=randn(size(zz)); RandTemp1=Random1.*sqrt(dt);I2=I2+t2*(RandTemp1*sigma0).*I1; %Outer integral I1=I1+(RandTemp1*sigma0).*(zz); %Inner integralzz = zz +RandTemp1*sigma0;end%My One Step method calculation startsRandom1=randn(size(zzT)); zzT=sigma0.^3* Random1.^3.*T.^2.5*sqrt(6)/2.*(1/4/2.5)-...3*sigma0.^3.*T.^2.5*Random1.^1.*(1/4/2.5)*sqrt(6)/2;I2_avg=sum(I2)./pathszzT_avg=sum(zzT)./pathsI2_avg2=sum(abs(I2))./pathszzT_avg2=sum(abs(zzT))./paths%%Below are second moments. To see if both are matchedI2_var=sum(I2.^2)./pathszzT_var=sum(zzT.^2)./pathsBinSize=.0001;%Here you can change the resolution of monte carlo density. You might have to look at the graph and if it is jagged and noisy, increase the bin sizeMaxCutOff=1000;%Max cut off is a value given to density generation program. [XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(I2,paths,BinSize,MaxCutOff );[XDensity2,IndexOutX2,IndexMaxX2] = MakeDensityFromSimulation_Infiniti(zzT,paths,BinSize,MaxCutOff );plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',IndexOutX2(1:IndexMaxX2),XDensity2(1:IndexMaxX2),'r')str=input('look at overlaid Graph comparison of traditional Monte carlo density of zdz(green line) with density of SDE generated from one step monte carlo method(red line)>');end

Last edited by Amin on April 24th, 2016, 10:00 pm, edited 1 time in total.

I will take a few more days to complete the paper but here are a few formulas for friends. Here I give the general algorithm for the type of integrals where only dz and dt integrals follow each other in any order. Like many other possibilities, one four order repeated integral with two dz uncertainties could be[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}\int _0^v\sigma \text{dz}(w)[$]The value of the above kind of integrals, as long as there is no explicit t, does not depend upon the sequence order(whether dz integral comes first or dt integral comes first does not change the value of the total integrals) so variance can be easily calculated by Ito Isometry. The order of hermite polynomial depends upon the number of normal uncertainties so is equal to number of dz integrals. Here is the simple formula to calculate the hermite polynomial representation of the integrals which equals[$] H_n(N)*\frac{\sqrt{\text{Variance}}}{\sqrt{n!}}[$]here n, which is the order of hermite polynomials, equals the number of normal uncertainties. I will give a few examples[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}=H_1(N)\frac{\text{$\sigma $t}^{1.5}}{\sqrt{1!}\sqrt{3}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)=H_2(N)\frac{\sigma ^2t}{\sqrt{2!}\sqrt{2}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_2(N)\frac{\sigma ^2t^2}{\sqrt{2!}\sqrt{12}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=H_1(N)\frac{\sigma ^2t^{2.5}}{\sqrt{1!}\sqrt{20}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_3(N)\frac{\sigma ^3t^{1.5}}{\sqrt{3!}\sqrt{6}}[$]Fourth order and higher integrals can be evaluated to perfect precision using this method so I am not giving examples for that.

Last edited by Amin on April 29th, 2016, 10:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**59685**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: AminI will take a few more days to complete the paper but here are a few formulas for friends. Here I give the general algorithm for the type of integrals where only dz and dt integrals follow each other in any order. Like many other possibilities, one four order repeated integral with two dz uncertainties could be[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}\int _0^v\sigma \text{dz}(w)[$]The value of the above kind of integrals, as long as there is no explicit t, does not depend upon the sequence order(whether dz integral comes first or dt integral comes first does not change the value of the total integrals) so variance can be easily calculated by Ito Isometry. The order of hermite polynomial depends upon the number of normal uncertainties so is equal to number of dz integrals. Here is the simple formula to calculate the hermite polynomial representation of the integrals which equals[$] H_n(N)*\frac{\sqrt{\text{Variance}}}{\sqrt{n!}}[$]here n, which is the order of hermite polynomials, equals the number of normal uncertainties. I will give a few examples[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}=H_1(N)\frac{\text{$\sigma $t}^{1.5}}{\sqrt{1!}\sqrt{3}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)=H_2(N)\frac{\sigma ^2t}{\sqrt{2!}\sqrt{2}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_2(N)\frac{\sigma ^2t^2}{\sqrt{2!}\sqrt{12}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=H_1(N)\frac{\sigma ^2t^{2.5}}{\sqrt{1!}\sqrt{20}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_3(N)\frac{\sigma ^3t^{1.5}}{\sqrt{3!}\sqrt{6}}[$]Fourth order and higher integrals can be evaluated to perfect precision using this method so I am not giving examples for that.what's N?I have having difficulty thinking how these formula can be used in applications.Is it not an idea to write it up in an (SSRN?) scientific paper from beginning to end?

Quotewhat's N?N is standard Normal random draw in case of monte carlo and you could also replace it with normal density for analytics.QuoteI have having difficulty thinking how these formula can be used in applications.I will be explaining that very carefully in the next post.QuoteIs it not an idea to write it up in an (SSRN?) scientific paper from beginning to end?Yes, I will be doing that in a few days. This is just preview.I will also be posting a small monte carlo program in a day that will be second order accurate in T. order(T^2) and could be extended to higher orders using the algorithm I describe here later today.

Last edited by Amin on May 1st, 2016, 10:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**59685**Joined:****Location:**Amsterdam-
**Contact:**

QuoteN is standard Normal random draw ?? cdf, pdf, or something else?Shreve has Hermite H(x,y), but now it is H(N()) where N is cumulative normal? Did I miss something? QuoteI will also be posting a small monte carlo program in a day that will be second order accurate in T. order(T^2) what's T?

Last edited by Cuchulainn on May 1st, 2016, 10:00 pm, edited 1 time in total.

QuoteQuoteN is standard Normal random draw ?? cdf, pdf, or something else?as in simple monte carlo like from box muller or any other normal rng. If you want to use analytics as I mentioned in my post you would have to use pdf.QuoteShreve has Hermite H(x,y), but now it is H(N()) where N is cumulative normal? Did I miss something? You could have used H(x,y) but I think they become equivalent the way I wrote variance out of H(x,y) and used only H(x). I thought it was simpler.QuoteQuoteI will also be posting a small monte carlo program in a day that will be second order accurate in T. order(T^2) what's T?T is the time step or final time.

Cuch, here you could see why those integrals are helpful. I have tried to explain in simple terms.We take a general stochastic differential equation.[$]x(t)=x(0)+\int _0^t\mu (x(s)) \text{ds}+\int _0^t\sigma (x(s)) \text{dz}(s)[$] Eq(1)Where [$]\text{$\mu $(x(s))}[$] and [$]\text{ $\sigma $(x(s))}[$] are x dependent functions and are stochastic due to stochastic nature of x. It is unfortunate that when we use the above notation, we always substitute [$]\text{$\mu $(x(0))}[$] for [$]\text{$\mu $(x(s))}[$] and [$]\text{ $\sigma $(x(0))}[$] for [$]\text{ $\sigma $(x(s))}[$] as is standard in monte carlo and other numerical analysis and little effort has been made to understand the time dependent nature and evolution of [$]\text{$\mu $(x(s))}[$] and [$]\text{ $\sigma $(x(s))}[$] in between the simulation intervals. As I earlier said that \text{$\mu $(x(s))} and [$]\text{ $\sigma $(x(s))}[$] are functions that depend on x and are stochastic due to stochastic nature of x. We define the evolution of these functions as[$]\mu (x(s))=\mu (x(0))+\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dv}+\int _0^s \frac{\partial \mu (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)+.5\int _0^s \frac{\partial ^2\mu (x(v))}{\partial x^2}\sigma ^2(x(v)) \text{dv}[$] Eq(2)[$]\sigma (x(s))=\sigma (x(0))+\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\mu (x(v))\text{dv}+\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)+.5\int _0^s \frac{\partial ^2\sigma (x(v))}{\partial x^2}\sigma ^2(x(v)) \text{dv}[$] Eq(3)We can easily make the observation that x dependent terms in the above integrals from 0 to s could be expanded by application of Ito formula into a time zero term and several other time dependent terms. We substitute equations 2 and 3 in Eq(1) to get[$]x(t)=x(0)+\int _0^t\left(\mu (x(0))+\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dv}+\int _0^s \frac{\partial \mu (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)+.5\int _0^s \frac{\partial ^2\mu (x(v))}{\partial x^2}.\sigma ^2(x(v)) \text{dv}\right) \text{ds}+[$][$]\int _0^t\left(\sigma (x(0))+\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\mu (x(v))\text{dv}+\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)+.5\int _0^s \frac{\partial ^2\sigma (x(v))}{\partial x^2}\sigma ^2(x(v)) \text{dv}\right) \text{dz}(s)[$] Eq(4)Simplifying, we get[$]x(t)=x(0)+\int _0^t\mu (x(0))\text{ds}+\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dvds}+\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)\text{ds}+.5\int _0^t\int _0^s \frac{\partial ^2\mu (x(v))}{\partial x^2}\left.\sigma ^2(x(v)) \text{dv}\right) \text{ds}+\int _0^t\sigma (x(0))\text{dz}(s)+\int _0^t\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\mu (x(v))\text{dvdz}(s)[$][$]+\int _0^t\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)\text{dz}(s)+.5\int _0^t\int _0^s \frac{\partial ^2\sigma (x(v))}{\partial x^2}\sigma ^2(x(v)) \text{dv}\text{dz}(s)[$] Eq(5) In the expression above only two single integrals [$]\int _0^t\mu (x(0))\text{ds}[$]and [$]\int _0^t\sigma (x(0))\text{dz}(s)[$] depend upon time zero value of x and the rest of the double integrals all depend upon forward time dependent values and have to be expanded again similarly and this process could continue until derivatives existed(though we would generally always stop at third or fourth level if we could atain desired precision). We expand just the first double integral in Eq(5) further as an example[$]\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dvds}=\int _0^t\int _0^s \frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\text{dvds}+\int _0^t\int _0^s\int _0^v\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\mu (x(u))\text{dudvds}[$][$]+\int _0^t\int _0^s\int _0^v\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\sigma (x(u))\text{dz}(u)\text{dvds}+\int _0^t\int _0^s\int _0^v.5\frac{\partial ^2}{\partial x^2} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\sigma ^2(x(u))\text{dudvds}[$] Eq(6)We could continue further but if further derivatives continue to exist, we have to ultimately truncate the expansion at some level and we could simply substitute time zero values of drift and volatility for forward time dependent values of drift and volatility. If we decide to stop the expansion at second level, we could write from Eq(5) (Please notice that 2nd equality in equation below is approximate equality)[$]x(t)=x(0)+\int _0^t\mu (x(0))\text{ds}+\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dvds}+\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)\text{ds}+.5\int _0^t\int _0^s \frac{\partial ^2\mu (x(v))}{\partial x^2}\left.\sigma ^2(x(v)) \text{dv}\right) \text{ds}[$][$]+\int _0^t\sigma (x(0))\text{dz}(s)+\int _0^t\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\mu (x(v))\text{dvdz}(s)+\int _0^t\int _0^s \frac{\partial \sigma (x(v))}{\partial x}\sigma (x(v))\text{dz}(v)\text{dz}(s)+.5\int _0^t\int _0^s \frac{\partial ^2\sigma (x(v))}{\partial x^2}\sigma ^2(x(v)) \text{dv}\text{dz}(s)[$][$]=x(0)+\int _0^t\mu (x(0))\text{ds}+\int _0^t\int _0^s \frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\text{dvds}+\int _0^t\int _0^s \frac{\partial \mu (x(0))}{\partial x}\sigma (x(0))\text{dz}(v)\text{ds}+.5\int _0^t\int _0^s \frac{\partial ^2\mu (x(0))}{\partial x^2}\left.\sigma ^2(x(0)) \text{dv}\right) \text{ds}[$][$]+\int _0^t\sigma (x(0))\text{dz}(s)+\int _0^t\int _0^s \frac{\partial \sigma (x(0))}{\partial x}\mu (x(0))\text{dvdz}(s)+\int _0^t\int _0^s \frac{\partial \sigma (x(0))}{\partial x}\sigma (x(0))\text{dz}(v)\text{dz}(s)+.5\int _0^t\int _0^s \frac{\partial ^2\sigma (x(0))}{\partial x^2}\sigma ^2(x(0)) \text{dv}\text{dz}(s)[$][$]=x(0)+\mu (x(0))\int _0^t\text{ds}+\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\int _0^t\int _0^s \text{dvds}+\frac{\partial \mu (x(0))}{\partial x}\sigma (x(0))\int _0^t\int _0^s \text{dz}(v)\text{ds}+.5\frac{\partial ^2\mu (x(0))}{\partial x^2}\sigma ^2(x(0)) \int _0^t\int _0^s \text{dv} \text{ds}[$][$]+\sigma (x(0))\int _0^t\text{dz}(s)+\frac{\partial \sigma (x(0))}{\partial x}\mu (x(0))\int _0^t\int _0^s \text{dvdz}(s)+\frac{\partial \sigma (x(0))}{\partial x}\sigma (x(0))\int _0^t\int _0^s \text{dz}(v)\text{dz}(s)+.5\frac{\partial ^2\sigma (x(0))}{\partial x^2}\sigma ^2(x(0))\int _0^t\int _0^s\text{dv}\text{dz}(s)[$]In case we wanted to go to fourth or hififth or even higher level, we could have continued the repeated expansion and finally replaced forward values at their time zero values but If we only wanted to continue to third level, we could have truncated the further series and replaced forward values by their time zero values and could have included integrals of the kind from Eq(6). Please recall that when we wrote Eq(6), we expanded just one integral and other integral terms have to be similarly expanded.(Please notice that 2nd equality in equation below is approximate equality)[$]\int _0^t\int _0^s \frac{\partial \mu (x(v))}{\partial x}\mu (x(v))\text{dvds}[$][$]=\int _0^t\int _0^s \frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\text{dvds}+\int _0^t\int _0^s\int _0^v\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\mu (x(u))\text{dudvds}+\int _0^t\int _0^s\int _0^v\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\sigma (x(u))\text{dz}(u)\text{dvds}+\int _0^t\int _0^s\int _0^v.5\frac{\partial ^2}{\partial x^2} \left[\frac{\partial \mu (x(u))}{\partial x}\mu (x(u))\right]\sigma ^2(x(u))\text{dudvds}[$][$]=\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\int _0^t\int _0^s \text{dvds}+\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\right]\mu (x(0))\int _0^t\int _0^s\int _0^v\text{dudvds}+\frac{\partial }{\partial x} \left[\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\right]\sigma (x(0))\int _0^t\int _0^s\int _0^v\text{dz}(u)\text{dvds}+.5\frac{\partial ^2}{\partial x^2} \left[\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\right]\sigma ^2(x(0))\int _0^t\int _0^s\int _0^v\text{dudvds}[$]Cuch, here is the answer to your question in an earlier post why these integrals are helpful. As you can see in the above equation, once we have taken initial time zero values of x dependent coefficients, the integrals to be evaluated are of the kind [$]\int _0^t\int _0^s\int _0^v\text{dudvds}[$],[$]\int _0^t\int _0^s\int _0^v\text{dz}(u)\text{dvds}[$], [$]\int _0^t\int _0^s\int _0^v\text{dz}(u)\text{dz}(v)\text{ds}[$] and [$]\int _0^t\int _0^s\int _0^v\text{dz}(u)\text{dz}(v)\text{dz}(s)[$] which I mentioned how to evaluate in a previous post using ito isometry and simple hermite polynomials.

Last edited by Amin on May 1st, 2016, 10:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**59685**Joined:****Location:**Amsterdam-
**Contact:**

Amin,I take the classic SDE/FDM approach and then use it in the MC simulator. I do not have enough knowledge of the integral approach to do MC with your formulae.I give the code in an attempt to possibly allow the Hermite stuff to be used in the path evolver stuff. Do I plug in your drift and diffusion functions into my code?You see, I'm a white belt here:)Thanks// Edited code. The loop is key!class SDE{ // Defines drift + diffusion + dataprivate: std::shared_ptr<OptionData> data; // The data for the optionpublic: SDE(const OptionData& optionData) : data(new OptionData(optionData)) {} double drift(double t, double S) { // Drift term return (data->r - data->D)*S; // r - D } double diffusion(double t, double S) { // Diffusion term return data->sig * S; }};int main(){ std::cout << "1 factor MC with explicit Euler\n"; // OptionData myOption { 100.0, 1.0, 0.06, 0.2, 0.03, 1 }; // Uniform initialisation OptionData myOption(( OptionParams::strike = 100.0, OptionParams::expiration = 1.0, OptionParams::volatility = 0.2, OptionParams::dividend = 0.03, OptionParams::optionType = 1, OptionParams::interestRate = 0.06)); SDE sde(myOption); // Initial value of SDE double S_0 = 100.0; // Variables for underlying stock double x; double VOld = S_0; double VNew; long NT = 100; std::cout << "Number of time steps: "; std::cin >> NT; // V2 mediator stuff long NSIM = 50000; std::cout << "Number of simulations: "; std::cin >> NSIM; double M = static_cast<double>(NSIM); double k = myOption.T / double (NT); double sqrk = std::sqrt(k); // Normal random number double dW; double price = 0.0; // Option price double payoffT; double avgPayoffT = 0.0; double squaredPayoff = 0.0; double sumPriceT = 0.0; // Normal (0,1) rng std::default_random_engine dre; std::normal_distribution<double> nor(0.0, 1.0); long coun = 0; // Number of times S hits origin for (long i = 1; i <= M; ++i) { // Calculate a path at each iteration if ((i/10000) * 10000 == i) {// Give status after each 1000th iteration std::cout << i << std::endl; } VOld = S_0; x = 0.0; for (long index = 0; index <= NT; ++index) { // Create a random number dW = nor(dre); // The FDM (in this case explicit Euler), equation (9.2) from the text VNew = VOld + (k * sde.drift(x, VOld)) +(sqrk * sde.diffusion(x, VOld) * dW); VOld = VNew; // Spurious values if (VNew <= 0.0) coun++; x += k; } // Assemble quantities (postprocessing) payoffT = myOption.myPayOffFunction(VNew); sumPriceT += payoffT; avgPayoffT += payoffT/M; avgPayoffT *= avgPayoffT; squaredPayoff += (payoffT*payoffT); } // Finally, discounting the average price price = std::exp(-myOption.r * myOption.T) * sumPriceT/M; std::cout << "Price, after discounting: " << price << ", " << std::endl; std::cout << "Number of times origin is hit: " << coun << std::endl; return 0;}

Last edited by Cuchulainn on May 2nd, 2016, 10:00 pm, edited 1 time in total.

Daniel, I will do simple math behind the sde you are simulating. Here is the equation that we expanded just one level more than traditional monte carlo. This is equation right after equation six in previous post. I will write it again.[$]x(t)=x(0)+\mu (x(0))\int _0^t\text{ds}+\frac{\partial \mu (x(0))}{\partial x}\mu (x(0))\int _0^t\int _0^s \text{dvds}+\frac{\partial \mu (x(0))}{\partial x}\sigma (x(0))\int _0^t\int _0^s \text{dz}(v)\text{ds}+.5\frac{\partial ^2\mu (x(0))}{\partial x^2}\sigma ^2(x(0)) \int _0^t\int _0^s \text{dv} \text{ds}[$][$]+\sigma (x(0))\int _0^t\text{dz}(s)+\frac{\partial \sigma (x(0))}{\partial x}\mu (x(0))\int _0^t\int _0^s \text{dvdz}(s)+\frac{\partial \sigma (x(0))}{\partial x}\sigma (x(0))\int _0^t\int _0^s \text{dz}(v)\text{dz}(s)+.5\frac{\partial ^2\sigma (x(0))}{\partial x^2}\sigma ^2(x(0))\int _0^t\int _0^s\text{dv}\text{dz}(s)[$] Eq(1)I have tried to use your notation.we use [$]\text{$\mu $(x(0))=(r-d)x(0)}[$][$]\frac{\partial \mu (x(0))}{\partial x}\text{=(r-d)}[$][$]\frac{\partial ^2\mu (x(0))}{\partial x^2}\text{=0}[$][$]\sigma\text{(x(0))= sig x(0)}[$][$]\frac{\partial \sigma (x(0))}{\partial x}\text{=sig}[$][$]\frac{\partial ^2\sigma (x(0))}{\partial x^2}\text{=0}[$]substituting the above values in Eq(1), we get[$]x(t)=x(0)+(r-d)x(0) \text{$\Delta $t}+(r-d)^2x(0)\frac{\text{$\Delta $t}^2}{2}+(r-d)\text{sig} x(0)\int _0^t\int _0^s \text{dz}(v)\text{ds}[$][$]+\text{sig} x(0)\text{nor}(0,1)\sqrt{\text{$\Delta $t}}+\text{sig} (r-d)x(0)\int _0^t\int _0^s \text{dvdz}(s)+\text{sig}^2 x(0)\int _0^t\int _0^s \text{dz}(v)\text{dz}(s)[$]We use the identities [$]\int _0^t\int _0^s \text{dz}(v)\text{ds}=\int _0^t\int _0^s \text{dvdz}(s)=\text{nor}(0,1)\frac{\text{$\Delta $t}^{1.5}}{\sqrt{3}}[$]and [$]\int _0^t\int _0^s \text{dz}(v)\text{dz}(s)=\left(\text{nor}(0,1)^2-1\right)\frac{\text{$\Delta $t}}{2}[$]Now we can get the monte carlo simulation scheme.[$]x(t+\text{$\Delta $t})=x(t)+(r-d)*x(t) \text{$\Delta $t}+(r-d)^2x(t)\frac{\text{$\Delta $t}^2}{2}+\text{sig}* x(t)\text{nor}(0,1)\sqrt{\text{$\Delta $t}}+2\text{sig}* (r-d)x(t)\text{nor}(0,1)\frac{\text{$\Delta $t}^{1.5}}{\sqrt{3}}[$][$]+\text{sig}^2x(t)\left(\text{nor}(0,1)^2-1\right)\frac{\text{$\Delta $t}}{2}[$]The above is order [$]\text{$\Delta $t}[$] accurate with some terms of higher order. Adding one more term that comes from further expansion will take the order to [$]\text{$\Delta $t}[$]^1.5. The new better scheme would be[$]x(t+\text{$\Delta $t})=x(t)+(r-d)*x(t) \text{$\Delta $t}+(r-d)(r-d)x(t)\frac{\text{$\Delta $t}^2}{2}+\text{sig}* x(t)\text{nor}(0,1)\sqrt{\text{$\Delta $t}}+2*\text{sig}* (r-d)x(t)\text{nor}(0,1)\frac{\text{$\Delta $t}^{1.5}}{\sqrt{3}}[$][$]+\text{sig}^2*x(t)\left(\text{nor}(0,1)^2-1\right)\frac{\text{$\Delta $t} }{2}+\text{sig}^3*x(t)\left(\text{nor}(0,1)^3-3\text{nor}(0,1)\right)\frac{\text{$\Delta $t}^{1.5}}{6}[$]

Last edited by Amin on May 2nd, 2016, 10:00 pm, edited 1 time in total.

- Cuchulainn
**Posts:**59685**Joined:****Location:**Amsterdam-
**Contact:**

Thanks. Ahsan, that's pretty clear.Let me get my head around to get it into the MC framework and I'll hopefully get back with some results and queries in a few days.

Last edited by Cuchulainn on May 3rd, 2016, 10:00 pm, edited 1 time in total.

Thank you Daniel. We all Wilmotters await the results of your monte carlo work. I am looking forward to your analysis, comments and questions.Though for mildly non-linear SDEs, we can take a quite large step, for many non-linear SDE equations, a universally much longer step becomes difficult especially close to zero where derivatives change extremely fast. I am presenting an improvement using Girsanov theorem so that a large enough step could be taken by even non-linear equations. I will upload that program here on Wilmott in a day hopefully that could do monte carlo with Girsanov theorem and we will be able to take large steps for non-linear SDEs.I had earlier said on the forum that we would be able to take large steps of even 5-10 years in monte carlo or other density construction. I stand behind my claim but that has to be done by using the concepts from my earlier work Stochastic Calculus of Standard Deviations. In stochastic calculus of standard deviations, we divide the initial delta density into 100-400 normal standard deviation subdivisions that preserve the mass in each subdivision and continue to advance ahead in time in the same normal standard deviation line. We can effortlessly take step size of roughly .25 years and continue to update the SDE derivatives of each subdivision in the direction of that particular normal SD. Unlike recombining tree, we move ahead in a non-recombining tree where number of tree branches remain exactly the same through delta origin (start of SDE simulation) to the end of the simulation and complexity remains exactly the same as we advance the density into time. It will be impossible for many non-linear SDEs to take the derivatives only at time zero and advance the solution till several years since it may become too difficult and unefficient to keep on taking derivatives of drift and volatility functions. However I am sure it can still be done extremely cheaply using the concepts from stochastic calculus of standard deviations.

- Cuchulainn
**Posts:**59685**Joined:****Location:**Amsterdam-
**Contact:**

Just to make sure; your SDE is?dS = (r-d)Sdt + sig S dW?It has an exact solution.

Last edited by Cuchulainn on May 4th, 2016, 10:00 pm, edited 1 time in total.

Yes, of course, we all know the lognormal density solution after transformation to log coordinates but my scheme has no knowledge of the exact solution. It is based on the naiive simulation technique used in your code and a comparison has to be made how better it is as compared to scheme you used in code.

Sorry about a careless mistake in Eq(6) in the fourth post from the bottom on previous page. I took a second derivative in the last term while I should have used Ito product rule. When I wrote the post in response to cuch's post, I was going to sleep in the evening after staying awake for more than 24 hours but decided to answer his post first before I sleep so I made the mistake. I think most forum members have a better mathematics than mine, so they would have understood the mistake on their own but if you did not notice it, please substitute Ito product rule instead of simple change of one variable Ito formula I carelessly used there.

GZIP: On