
amin


Total Posts: 279 
Joined: Aug 2005 


I will try to give NuclearPhynance 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
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
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 nonlinear stochastic differential equations to arbitrary accuracy. The method is extremely simple and is based on the simple understanding of how the nonlinear coefficient functions are changing in time. I have tried to give some intuition about it in the exposition below.
Before, we embark on a fullfledged 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
Eq(2) or Eq(3)
We would like to start by learning how to calculate the following time integral with extreme precision and arbitrary accuracy Eq(4)
We would also like to explain how to calculate the following stochastic integral with extreme precision and arbitrary accuracy
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).
Eq(6)
Expanding the second term on RHS
Eq(7)
We know that in first term of Eq(7) is a function of x(s) and we can apply Ito Change of variable formula on this as
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.
Eq(9)
So we have a few nested integrals and we can easily deal with them. We could continue the expansions of and 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 and 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
Eq(10)
Now we are left to evaluation of intgrals
Eq(11)
and double nested integrals
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 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(1113). 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.





amin


Total Posts: 279 
Joined: Aug 2005 


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
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.
can be simulated as
where N(0,1) indicates a single standard normal draw in simulation.
We give another general simulation expansion as
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 510 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.
while the general version of above integral can be solved as
I give the last integral for today here as Eq(13) which can be simulated as
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


Total Posts: 279 
Joined: Aug 2005 


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. Your feedback is requested.
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*4
T_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=(nn1)*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))./paths
BinSize=.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 size MaxCutOff=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
%%%%%%%%%%%%%%Above function has ended<<<<<< Here is the helper function to graph from monte carlo
function [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 Xmax=real(X(p)); end end IndexMax=floor((XmaxXmin)/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:(IndexMax1))*BinSize; end 




amin


Total Posts: 279 
Joined: Aug 2005 


For the interest of NP friends I will post slightly modified version of the previous program that computes the integral of tzdz or t(zz0)dz by putting z0=0. If my one step method has a different shape, it is more accurate than the shape of density from very short stepping method. You can look at mean and second moment numbers that they match very well. You will need the helper matlab function from previous post to plot the density. I will be posting a full fledged simulation of CEV model option pricing program in 24 days.
function [] = SDE_Compare_Zintegrals_Infiniti_public02( ) s1 = RandStream('mt19937ar'); RandStream.setDefaultStream(s1); savedState = s1.State; z0=4; %z0=0; %z0=10; %z0=.5;
sigma0=.710/2; %sigma0=0.710/4 %sigma0=0.710*4
T_index=400;% 00; %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 tzdz over time by traditional short stepped monte carlo. zzT(1:paths)=0.0;%One step integral tzdz calculated by my method for nn=1:T_index %Loop over time. t=(nn)*dt; t2=(nn1)*dt;
Random1=randn(size(zz)); RandTemp1=Random1.*sqrt(dt);
I2=I2+t2*(RandTemp1*sigma0).*(zz); zz = zz +RandTemp1*sigma0; end
%My One Step method calculation starts Random1=randn(size(zzT));
zzT=z0*sigma0.*Random1*T.^1.5*.5760+... sigma0.^2* Random1.^2.*T.^2*sqrt(2).*(1/4)... .5*sigma0.^2.*T.^2*(1/2)*sqrt(2);
I2_avg=sum(I2)./paths zzT_avg=sum(zzT)./paths I2_avg2=sum(abs(I2))./paths zzT_avg2=sum(abs(zzT))./paths
%%Below are second moments. To see if both are matched I2_var=sum(I2.^2)./paths zzT_var=sum(zzT.^2)./paths
BinSize=.01;%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 size MaxCutOff=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 



amin


Total Posts: 279 
Joined: Aug 2005 


Though I just gave a detailed scheme for a double nested dz integral, many intelligent and smart mathematicians would have guessed that a doubly nested dz integral can be described in terms of hermite polynomials of degree up to 2 and similarly the integrals with nesting degree 3 could be described in terms of hermite polynomials of degree three and less and so on for higher degree of nesting. Of course, you have to make adjustments due to volatility and time multipliers that affect the total variance. 




amin


Total Posts: 279 
Joined: Aug 2005 


For a reference to above formulas with nested dzintegrals when you have calculated the appropriate variance, you can calculate them by a slight modification after following the reference below.
"Brownian Motion and Stochastic Calculus" by Karatzas and Shreve. You have to follow the Exercise 3.31 on page 167. 



amin


Total Posts: 279 
Joined: Aug 2005 


Most friends would realize that we can do a much high order monte carlo than existing practices. Most of the integrals involve simple and basic arithmetic with a minor application of integral calculus. If you take high enough order in expansion of drift and volatility part of SDE, you could do a 5yr monte carlo with desired accuracy in just one step. In the context of my first post, this technique has never been applied to monte carlo before.
My upcoming paper also has a portion devoted to very interesting applications of Girsanov theorem to nonlinear SDEs with drift.
I will also try to write a C++ code for NP friends because that is what I will do myself finally.
A word of caution, Karatzas and Shreve reference I mentioned will have to be slightly modified but you could see that different nested integrals in our case remain orthogonal to each other which is probably the only requirement for a general application.





amin


Total Posts: 279 
Joined: Aug 2005 


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.
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*4
T_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=(nn1)*dt;
Random1=randn(size(zz)); RandTemp1=Random1.*sqrt(dt); I2=I2+t2*(RandTemp1*sigma0).*I1; %Outer integral I1=I1+(RandTemp1*sigma0).*(zz); %Inner integral zz = zz +RandTemp1*sigma0; end
%My One Step method calculation starts Random1=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)./paths zzT_avg=sum(zzT)./paths I2_avg2=sum(abs(I2))./paths zzT_avg2=sum(abs(zzT))./paths
%%Below are second moments. To see if both are matched I2_var=sum(I2.^2)./paths zzT_var=sum(zzT.^2)./paths
BinSize=.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 size MaxCutOff=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 



