I will try to give wilmott friends an idea about my new research on stochastic differential equations and their simulation. My formal research paper will be ready in a week but I will present some preview details about the research for the friends for discussion.Let us write a simple stochastic differential equation in the general form with a finite step here as[$]x(t)=\int _0^t\mu (x(t)) \text{dt}+\int _0^t\sigma (x(t)) \text{dz}(t)\text{ }\text{Eq (1)}[$]When we simulate the above stochastic differential equation, the simplest method is to freeze the drift and volatility coefficients at their initial value, and multiply with dt and Normal(sqrt(t)) as[$]x(t)=\mu (x(0)) \text{dt}+\sigma (x(0)) \text{sqrt}(t)z.[$]Though many better techniques exist to simulate the stochastic differential equations for a finite step or to find their density, most of them lack generalization and have to be applied in special settings. We derive a generic new method that can simulate even non-linear stochastic differential equations to arbitrary accuracy. The method is extremely simple and is based on the simple understanding of how the non-linear coefficient functions are changing in time. I have tried to give some intuition about it in the exposition below.Before, we embark on a full-fledged scheme for the above general SDE, we would like to give the reader some crucial insights. To give the proper insights, we would like to start by taking a normal diffusion so that we could proceed in a step by step manner [$]\text{dx}(s)=\sigma \text{dz}(s)[$] Eq(2)or[$]x(t)=x(0)+\int _0^t\sigma \text{dz}(s)[$] Eq(3)We would like to start by learning how to calculate the following time integral with extreme precision and arbitrary accuracy[$]I(t)=\int _0^t\mu (x(s)) \text{ds}[$] Eq(4)We would also like to explain how to calculate the following stochastic integral with extreme precision and arbitrary accuracy[$]S(t)=\int _0^t\sigma (x(s)) \text{dz}(s)[$] Eq(5)Since we want to proceed step by step, we will first calculate these integrals for a normal brownian motion as given in Eq (2) and after gaining intuition , we later generalize the SDE from Eq(2) to Eq(1) and calculate I(t) and S(t) with the most general dynamics as given by Eq(1).We start by calculating I(t) as given by Eq(4) and Eq(2).[$]\int _0^t\mu (x(s))\text{ds}=t \mu (x(t)) -\int _0^ts d[\mu (x(s))][$] Eq(6)Expanding the second term on RHS[$]\int _0^ts d[\mu (x(s))]=\int _0^t \frac{\partial \mu (x(s))}{\partial x}s\text{$\sigma $dz}(s)+.5\int _0^t \frac{\partial ^2\mu (x(s))}{\partial x^2}s\sigma ^2 \text{ds}[$] Eq(7)We know that [$]\frac{\partial \mu (x(s))}{\partial x}[$] in first term of Eq(7) is a function of x(s) and we can apply Ito Change of variable formula on this as [$]\frac{\partial \mu (x(s))}{\partial x}=\frac{\partial \mu (x(0))}{\partial x}+\int _0^s \frac{\partial ^2\mu (x(u))}{\partial x^2}\text{$\sigma $dz}(u)+.5\int _0^s \frac{\partial ^3\mu (x(u))}{\partial x^3}\sigma ^2 \text{du}[$] Eq(8)Now we come to the most important step towards understanding of this expansion. We expand the first noise term on RHS of Eq(7) by substituting Eq(8) in it. [$]\int _0^t \frac{\partial \mu (x(s))}{\partial x}s\text{$\sigma $dz}(s)=\int _0^t \frac{\partial \mu (x(0))}{\partial x}s\text{$\sigma $dz}(s)+\int _0^t s \text{$\sigma $dz}(s)\int _0^s \frac{\partial ^2\mu (x(u))}{\partial x^2}\text{$\sigma $dz}(u) +\int _0^t s\text{$\sigma $dz}(s)\int _0^s .5 \frac{\partial ^3\mu (x(u))}{\partial x^3}\sigma ^2 \text{du}\text{ }[$] Eq(9)So we have a few nested integrals and we can easily deal with them. We could continue the expansions of [$]\frac{\partial ^2\mu (x(u))}{\partial x^2}[$] and [$]\frac{\partial ^3\mu (x(u))}{\partial x^3}[$] to further order just by following Eq(8) but we stop the order here for initial simplicity and explanation of the basic example. We freeze [$]\frac{\partial ^2\mu (x(u))}{\partial x^2}[$] and [$]\frac{\partial ^3\mu (x(u))}{\partial x^3}[$] in the above equation to their initial time zero values (we could have continued to expand them to desired accuracy if we wished) and write Eq(9) as[$]\int _0^t \frac{\partial \mu (x(s))}{\partial x}s\text{$\sigma $dz}(s)=\frac{\partial \mu (x(0))}{\partial x}\int _0^t s\text{$\sigma $dz}(s)+\frac{\partial ^2\mu (x(0))}{\partial x^2}\int _0^t s \text{$\sigma $dz}(s)\int _0^s \text{$\sigma $dz}(u) +.5 \frac{\partial ^3\mu (x(0))}{\partial x^3}\int _0^t s\text{$\sigma $dz}(s)\int _0^s \sigma ^2 \text{du}\text{ }[$] Eq(10)Now we are left to evaluation of intgrals[$]\int _0^t \text{s$\sigma $dz}(s)[$] Eq(11)and double nested integrals[$]\int _0^t s \text{$\sigma $dz}(s)\int _0^s \text{$\sigma $dz}(u)=\int _0^t s \sigma ^2\text{ }(z(s)-z(0))\text{dz}(s)[$][$]\int _0^t s\text{$\sigma $dz}(s)\int _0^s \sigma ^2 \text{du}=\int _0^t s\sigma \left(\left.\sigma ^2 s\right)\text{dz}(s)\right.[$]None of the above nested integrals vanishes and their continuations do not go to zero. I have devised schemes in my research that can calculate all different types of above integrals with arbitrary complexity of nesting to perfect precision using very simple integral calculus in a form that is a function of only time t, and a unit gaussian.We would similarly expand the second term of Eq(7) on RHS using Ito expansion of [$]\frac{\partial ^2\mu (x(s))}{\partial x^2}[$] and substituting it in Eq(7) . We will freeze the coefficients of above expansion at time zero like we did to get Eq(10) and get different integrals like Eq(11-13). All these stochastic integrals can be trivially solved by scheme I have devised in my research.To be continued tomorrow. Sorry, tired of writing more at the moment.

