
amin


Total Posts: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
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: 259 
Joined: Aug 2005 

 