amin


Total Posts: 279 
Joined: Aug 2005 


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
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 unit 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
here n, which is the order of hermite polynomials, equals the number of normal uncertainties. I will give a few examples
Fourth order and higher integrals can be evaluated to perfect precision using this method so I am not giving examples for that. 




amin


Total Posts: 279 
Joined: Aug 2005 


We take a general stochastic differential equation.
Eq(1)
Where and 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 for and for 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 and in between the simulation intervals. As I earlier said that and are functions that depend on x and are stochastic due to stochastic nature of x. We define the evolution of these functions as Eq(2)
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 Eq(4)
Simplifying, we get
Eq(5)
In the expression above only two single integrals and 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
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)
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)
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 , , and which I mentioned how to evaluate in a previous post using ito isometry and simple hermite polynomials 



amin


Total Posts: 279 
Joined: Aug 2005 


Guys, very sorry. But you have to follow the missing discussion on W****tt. I feel truly bad mentioning W****tt here but due to so many different thing happening, I could not keep up with posting here regularly. My apologies.
http://W****tt.com/messageview.cfm?catid=4&threadid=99702
Here I restart the discussion. The method can be applied to ODEs as well.
Though I think most intelligent friends would already have realized how my expansion of variable coefficients into iterated integrals with constant integrands in the case of SDEs as I have shown in the previous discussion applies to ODEs, I will still like to give a general example.
we have an ODE of the form Eq(1)
where a(X) and b(X) are totally general variable coefficients. In order to make exposition more clear, we emphasize the implicit dependence of X on t. X depends on t but not necessarily explicitly and a(X) and b(X) are coefficients that depend explicitly on X but would usually have dependence on t only through implicit dependence of X on t.
Eq(2)
We start at time t=t0 and want to find the value of X at a t=t1. We have for variable coefficient a(X) Eq(3)
Let us see how we can improve on this expansion
Eq(4)
substituting eq(4) in Eq(3), we get
Similarly, we can expand again to find.
Since and are constants, we can write
We have to realize that the higher order terms(>3 in the above equation) continue and their addition will continue to add more accuracy.
Many mathematicians would have realized that dt integral terms when considered alone without constant coefficients converge to exponential of dt. But we have now constant weights on each term of the Taylor expansion of the exponential and this expansion with constant weights gets equated with the solution of ODEs with variable coefficients.
Please note that this method can be used for other kind of ODES like 2nd order with variable coefficients and many other ODEs and so on. You could probably use it for many PDEs as well. You have to be ingenious.
Is it(for ODEs) not exactly the same methodology and same series of steps what we have been doing for SDEs previously.
dX/dt = f(X,t) Eq(A)
Intelligent friends would know that the above most general ODE in Eq(A) can very easily be solved by the above method of iterated integrals that I proposed. The general ODE solution expansion just would not have constant integrands as we see for the simplest type for which I did the solution exercise. In the solution of most general ODE as in Eq(A), the iterated integrals solution would possibly keep some terms in t that would have to be integrated in t.





amin


Total Posts: 279 
Joined: Aug 2005 