It might be you lost term x ( 0 ) on the right hand side (1). The simplest discrete scheme is Euler's one. You can find a stepwise approximation of the solution of the nonlinear SDE. You can present [$]\Delta t[$] order convergence of the approximate to exact solution. Your scheme looks to be assigned to only first step on [ 0 , [$]\Delta t[$] ] interval. In other words you should make first a stepwise recurrent continuation of the finite scheme on arbitrary interval [0 , T ]. Next it looks reasonable to present estimate [$](\Delta t )^n[$] of the convergence of the approximate solution to exact solution.

QuoteOriginally posted by: list1It might be you lost term x ( 0 ) on the right hand side (1). Sorry it is a mistake.Quote Your scheme looks to be assigned to only first step on [ 0 , [$]\Delta t[$] ] interval. In other words you should make first a stepwise recurrent continuation of the finite scheme on arbitrary interval [0 , T ].No it is totally general. Once you understand my logic, you can easily work out the mathematics.

Quoteperfect precision There is no such thing as perfect precision in (numerical) mathematics. Accuracy is a polynomial function of dt based on continuity of the unknown exact solution. Which is what @list has also said, when I look at it. Is convergence monotonic?

QuoteOriginally posted by: AminQuoteOriginally posted by: list1It might be you lost term x ( 0 ) on the right hand side (1). Sorry it is a mistake.Quote Your scheme looks to be assigned to only first step on [ 0 , [$]\Delta t[$] ] interval. In other words you should make first a stepwise recurrent continuation of the finite scheme on arbitrary interval [0 , T ].No it is totally general. Once you understand my logic, you can easily work out the mathematics.In all equation s you used 0 on the low bounds of the integral besides you used x ( 0 ) in eq (10). It suggests that you consider only first step on finite time partition [$] { t\,_k } [$] of the arbitrary interval [ 0 , T ]. You should consider a partition of the [ 0 , T ] and then to apply your method on each subinterval [$] [\, t\,_k , t\,_{k + 1}\, ) [$]

QuoteOriginally posted by: CuchulainnQuoteperfect precision There is no such thing as perfect precision in (numerical) mathematics. Accuracy is a polynomial function of dt based on continuity of the unknown exact solution. Which is what @list has also said, when I look at it. Is convergence monotonic?why do you put "numerical" in brackets, makes me nervous