Please check this attached program for simulation of CEV noise using SCSD technology. Attached File: SDEXticsInfinitiTechnologiesCEVNoise02_Public.zip In case you cannot download it. Here is the code.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [] = SDEXticsInfinitiTechnologiesCEVNoise02_Public( ) X0_max=800;%Total number of SD Fractions. X0_mid=400.5;%Mid of the SD Grid. Lies at z0 by definition. X0_mid0=400; X0_start=1; dX0=.0125;%*8.0;%Width of SD cells in terms of Standard Deviations.
T_index=100*1; dt=.01/1.0;
epsilon=0.70 gamma=.650;
z0=1; dz_integral2(1:X0_max)=z0; iivector(1:X0_max)=1:X0_max; paths=1000000; %No. of Paths for creation of Monte Carlo density by simulating z, F(z); zz4(1:paths)=z0;
const=.398942*exp(dX0^2.*(.25+(iivectorX0_mid).^2)).*(exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2)+exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2));;
const2=.398942*(dX0.*exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*(.5(iivectorX0_mid))+... dX0.*exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*(.5+(iivectorX0_mid))... 1.25331*erf(dX0.*(.353553+.707107*(iivectorX0_mid)))+... 1.25331*erf(dX0.*(.353553+.707107*(iivectorX0_mid))));
const3=.398942.^1*(exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*(2+(.5(iivectorX0_mid)).^2.*dX0^2)+... exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*(2(.5+(iivectorX0_mid)).^2*dX0^2));
const4=.398942.^1*(exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*dX0.*(1.53*(iivectorX0_mid)+(.5(iivectorX0_mid)).^3.*dX0^2)+... exp(.5*dX0^2.*(.5+(iivectorX0_mid)).^2).*dX0.*(1.5+3*(iivectorX0_mid)+(.5+(iivectorX0_mid)).^3.*dX0^2)... 3.75994*erf(dX0.*(.353553+.707107*(iivectorX0_mid)))+... 3.75994*erf(dX0.*(.353553+.707107*(iivectorX0_mid))));
prob2=exp(.5*((iivectorX0_mid).*dX0).^2)./sqrt(22/7.0*2.0)*dX0; prob2=normcdf(((iivectorX0_mid).*dX0)+.5*dX0)normcdf(((iivectorX0_mid).*dX0).5*dX0);
% plot((1:X0_max),const(1:X0_max),'g',(1:X0_max),const2(1:X0_max),'r',(1:X0_max),const3(1:X0_max),'b',(1:X0_max),const4(1:X0_max),'y') % str=input('enter a key'); % % plot((1:X0_max),const(1:X0_max)./prob2(1:X0_max),'g',(1:X0_max),const2(1:X0_max)./prob2(1:X0_max),'r',(1:X0_max),const3(1:X0_max)./prob2(1:X0_max),'b',(1:X0_max),const4(1:X0_max)./prob2(1:X0_max),'y') % str=input('enter a key');
for nn=1:T_index %Loop over time. t=(nn)*dt; t0=(nn1)*dt; prob=exp(.5*((iivectorX0_mid).*dX0).^2)./sqrt(22/7.0*2.0)/(sqrt(t)*epsilon); %Probability in each SD Grid cell.
dt=tt0;
dz_integral2=dz_integral2+... epsilon.*dz_integral2.^gamma.*(sqrt(t)sqrt(t0)).*(const./prob2)... +1.0*epsilon.*gamma.*dz_integral2.^(gamma1).*(epsilon.*dz_integral2.^(gamma)).*(tt0)/2.0.*(((const2)1))... +epsilon.^3.*((gamma*(gamma1).*dz_integral2.^(gamma2+2*gamma))+gamma.^2.*dz_integral2.^(2*gamma2+gamma)).*(t.^1.5t0.^1.5)./6.*((const3)const.*3)... +0*epsilon.^4.*(gamma*(gamma1)*(gamma2).*dz_integral2.^(gamma3+3*gamma)+4*gamma.^2*(gamma1).*dz_integral2.^(gamma2+gamma1+2*gamma)+... gamma.^3*dz_integral2.^(3*gamma3+gamma)).*(t.^2t0.^2)./24.*(const46*const2+3)+... +.5*(epsilon.^3.*(gamma*(gamma1).*dz_integral2.^(gamma2+2*gamma))).*(t.^1.5/1.5t0.^.5*t+.5*t0.^1.5/1.5).*(const);
% CONST0=epsilon.*Z_power1.*(sqrt(t)sqrt(t0)).*(const./prob2); % CONST1=epsilon.*dz_integral2.^gamma.*(sqrt(t)sqrt(t0)).*(const./prob2); % CONST2=epsilon.*gamma.*dz_integral2.^(gamma1).*(epsilon.*dz_integral2.^(gamma)).*(tt0)/2.0.*(((const2)1)); % CONST3=epsilon.^3.*((gamma*(gamma1).*dz_integral2.^(gamma2+2*gamma))+gamma.^2.*dz_integral2.^(2*gamma2+gamma)).*(t.^1.5t0.^1.5)./6.*((const3)const.*3); % CONST4=epsilon.^4.*(gamma*(gamma1)*(gamma2).*dz_integral2.^(gamma3+3*gamma)+4*gamma.^2*(gamma1).*dz_integral2.^(gamma2+gamma1+2*gamma)+... % gamma.^3*dz_integral2.^(3*gamma3+gamma)).*(t.^2t0.^2)./24.*(const46*const2+3); % % CONST4=.5*(epsilon.^3.*(gamma*(gamma1).*dz_integral2.^(gamma2+2*gamma))).*(t.^.5t0.^.5).*sqrt(dt).*(const./prob2); % %plot(dz_integral2(1:X0_max),CONST0(1:X0_max),'k',dz_integral2(1:X0_max),CONST1(1:X0_max),'g',dz_integral2(1:X0_max),CONST2(1:X0_max),'r',dz_integral2(1:X0_max),CONST3(1:X0_max),'b',dz_integral2(1:X0_max),CONST4(1:X0_max),'y') %str=input('enter a key');
%Naive brute force Monte carlo part follows. Since the steps are extremely %short the monte carlo is uite close to the right density for most of the cases %but may be off for a very small number of cases but it is a good initial %check on our analytic density
Random1=randn(size(zz4)); RandTemp1=Random1.*sqrt(dt); zz4 =zz4+zz4.^gamma.*RandTemp1*epsilon; %dz_integral2(dz_integral2<0)=0.0; zz4(zz4<0)=0; end effective_sigma=epsilon; Jacobian_dz1(1:X0_max)=0; Jacobian_dz2(1:X0_max)=0; %dz_integral2(dz_integral2<0)=0.0; dz_integral2(dz_integral2<0)=0.0; for ii=2:X0_max1
%Jacobian_dz2(ii)=((dz_integral2(ii2))8*(dz_integral2(ii1))+8*(dz_integral2(ii+1))(dz_integral2(ii+2)))/(12*dX0*sqrt((T_index1+1)*dt)*effective_sigma); Jacobian_dz2(ii)=((dz_integral2(ii+1))(dz_integral2(ii1)))/(2*dX0*sqrt((T_index1+1)*dt)*effective_sigma);
if(Jacobian_dz2(ii)~=0.0) Jacobian_dz2(ii)=1/Jacobian_dz2(ii); end end
dz_integral2_avg=sum(dz_integral2.*prob2)
zz4_avg=sum(zz4)./paths
BinSize=.0025;%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 size %If it is made of straight line increments, decrease the bin size. MaxCutOff=20; [XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(zz4,paths,BinSize,MaxCutOff ); str=input('Press any key to see the graph, Please make sure to rescale the graph from Edit menu on the matlab graph when diffusion reaches zero.') plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',(dz_integral2(X0_start+2:X0_max2)),prob(X0_start+2:X0_max2).*Jacobian_dz2(X0_start+2:X0_max2),'r') str=input('look at overlaid Graph comparison of Monte carlo density of SDE(green line) with density of SDE generated from analytic method(red line)>');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



amin


Total Posts: 279 
Joined: Aug 2005 


Here are some quick notes about the previous program. The diffusion I am analytically simulating is dX=epsilon X^gamma dz(t). The results are very encouraging. However the accuracy sharply deteriorates when significant mass of the diffusion gets absorbed in zero. If diffusion reaches zero but significant mass does not get absorbed in zero, the results are still very good. The reason is that we have divided the total density into standard deviation fractions. The most negative standard deviations expand the fastest in negative direction and the most positive standard deviation expand fastest in the positive direction at a rate proportional to standard deviation count and so on for intermediate values. When significant mass is absorbed at zero, a large fraction of the negative standard deviations is absorbed in zero but remaining more positive standard deviation fractions continue to advance. By contrast, when we simulate by monte carlo, the mass that is away from zero contains all standard deviations and hence a method that mimicks Gaussian must contain all standard deviations. Our analytic method works by advancing all standard deviations independently. When significant standard deviation fractions are absorbed at zero, we must reassign standard deviation fractions to all the probability mass that has not reached zero so that it contains all standard deviation fractions. This is however not very simple or straightforward to do and would require an algorithmic effort to do it. Another way to approach the problem would be to reflect the negative standard deviation fractions that reach zero back into positive domain with negative amplitude that cancels the increased amplitude due to lack of negative standard deviation. This will also require some effort and I have not worked on it. Close vicinity of zero is not a problem for the algorithm when significant mass is not absorbed. The only problem is that negative standard deviation fractions that get absorbed in zero are lacked in the positive body/domain of the noise diffusion where all standard deviation fractions must be present for a perfectly Gaussian algorithm. Some way around it can be worked out with a bit of effort. After a few years (>14)the accuracy starts to decrease. Here I will like to mention that my estimate of mean is only based on a rough summation so you will have to write your own integration routine on subdivisions each with different width. Here I would like to remind readers of my density correction technique that changes the density so as to make its mean equal to a given number while keeping the mass equal to unity. I presented it in an old program in which I calculated transition densities with Girsanov. I will find that out again and explain it here in a day or two. That routine can be very handy for removing small biases in the density in a way that density moves in the desired direction and its mass remains unity and average over the density becomes equal to the desired average. With the help of that small density correction routine, you can freely use the above method even where large accuracy is required for the CEV noises where significant mass is not absorbed at zero for short to medium time horizons. I will explain that routine here probably in my next post. 




amin


Total Posts: 279 
Joined: Aug 2005 

 

amin


Total Posts: 279 
Joined: Aug 2005 


In light of my post on this thread dated 20160502 18:31 and another post dated: 20160430 15:23, I want to mention that I had fully specified the method of iterated integrals solution to SDEs as early as May 2016. Though I have written a formal paper on iterated integrals solution to ODEs, I have not written a formal paper on iterated integrals solution to SDEs. If any friend wants to do further research and discover new things using method of iterated integrals beyond what I have done or include method of iterated integrals in a book, I would request them to please give reference to my research that I deserve. I would be very thankful to friends who would give me credit that I deserve due to my original research.

Below is my post on this thread Posted: 20160502 18:31
We take a general stochastic differential equation.
Eq(1)
Where and 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 for and for 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 and in between the simulation intervals. As I earlier said that and are functions that depend on x and are stochastic due to stochastic nature of x. We define the evolution of these functions as Eq(2)
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 Eq(4)
Simplifying, we get
Eq(5)
In the expression above only two single integrals and 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
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)
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)
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 , , and which I mentioned how to evaluate in a previous post using ito isometry and simple hermite polynomials

Here is another post on this forum Posted: 20160430 15:23
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
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 unit 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
here n, which is the order of hermite polynomials, equals the number of normal uncertainties. I will give a few examples
Fourth order and higher integrals can be evaluated to perfect precision using this method so I am not giving examples for that. 




amin


Total Posts: 279 
Joined: Aug 2005 


Here I uploaded my paper https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3133959
The title is: Method of Iterated Integrals for Solution of Stochastic Integrals with Applications to Monte Carlo Simulation of Stochastic Differential Equations.
Abstract: In this article we suggest a new method for solutions of stochastic integrals where the dynamics of the variables in integrand are given by some stochastic differential equation. We also propose numerical simulation of stochastic differential equations which is based on iterated integrals method employing Ito change of variable expansion of the drift and volatility function of the SDE. We repeatedly apply the Ito change of variable formula or Ito product rule to integrand and find its higher order derivatives and later we find the contribution from these higher order derivatives to the integrand using higher order iterated integrals. First order derivatives of the integrand, for example, would employ second order iterated integrals to find their contribution to the solution. We evaluate these iterated integrals terms at the initial time of the simulation along with the solution of the stochastic integrals which is found in terms of Hermite polynomials and variance of the integrals. We apply the method of iterated integrals to simulation of various stochastic differential equations. We find that the accuracy of simulation sharply increases using the method of iterated integral as compared to naïve and simple monte carlo simulation methods. 



amin


Total Posts: 279 
Joined: Aug 2005 