QuoteOriginally posted by: pcaspersQuoteOriginally posted by: CuchulainnQuoteperfect precision There is no such thing as perfect precision in (numerical) mathematics. Accuracy is a polynomial function of dt based on continuity of the unknown exact solution. Which is what @list has also said, when I look at it. Is convergence monotonic?why do you put "numerical" in brackets, makes me nervousGood point. I was thinking that someone gives a closed/exact solution which they might think is perfectly precise until you want numbers:)

In the previous post, I mentioned a general technique of expanding stochastic integrals in a way that we can utilize higher derivatives(than the two derivatives of stochastic integrals as is the current best practice) in evaluation of stochastic integrals and in simulation . We can continue to take higher derivatives, and freeze the derivatives in the tail at their zero value and later evaluate the stochastic integrals composed merely of z(t) and t and expand them in powers of standard normal N(0,1) and a function of time and volatility.I will mention a caveat that in following exposition normal SDE is [$]z(t)=z(0)+\int _0^t\sigma \text{dz}(s)[$]I will explain with some basic examples how nested integrals mentioned in the previous post can be trivially simulated with perfect accuracy to any length of time. The detailed evaluation methodology will follow in a later post.[$]\int _0^t \text{s$\sigma $dz}(s)[$] Eq(11)can be simulated as [$]\int _0^t \text{s$\sigma $dz}(s)=\frac{2.175}{\sqrt{2\text{Pi}}}\sigma \frac{t^{1.5}}{1.5} N(0,1)[$]where N(0,1) indicates a single standard normal draw in simulation.[$]\int _0^t z(s)\text{ds}=z(0)t+\frac{2.175}{\sqrt{2\text{Pi}}}\sigma \frac{t^{1.5}}{1.5} N(0,1)[$]We give another general simulation expansion as[$]\int _0^t z(s)\text{$\sigma $dz}(s)=\frac{2.5}{\sqrt{2\text{Pi}}}z(0)\sigma \sqrt{t}N(0,1) +.5\frac{2.5}{\sqrt{2\text{Pi}}}\sigma ^2t\text{ }N(0,1)^2-.5\sigma ^2t[$]The above simulation for the stochastic integral is valid for any length of time as is the case with simulation of all other nested integrals that I will mention. For example, the above simulation equation can be simply verified by most members of the forum by setting up a simulation in 5-10 minutes. Similarly taking Equation 12 in the first post, we have a stepping technique that remains perfectly accurate for all time horizons.C1, C2, C3 in the below are fixed constants.[$]\int _0^t s \text{$\sigma $dz}(s)\int _0^s \text{$\sigma $dz}(u)=\int _0^t s \sigma \text{ }(z(s)-z(0))\text{dz}(s)=.5 \text{C2} \sigma ^2\frac{t ^2}{2} N(0,1)^2-.5\sigma ^2\frac{t^2}{2}[$]while the general version of above integral can be solved as[$]\int _0^t s \sigma \text{ }z(s)\text{dz}(s)=\text{C1} z(0)\sigma \frac{t^{1.5}}{1.5}N(0,1)+.5 \text{C2} \sigma ^2\frac{t ^2}{2} N(0,1)^2-.5\sigma ^2\frac{t^2}{2}[$]I give the last integral for today here as[$]\int _0^t s\text{$\sigma $dz}(s)\int _0^s \sigma ^2 \text{du}=\int _0^t s\sigma \left(\left.\sigma ^2 s\right)\text{dz}(s)\right.[$] Eq(13)which can be simulated as[$]\int _0^t s\text{$\sigma $dz}(s)\int _0^s \sigma ^2 \text{du}=\int _0^t s\sigma \left(\left.\sigma ^2 s\right)\text{dz}(s)=.5 \text{C3} \sigma ^3\frac{t^{2.5}}{2.5}N(0,1)\right.[$]Most friends would have realized that in a complex simulation we can put together the coefficients of powers of standard Normal and simply multiply those coefficients with the N(0,1) powers and add to get the result.

Amin , your formulas might make sense only for [ 0 , t ] where t is a small number. There are a large number of approaches have been developed in the theory of SDEs for the construction of the DISCRETE time approximations of the SDEs solutions. It might make sense to look through them and might be contact to some of the authors introducing yours approach.

list1, please read what I have stated. These formulas are valid for all time. Please do not give any comments out of your usual intuition. I will accept your comments if you have spent 5 minutes in setting up a small simulation to test my formulas and then say something meaningful about the results you got. Now do not write a word here without setting up the simulation. You will be most welcome to comment after you have actually done that.