Sorry, I made minor changes to paper. The legends for some graphs were bad and one graph was misplaced and I had to change that.
I want to ask friends what is the right reference for section 3 of my paper. When I did it, I followed a problem in "Brownian motion and stochastic calculus" in exercises where the solution was given. I took cue from it, played slightly around and came up with the formula in section 3. But I am sure the precise problem must have been solved decades ago. Can somebody point me to appropriate reference and I will add it to my paper. 




amin


Total Posts: 279 
Joined: Aug 2005 


Ok, this is an older post for a related reference.
For a reference to above formulas with nested dzintegrals when you have calculated the appropriate variance, you can calculate them by a slight modification after following the reference below.
"Brownian Motion and Stochastic Calculus" by Karatzas and Shreve. You have to follow the Exercise 3.31 on page 167.
Can somebody point to appropriate reference.




amin


Total Posts: 279 
Joined: Aug 2005 


I made some changes in the paper. Some wording in the abstract and the paper was ambiguous and could be wrongly interpreted. Here is the link to paper. https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3133959
New Abstract:
In this article we suggest a new method for solutions of stochastic integrals where the dynamics of the variables in integrand are given by some stochastic differential equation. We also propose numerical simulation of stochastic differential equations which is based on iterated integrals method employing Ito change of variable expansion of the drift and volatility function of the SDE. Method of iterated integrals is a method of solution of integrals in which we repeatedly expand different terms in the integral as an initial value and an integrated derivative so that their sum would equal the original term in the integrand that was expanded. We continue to repeatedly expand the integrals until we have enough accuracy in the terms evaluated at initial values so as to be able to neglect the highest derivative term. We evaluate all these iterated integrals terms at the initial time of the simulation along with the solution of the stochastic integrals which is found in terms of Hermite polynomials and variance of the integrals. We apply the method of iterated integrals to simulation of various stochastic differential equations. We find that the accuracy of simulation sharply increases using the method of iterated integral as compared to naïve and simple monte carlo simulation methods. 




amin


Total Posts: 279 
Joined: Aug 2005 


I am giving mathematica with which you can now solve any SDE or their stochastic integrals to arbitrary order. I will follow soon with examples and graphs. Restriction is that you cannot have an explicit t or its powers in the SDE and it can solve the most of the SDEs of the type that usually appear in math finance usually with only dt. If you want a program that also calculates the solution with explicit t in drift or volatility, you have to contact me or change this program. I have posted explanation of code on W****tt technical forum.
Clear[x, x0, Y, p, l1, l2, Z, mu, sigma, mu0, sigma0, k, m, ll1, ll2, YAns, H]; p = 4; ll1 = 0; ll2 = 0; Array[Y, {p + ll1 + 1, p + ll2 + 1}, {0, 0}]; For[k = 0, k <= p + ll1, k++, (For[m = 0, m <= p + ll2, m++, (Y[k, m] = 0);]);]; Y[ll1, ll2] = x; mu := mu0 *x; sigma := sigma0*x; For[k = 0, k Y[l1 + 1, l2] = Y[l1 + 1, l2] + D[Y[l1, l2], x]*mu + .5*D[D[Y[l1, l2], x], x]*sigma^2; Y[l1, l2 + 1] = Y[l1, l2 + 1] + D[Y[l1, l2], x]*sigma;)];)]; Array[H, p + ll2 + 1, 0]; YAns = 0; For[k = 0, k <= p, k++, (For[m = 0, m <= k, m++, (l1 = ll1 + k  m; l2 = ll2 + m; YAns = YAns + Y[l1, l2]*Sqrt[(t^(2*l1 + l2)*(2*l1)!)/((2*l1 + l2)! (l1!)^2)]*H[l2]/Sqrt[l2!];)];)]; YAns
The output again is given as x H[0] + Sqrt[t^2] (0. + mu0 x) H[0] + 1/2 Sqrt[t^4] (0. + mu0^2 x) H[0] + 1/6 Sqrt[t^6] (0. + mu0^3 x) H[0] + 1/24 Sqrt[t^8] (0. + mu0^4 x) H[0] + sigma0 Sqrt[t] x H[1] + ( Sqrt[t^3] (0. + 2 mu0 sigma0 x) H[1])/Sqrt[3] + ( Sqrt[t^5] (0. + 3 mu0^2 sigma0 x) H[1])/(2 Sqrt[5]) + ( Sqrt[t^7] (0. + 4 mu0^3 sigma0 x) H[1])/(6 Sqrt[7]) + 1/2 sigma0^2 Sqrt[t^2] x H[2] + ( Sqrt[t^4] (0. + 3 mu0 sigma0^2 x) H[2])/(2 Sqrt[6]) + ( Sqrt[t^6] (0. + 6 mu0^2 sigma0^2 x) H[2])/(4 Sqrt[15]) + 1/6 sigma0^3 Sqrt[t^3] x H[3] + ( Sqrt[t^5] (0. + 4 mu0 sigma0^3 x) H[3])/(6 Sqrt[10]) + 1/24 sigma0^4 Sqrt[t^4] x H[4]




amin


Total Posts: 279 
Joined: Aug 2005 


I am giving mathematica with which you can now solve any SDE or their stochastic integrals to arbitrary order. I will follow soon with examples and graphs. Restriction is that you cannot have an explicit t or its powers in the SDE and it can solve the most of the SDEs of the type that usually appear in math finance usually with only dt. If you want a program that also calculates the solution with explicit t in drift or volatility, you have to contact me or change this program. I have posted explanation of code on W****tt technical forum.
Clear[x, x0, Y, p, l1, l2, Z, mu, sigma, mu0, sigma0, k, m, ll1, ll2, YAns, H]; p = 4; ll1 = 0; ll2 = 0; Array[Y, {p + ll1 + 1, p + ll2 + 1}, {0, 0}]; For[k = 0, k <= p + ll1, k++, (For[m = 0, m <= p + ll2, m++, (Y[k, m] = 0);]);]; Y[ll1, ll2] = x; mu := mu0 *x; sigma := sigma0*x; For[k = 0, k Y[l1 + 1, l2] = Y[l1 + 1, l2] + D[Y[l1, l2], x]*mu + .5*D[D[Y[l1, l2], x], x]*sigma^2; Y[l1, l2 + 1] = Y[l1, l2 + 1] + D[Y[l1, l2], x]*sigma;)];)]; Array[H, p + ll2 + 1, 0]; YAns = 0; For[k = 0, k <= p, k++, (For[m = 0, m <= k, m++, (l1 = ll1 + k  m; l2 = ll2 + m; YAns = YAns + Y[l1, l2]*Sqrt[(t^(2*l1 + l2)*(2*l1)!)/((2*l1 + l2)! (l1!)^2)]*H[l2]/Sqrt[l2!];)];)]; YAns
The output again is given as x H[0] + Sqrt[t^2] (0. + mu0 x) H[0] + 1/2 Sqrt[t^4] (0. + mu0^2 x) H[0] + 1/6 Sqrt[t^6] (0. + mu0^3 x) H[0] + 1/24 Sqrt[t^8] (0. + mu0^4 x) H[0] + sigma0 Sqrt[t] x H[1] + ( Sqrt[t^3] (0. + 2 mu0 sigma0 x) H[1])/Sqrt[3] + ( Sqrt[t^5] (0. + 3 mu0^2 sigma0 x) H[1])/(2 Sqrt[5]) + ( Sqrt[t^7] (0. + 4 mu0^3 sigma0 x) H[1])/(6 Sqrt[7]) + 1/2 sigma0^2 Sqrt[t^2] x H[2] + ( Sqrt[t^4] (0. + 3 mu0 sigma0^2 x) H[2])/(2 Sqrt[6]) + ( Sqrt[t^6] (0. + 6 mu0^2 sigma0^2 x) H[2])/(4 Sqrt[15]) + 1/6 sigma0^3 Sqrt[t^3] x H[3] + ( Sqrt[t^5] (0. + 4 mu0 sigma0^3 x) H[3])/(6 Sqrt[10]) + 1/24 sigma0^4 Sqrt[t^4] x H[4]





amin


Total Posts: 279 
Joined: Aug 2005 


Sorry, I had earlier promised to post a matlab program for calculation of solution of distribution of SDEs and their stochastic integrals. I wrote the program a few days ago with monte carlo and it works well on many SDEs but has problems when I apply it to some other SDEs. One reason I think is that instead of integrating over normal distribution, I started by using randomly drawn normal numbers. And with larger times and higher order of hermrite polynomials, I continued to see problems due to bad higher order distributional properties of pseudo random numbers. If somebody wishes, I would be willing to share my program with monte carlo. I would be writing and testing another matlab program to integrate the hermite polynomials with normal density and hopefully share it after the weekend. Sometimes I think that such a method should probably be called hermite transform method to solve the SDEs because I think that they are in some sense analogous to what sine and cosine with different frequencies are to fourier transforms and I am sure some of the similar problems might arise when we integrate higher order hermite polynomials just like the sort of problems arise when high frequencies are integrated in fourier transforms. We can of course know that hermite polynomials are orthogonal basis functions just like sines and cosines with different frequencies are basis functions for fourier transforms. An obvious difference is that high frequencies in fourier are mostly seen in short term options. But comparing with fourier transforms, we see that as fourier with zero frequency gives us the mean, here hermite polynomial of order zero gives us the mean. It seems that many times it could possibly be more important to get lower order hermite polynomials right as compared to adding larger order hermite polynomials. May be we can take series in hermite polynomials upto possibly eighth order and take larger number of terms so that coefficient terms with each of these low order hermite polynomials converge. Taking the extreme case, we can see if series of terms with hermite polynomial of order zero does not coverge we will never have a right distribution since mean will be off even if we have included the 32nd order hermite polynomials in the expansion. But this would have to be very carefully seen and properly researched. But we can probably do a rough analysis of convergence once we know how and when the time coefficients that arise from the integration and calculation of variance converge and when factorial takes over the higher order powers of terminal time and based on that and combined with already calculated coefficients of hermite polynomials a simple analysis of convergence could be possible. Though lack of convergence in coefficient terms of higher order hermite polynomials would be a problem but lack of convergence in the series of terms which are coefficients of lower order hermite polynomials would be a far larger problem. With the current expansion, I do not believe that it would be possible to solve all the SDEs universally with something like hermite polynomials to 16th order which could possibly work only till a year in some of my simulations. Really it could be a possibility to use higher order hermite polynomials only if we use a carefully chosen numerical discretization scheme instead of doing monte carlo simulations. Ok, I have really not thought about it carefully but can we somehow calculate delta and gamma using this method as a possible byproduct of the solution technique? I will try to think about it later but any thoughts from friends? I also plan to work on some simple algorithm to add term structure of volatilities to existing method. 



amin


Total Posts: 279 
Joined: Aug 2005 