QuoteOriginally posted by: Aminlist1, please read what I have stated. These formulas are valid for all time. Please do not give any comments out of your usual intuition. I will accept your comments if you have spent 5 minutes in setting up a small simulation to test my formulas and then say something meaningful about the results you got. Now do not write a word here without setting up the simulation. You will be most welcome to comment after you have actually done that.Wrong answer IMO. Imagine if you said that at a training session. That's the way @Quantigic reacts.list has a point. It is up to the author to prove the benefits of a new method, not the reader! All we see is a number of equations which (at least for me) have to be taken on blind faith.Take a concrete case would be a good idea. And where can these formulae be used?If the author spends 5 minutes => save the 90K forum members each spending 5 minutes. Quotewe can put together the coefficients of powers of standard Normal and simply multiply those coefficients with the N(0,1) powers and add to get the result.I wonder how this pans out. I smell a series..

It might be I am wrong but let check your approximation in a trivial case [$]\mu\,=\, 0\, , \, \sigma\,=\,1[$] and t = 10.

Sorry about being a bit curt. Thank you cuch for your suggestion. I am quickly pasting the code for one example integral zdz. Please copy it in your matlab editor or simply clone in any other language. As I have set up the simulation, you can trivially modify the simple monte carlo using equations I gave in previous post vs how we traditionally do an accurate short stepping monte carlo. Excuse any possible minor errors as I am not being formal. @Cuch and others, your feedback is requested.[$]\int _0^t z(s)\text{$\sigma $dz}(s)=\frac{2.5}{\sqrt{2\text{Pi}}}z(0)\sigma \sqrt{t}N(0,1) +.5\frac{2.5}{\sqrt{2\text{Pi}}}\sigma ^2t\text{ }N(0,1)^2-.5\sigma ^2t[$]Here is the main function.function [] = SDE_Compare_Zintegrals_Infiniti_public( )s1 = RandStream('mt19937ar'); RandStream.setDefaultStream(s1); savedState = s1.State; z0=1;%z0=10;%z0=.5;sigma0=0.710%sigma0=0.710/4%sigma0=0.710*4T_index=200; %T_index=2000; %T_index=20; dt=.005/1.0;T=(T_index)*dt;paths=1000000; zz(1:paths)=z0; I2(1:paths)=0.0;%integral zdz over time by traditional short stepped monte carlo.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+(RandTemp1*sigma0).*(zz); zz = zz +RandTemp1*sigma0; end %My One Step method calculation starts Random1=randn(size(zzT)); zzT=z0*sigma0.*Random1*T.^.5*2.5/sqrt(2*pi)+... %Excuse the coefficients 2.5/sqrt(2pi). It is unity. I followed as I wrote on internet. Please modify. .5*sigma0.^2* Random1.^2.*T*2.5/sqrt(2*pi)-... .5*sigma0.^2.*T; I2_avg=sum(I2)./paths zzT_avg=sum(zzT)./paths I2_avg2=sum(abs(I2))./paths zzT_avg2=sum(abs(zzT))./pathsBinSize=.05;%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)>'); endHere is the helper function to graph from monte carlofunction [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti(X,Paths,BinSize,MaxCutOff )%Processes monte carlo paths to return a series Xdensity as a function of IndexOut. IndexMax is the maximum value of index.% Xmin=0; Xmax=0; for p=1:Paths if(X(p)>MaxCutOff) X(p)=MaxCutOff; end if(Xmin>real(X(p))) Xmin=real(X(p)); end if(Xmax<real(X(p))) Xmax=real(X(p)); end end IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1 XDensity(1:IndexMax)=0.0; for p=1:Paths index=floor(real(X(p)-Xmin)/BinSize+.5)+1; if((index)<1) index=1; end if(real(index)>IndexMax) index=IndexMax; end XDensity(index)=XDensity(index)+1.0/Paths/BinSize; end IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize; end

Amin,I would be interested in testing your formula in a 1-factor MC option pricer. if possible, could you delineate the steps? Thx.

Yes, Cuch, I will post it here for a lognormal option pricer or CEV option pricer. But please give me two or three days since I am busy working out a few last details of the method and its derivation. I hope you do not mind a few days. In the meantime, you could play around with the integrals I mentioned in the previous post using the matlab code I pasted.