I have worked with a new ItoTaylor tree method for calculation of densities of SDEs and their path integrals. In this method, we start an ItoTaylor tree for the SDE from origin. Each branch of ItoTaylor tree is like an ODE for a particular value of the standard normal variable. And we recompute the ItoTaylor coefficients given the new initial values of X for every branch along the ItoTaylor tree for every time grid. We have to account for variances.Basically ItoTaylor tree method is just converting an SDE evolution into ODElike evolution where each ODE evolves with the particular value of normal assigned to the specific branch of ItoTaylor tree.
The method is lightning fast since the tree branches once at origin and never branches after that. I have taken only 81 branches and there are only 64 time steps per year and one step computational effort for every branch is roughly comparable to monte carlo simulation computational effort for one single path. So you can yourself imagine how superfast the method is. Its speed is roughly comparable to monte carlo with only 81 paths.
I will soon be writing a formal research paper about this method.
We analytically compute the hermite polynomials weights for a particular SDE in terms of initial value parameter i.e starting value of the SDE variable, X0.
I have lumped together all time calculations related to variances and drifts into calculation of hermite polynomial weights and later, in the evolution of Itotree, I solely use variances of standard normal random variable without any further consideration of variance or drift dependent time values.
Here is an ItoTaylor tree matlab program based on only first two hermite polynomials. This is just a first initial version but I decided to post it. You can change the program in all sort of interesting ways. And you will have to rescale the axes in matlab many times. I have compared the output from ItoTaylor tree to analytical solution for the SDE dX=mu0 * X(t)^beta *dt + sigma0 * X(t)^gamma *dz(t) Please note that analytical solution(from Schroeder's formula/its time changed version) given in the program can be compared with ItoTaylor solution only when beta=1 in the program. I have taken a very small step and used 64 time intervals per year. Green line in graph is ItoTaylor solution and red line is analytical solution. Sometimes, when parameters of SDE are extreme, you might want a ItoTaylor tree with first four hermite polynomials and you could then also take a larger step. Larger step might not have much of an advantage since we are working with a hundred or so paths. Almost the density of any SDE can be found by this method. I have written the following program for CEVlike diffusion and lognormallike drift because it can be analytically evaluated with Schroeder's formula. But you can find the density of almost any SDE with meanreversion and so many other features. You could also calculate the densities of path integrals of these stochastic volatility SDEs and plug those values in lognormal or CEV formula. Here is the matlab program. Please note that program execution speed involves the calculation of analytic density and I did not work on optimizing that. I have taken like 2000 integration steps while anybody could easily do it with only 200 steps.  function [] = SDEdensity() x0=.50; mu0=.15; sigma0=.55; beta=1.0; gamma=.85; dt=.0625/4; Tt=64*2; T=Tt*dt; Nn=81; IdtdzN=1/sqrt(3.0)*dt^1.5; IdzdtN=(11/sqrt(3.0))*dt^1.5; Idtdt0=dt^2/2.0; IdzdzN2=.5*dt; Idzdz0=.5*dt; Idt=dt; Idz=sqrt(dt); for nn=1:Nn X(nn)=x0; end for tt=1:Tt for nn=1:Nn zz(nn) = (nn*.1  4.0); HCoeff(nn,1)=mu0*X(nn)^beta*Idt+mu0^2.0*beta*X(nn)^(2*beta1)*Idtdt0+ ... sigma0^2*gamma*X(nn)^(2*gamma1)*Idzdz0+.5*beta*(beta1)*mu0*sigma0^2*X(nn)^(beta+2*gamma2)*Idtdt0; HCoeff(nn,2)=sigma0*X(nn)^gamma*Idz+mu0*sigma0*beta*X(nn)^(gamma+beta1)*IdzdtN+ ... sigma0*mu0*gamma*X(nn)^(gamma+beta1)*IdtdzN+ .5*gamma*(gamma1)*sigma0^3*X(nn)^(3*gamma2)*IdtdzN; HCoeff(nn,3)=sigma0^2*gamma*X(nn)^(2*gamma1)*IdzdzN2; X(nn)=X(nn)+HCoeff(nn,1)+HCoeff(nn,2)*1/sqrt(Tt)*zz(nn)+HCoeff(nn,3)*(1/sqrt(Tt)*zz(nn))^2.0; end end CFlag = 1; for nn = Nn  1:1:2 if((X(nn) > X(nn  1)) && (CFlag == 1)) nnStart = nn; else CFlag = 0; end end for nn=nnStart:Nn1 zz(nn) = (nn*.1  4.0); Dfx(nn)=(X(nn+1)X(nn1))/(zz(nn+1)zz(nn1)); end for nn=nnStart:Nn1 fx(nn)=normpdf(zz(nn),0,1)/Dfx(nn); end
Nn2 = 2000; gamma0 = gamma  1; X0 = x0; if(mu0~=0) tau=(exp(2*mu0*gamma0*T)  1)/(2*mu0*gamma0); else tau=T; end for n = 1:Nn2 yy(n) = (.0002 + n*.0050); yy1(n) = (.0002 + n*.005)*exp(mu0*T); fyy(n) = exp(mu0*T)*yy1(n)^(2*gamma0  3/2)*x0^(.5)/(sigma0^2*abs(gamma0)*tau)* ... besseli(abs(1/(2*gamma0)), (yy1(n)^(gamma0)*x0^(gamma0)/(sigma0^2*gamma0^2*tau)))* ... exp(1*(x0^(2*gamma0) + yy1(n)^(2*gamma0))/(2*sigma0^2*gamma0^2*tau)); end fyy; plot(X(nnStart:Nn1),fx(nnStart:Nn1),'g',yy(1:Nn2),fyy(1:Nn2),'r'); end

Here is an output of the matlab program from previous post. Here X0=1, mu0=.15, sigma0=.65, gamma=.6 and T=2.0, I have used just two hermite polynomials with 64 steps per year. Green line is ItoTaylor density and red line is analytical density.
Here is another comparison of ItoTaylor tree density and density from analytical formula. Here X0=1, mu0=.15, sigma0=.65, gamma=.8 and T=5.0 I have rescaled the graph to show the density area. And I have used just two hermite polynomials.
In the previous post, some graphs diverge from true solution when certain treepaths (or ODE paths) get closer to zero and this can be easily tackled by selectively decreasing the step size for ODEs(for the specific branch of ItoTaylor tree that has gotten close to zero) in the ItoTaylor tree when they reach close to zero. Paths diverge from true solution mostly for diffusions that have a CEV parameter gamma close to .5
I have also discussed the method in more detail at W****tt.com from post 703 to post 712. You can view that discussion here. https://forum.W****tt.com/viewtopic.php?f=4&t=99702&start=705





amin


Total Posts: 279 
Joined: Aug 2005 


Though I hope there will be interesting uses of the new ItoTaylor tree method I mentioned in the previous post, some new possibilities also open up with this method. For example, in SDEs literature, it is almost always assumed that noise scales as square root of time as in Brownian motion. Though this might be intelligent in most of the situations, there is nothing there to say that this must always be the case. In many SDEs modelling physical phenomenon, noise might as well scale as .25 power of time interval or possibly .75 power of time interval and we can very easily solve those SDEs with different time scaling with this new method. Time scaling must not necessarily be related to Brownian motion and when it is actually not, we can use the new method to solve other SDEs that arise with different time scaling in the random noise. I am sure many smart finance professionals would like to use this method with time scaling different from squreroot of t. 



amin


Total Posts: 279 
Joined: Aug 2005 


Here I made a post on linkedin about this new method. I have also mentioned nuclearphynance.com there. Here is the link to post. https://www.linkedin.com/pulse/convertingsdeautonomousodesitotaylortreemethodahsanamin/
Here I copy the post.
Converting An SDE into Autonomous ODEs: The ItoTaylor Tree Method for Solution of SDEs.
In my new research on the solution of SDEs and their stochastic integrals, I have worked with a new ItoTaylor tree method for calculation of densities of SDEs and their path integrals. In this method, we start an ItoTaylor tree for the SDE with a few branches and each branch is associated with a certain value of standard normal random variable and the tree branches just once at origin/starting value of SDE . Each branch of ItoTaylor tree is like an ODE for a particular value of the standard normal variable and each branch of the tree evolves in time autonomously and never rebranches save once at origin. And we recompute the ItoTaylor coefficients given the new initial values of X for every branch along the ItoTaylor tree for every time level. Basically, ItoTaylor tree method is just converting an SDE evolution into ODElike evolution where each ODE evolves autonomously with the particular value of normal assigned to the specific branch of ItoTaylor tree.
The method is lightning fast since the tree branches once at origin and never branches after that. In the accompanying matlab program, I have taken only 81 branches that span a standard normal between 4 and +4 with a normal gird spacing of .1 and there are only 64 time steps per year and one step computational effort for every branch is roughly comparable to monte carlo simulation computational effort for one single path. So you can yourself imagine how superfast the method is. Its speed is roughly comparable to monte carlo with only 81 paths.
I have posted an initial version of ItoTaylor tree matlab program based on only first two hermite polynomials that you can download. You can change the program in all sort of interesting ways. And you will have to rescale the axes in matlab many times. I have compared the output from ItoTaylor tree to analytical solution for the SDE dX=mu0 * X(t)^beta *dt + sigma0 * X(t)^gamma *dz(t) Please note that analytical solution(from Schroeder's formula/its time changed version) given in the program can be compared with ItoTaylor solution only when beta=1 in the program. I have taken a very small step and used 64 time intervals per year. Green line in graph is ItoTaylor solution and red line is analytical solution.
You can download and freely distribute the matlab program for calculation of density of SDEs using this new method. There is a lot of information I posted about the new method on W****tt.com math finance web site in past few days. I have posted a very simple matlab program that you can download from post 709 on W****tt.com mathematical finance forum here: https://forum.W****tt.com/viewtopic.php?f=4&t=99702&start=705 or from nuclearphynance.com forum here: http://nuclearphynance.com/Show%20Post.aspx?PostIDKey=180518
Sometimes, when parameters of SDE are extreme, you might want an ItoTaylor tree with first four hermite polynomials instead of just two that I used in the program and you could then also take a larger step. Larger step might not have much of an advantage since we are working with only a hundred or so paths. And sometimes, when stability issues arise, you could use the method with just one polynomial.
Almost the density of any SDE can be found by this method. I have written the following program for CEVlike diffusion and lognormallike drift because it can be analytically evaluated with Schroeder's formula. But you can find the density of almost any SDE with meanreversion and so many other features. You could also calculate the densities of path integrals of these stochastic volatility SDEs and plug those values in lognormal or CEV formula.
Here is a graphical comparison of the density from a simple matlab program and the analytical density. Here X0=1, mu0=.15, sigma0=.65, gamma=.6 and T=2.0, I have used just two hermite polynomials with 64 steps per year. Green line is ItoTaylor density and red line is analytical density. On the link below, you can view the graph comparison. https://s26.postimg.cc/cdnmh80y1/Matlab_Graph01.jpg
Here is another comparison of ItoTaylor tree density and density from analytical formula. Here X0=1, mu0=.15, sigma0=.65, gamma=.8 and T=5.0, I have rescaled the graph to show the density area. And I have used just two hermite polynomials. You can view the graph comparison on link below. https://s26.postimg.cc/gzjqprwtl/Matlab_Graph02.jpg









