SERVING THE QUANTITATIVE FINANCE COMMUNITY

 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 10th, 2020, 7:05 am

Sorry, I was too sleepy last night to express myself clearly and I was being very incoherent. But the true idea that I had yesterday was that we should freeze the expectations made at earlier times and simply forward transport them along the Z-standard deviation grid associated with standard normal probability mass. By Z-standard deviation grid, I mean an expanding or contracting grid where grid cells are divided so that probability mass within each cell remains constant over time. So if there are payoffs at t1, t2, t3 and t4. We could freeze the expected payoffs at t1, on expanding/contracting standard deviations grid that preserves the probability mass in them and simply transport the expectation at t1 to respective cells on standard deviations grid at time t2 if we need to analyze the probability distribution at time t2 otherwise we could simply transport the frozen expectations to corresponding standard deviations grid at terminal time t4. 
This could work for many types of "payoffs". But it still remains to be seen if the idea would work for all "payoffs". But I am sure exploring it will help solve the most general problem.
I am trying to superimpose a standard deviations grid on regular fixed grid and trying to learn "calculus" to convert expectations across two grids.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 13th, 2020, 4:24 pm

Here is why the expectation on transition probabilities grid maintains its mean and its variance changes as the probability density evolves. First of all, let us suppose we have evolved the expectation using appropriate diffusion till time [$]t_{1}[$] or alternatively calculated it at [$]t_1[$] particular point in time using the grid, this expectation is given as [$]E[f(Y(t_{1}))]=f(Y(t_{1}))P(Y(t_{1}))[$]. This expectation stops changing even after it continues to evolve to step [$]t_2[$] after we do not add an appropriate diffusion generator for the next step. At time [$]t_{2}[$], this expectation changes to
[$]E[f(Y_1(t_{2}))]=\frac{E[f(Y(t_{1}))]}{P(Y(t_{2}))} P(Y(t_{2})) =\frac{f(Y(t_{1}))P(Y(t_{1}))}{P(Y(t_{2}))} P(Y(t_{2}))=f(Y(t_{1}))P(Y(t_{1}))[$]
Here I have used the term [$]E[f(Y_1(t_{2}))][$] for an expectation that was calculated at time [$]t_1[$] and was diffused/evolved to time [$]t_2[$] without its appropriate Ito-generator where we needed to calculate the expectation with probability [$]P(Y(t_{2}))[$]. As we saw in the above calculating the expectation with a different probability did not change its expectation. However repeating the above step for variance gives us the new equations as (supposing the mean is zero)
[$]Var[f(Y_1(t_{2}))]={[\frac{E[f(Y(t_{1}))]}{P(Y(t_{2}))}]}^2 P(Y(t_{2}))[$]
[$] ={\frac{[f(Y(t_{1}))P(Y(t_{1}))}{P(Y(t_{2}))}]}^2 P(Y(t_{2}))=f(Y(t_{1}))P(Y(t_{1}))[$]
[$]=\frac{[P(Y(t_{1}))}{P(Y(t_{2}))} [{f(Y(t_{1}))]}^2 P(Y(t_{1}))][$]
[$]=K(t_1,t_2)*Var(Y(t_1))[$]

where [$]K(t_1,t_2)=\frac{P(Y(t_{1}))}{P(Y(t_{2}))}[$] is the one step variance multiplier and  true variance is given as [$]Var(Y(t_1))={f(Y(t_{1}))}^2 P(Y(t_{1}))[$]
So the variance of a [$]t_1[$] expectation calculated with time [$]t_2[$] probabilities is given as
[$]Var[f(Y_1(t_{2}))]=K(t_1,t_2)*Var(Y(t_1))[$]
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 17th, 2020, 1:41 pm

Friends, there is more progress on the work regarding densities of arithmetic path integrals of arbitrary SDEs. In our ito-taylor expnasion of the volatility, some terms are of order sqrt(dt) other are of order dt and even other are of order dt^1.5. The noise of order sqrt(dt) is transported forward on step by multipliplying the noise with sqrt(t2/t1)=sqrt((t1+dt)/t1). while noise terms of order dt are transported forward by multiplication with (t2/t1)=((t1+dt)/t1) and noise terms of order dt^1.5 are transported forward by multiplication with (t2/t1)^1.5=((t1+dt)^1.5/(t1^1.5); When noises of different order are separately evolved, it is very easy to forward transport them by appropriate multiplication ratio for relevant "time order" of the noise. This was a very needed discovery to forward transport the "payoffs" or "data" on our transition probabilities grid for non-linear SDEs. Having done this, I am very confident that I will very soon be able to write and upload a complete program for the densities of arithmetic path integrals of functions of arbitrary SDEs. If you try this on your own, please know that sign on mean-reverting noise terms is negative.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 28th, 2020, 8:39 am

I have been able to concentrate on my research for past few days. While working on the arithmetic path integrals in transition probabilities framework project, I had some ideas and I was able to find a precise hermite polynomial representation for a general (existing) density in bessel coordinates. I found an algorithm using which we can easily convert any general bessel density(that has been arbitrarily evolved in time already) into a hermite polynomial representation. This was earlier not possible. When I found this algorithm, I digressed from the arithmetic path integrals project for a bit and decided to switch for a very few days towards the density evolution algorithm using SD fractions that expand or contract with time. Ok, here is the good news that I was able to refine my earlier algorithm and make it quite sharply precise for the evolution of densities of the SDEs in SD fractions grid(moving and expanding grid so as to keep probability mass in each SD fraction constant) framework. Here are the features of the new algorithm.
1. It is far more accurate than the earlier version. The accuracy is excellent everywhere.

2. It is very accurate for mean reverting SDEs  with smaller values of gamma when SDE mean is very close to zero where the earlier versions were being very inaccurate. For example with gamma=.75 or gamma=.65, the density was far too peaked as compared to true density when SDE mean was close to zero. This problem has been perfectly corrected in the new version since we are using an exact hermite representation of the bessel density.

3. The same density evolution algorithm works for mean reverting, partially mean-reverting,  non-mean reverting, zero mean and all other types of SDE densities. There are no arbitrary multiplications with parameters like gamma or kappa like I had in the previous version. I am sure that algorithm would work for all possible SDEs that are well behaved whether mean-reverting or partially mean-reverting or otherwise.

I am testing the algorithm for all possible errors and will be posting it in another two to three days. This new program became possible only because I was able to find a hermite representation of the arbitrary SDE density in bessel coordinates. I will be psoting the code here with explanations very soon. 
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 30th, 2020, 2:18 pm

Dear friends, in this current post, I want to describe how to convert an existing density that is in bessel coordinates into hermite polynomial representation. There could be several ways of doing this. One possible way could follow the route using the concept of traditional standard deviations but this way does not work everywhere in stochastics when we are dealing with expansions of SDEs. As we learnt earlier when we were expanding to find higher order terms in the expansion of SDEs, all the expansions have to be done in terms of the standard gaussian and hermite polynomials. So the ideas I am going to present are extremely simple but could be quite far-reaching for the stochastics, therefore I decided to formally explain the ideas in this post.
Suppose we are given an existing density of an SDE w that is in bessel coordinates and we want its expansion in terms of hermite polynomials. We suppose that the density of w is spread on N discrete grid points. The value of w on nth grid point is [$]w(n)[$] and the associated value of standard normal on that nth grid point is [$]Z(n)[$]. We associate a value of volatility [$]\sigma_{0}(n)[$] with each of the n grid points. This  [$]\sigma_{0}(n)[$] is calculated simply as
[$] \sigma_{0}(n)=\frac{(w(n)-w_{mean})}{(Z(n)-Z_{mean})}[$]   Equation(1)

Here [$]Z_{mean}[$] is the value of Z mapped on [$]w_{mean}[$]
So we can find a first order representation of the SDE variable w as
[$]w(n)=w_{mean}- \sigma_{0}(n) Z_{mean}+ \sigma_{0}(n) Z(n)[$]    Equation(2)
In the expansion of bessel SDEs, the mean is very close to median where median is associated with [$]Z(n)=0[$]. We have to do some extra calculations to account for the difference between expansion from median or mean. I want to just highlight the expansion method so we suppose for now that difference between [$]Z_{Mean}[$] and [$]Z_{median}=Z(n_{Median})=0[$] is negligible. We will over come this later.
So we suppose we can write 
[$]w(n)=w_{mean}+ \sigma_{0}(n)  Z(n)[$]               Equation(3)
but we can also expand state dependent [$]\sigma_{0}(n)[$] around the median expansion point as
[$]\sigma_{0}(n)=\sigma_{0}(n_{Median})+\int_{Z(n_{Median})}^{Z(n)} \frac{d\sigma_{0}(n)}{dZ(n)} dZ(n)[$]  Equation(4)
which can be written as
[$]\sigma_{0}(n)=\sigma_{0}(n_{Median})+\sigma_1(n) Z(n)[$]  Equation(5)
where [$]\sigma_1(n)=\frac{\int_{Z(n_{Median})}^{Z(n)} \frac{d\sigma_{0}(n)}{dZ(n)} dZ(n)}{\int_{Z(n_{Median})}^{Z(n)} dZ(n)}[$]
In what follows, we use notation [$]\sigma_{00}= \sigma_{0}(n_{Median})[$] for the constant value.
using above equation(5) in equation(3), we can write after a few calculations as
[$]w(n)=w_{mean}+ \sigma_{00}  Z(n) + \sigma_1(n) (Z(n)^2-1)[$] Equation(6)

Now similarly to previous calculations, we expand state dependent [$]\sigma_{1}(n)[$] around the median expansion point as
[$]\sigma_{1}(n)=\sigma_{1}(n_{Median})+\int_{Z(n_{Median})}^{Z(n)} \frac{d\sigma_{1}(n)}{dZ(n)} dZ(n)[$]  Equation(7)

which can be written as
[$]\sigma_{1}(n)=\sigma_{1}(n_{Median})+\sigma_2(n) Z(n)[$]  Equation(8)
where [$]\sigma_2(n)=\frac{\int_{Z(n_{Median})}^{Z(n)} \frac{d\sigma_{1}(n)}{dZ(n)} dZ(n)}{\int_{Z(n_{Median})}^{Z(n)} dZ(n)}[$]
using the notation [$]\sigma_{11}= \sigma_{1}(n_{Median})[$] we can write after forcing equation(8) in equation (6) as

[$]w(n)=w_{mean}+ \sigma_{00}  Z(n) + \sigma_{11} (Z(n)^2-1)+ \sigma_2(n) (Z(n)^3-3Z(n))[$] Equation(9)
So if we can continue to expand the density into higher hermite polynomials. However the procedure would require care for accuracy in numerical discrete differentiation and integrations and other errors when values associated with higher hermite polynomials become incorrect garbage.

Ok, I have not fully developed the ideas yet and these are very recent thoughts but I decided to share them with friends here to highlight the importance of above expansions.
Let us suppose, we have an SDE of the form in bessel coordinates given as
[$]dw(t)=\mu(w(t))dt + \sigma dz(t) [$]
we can expand both [$]w(t)[$] and the term [$]\mu(w(t))dt[$] in terms of hermite polynomials using the above hermite representation procedure and we could write for example as
[$]w(t)=\mu_{w}+ a0* Z(n) + a1* ((Z(n)^2-1)[$]
[$]+ a2* ((Z(n)^3-3 Z(n)) + a3(n)* (Z(n)^4-6 Z(n)^2+3) [$] 
and similarly for  
[$]drift(t)=\mu_{drift}+ b0* Z(n) + b1* ((Z(n)^2-1)[$]
[$]+ b2* ((Z(n)^3-3 Z(n)) + b3(n)* (Z(n)^4-6 Z(n)^2+3) [$]  
where drift(t) represents the term [$]\mu(w(t))dt[$]
and use some straightforward way to add across the orthogonal coordinates. Of course there are going to be some slight complexities but I am sure we can find some quite simple ways of advancing the densities very soon.
I will also, possibly tomorrow, share the very simple program that converts the densities into a hermite polynomial representation. Over the weekend, I will be sharing a lot of programs also about the SDE density evolution in moving and expanding constant mass SD fractions. But I think the ideas presented in the above post would have very interesting implications when adding the two arbitrary SDE densities. Apart from so many other possible applications, we could use them in hybrid fX/equity/IR models where so many different SDEs are in action together.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

January 31st, 2020, 5:25 pm

This program find hermite representation of the SDE variable w to Second Order as a function of two hermite polynomials. I will soon be posting another program in a few days that will find the hermite representation to arbitrary order specified by the user with interesting applications. Lot of interesting stuff is due in next few days. This program should become easily familiar for friends who have been downloading the earlier programs and understand them.
I expanded the SDE variable w around the mean as mean is usually very convenient but you can expand around median or any arbitrary point on the SDE grid.
 
function [wMean,Sigma00,Sigma1,Sigma0] = CalculateHermiteRepresentationOfDensityOrderTwo(w,Z,Prob,wnStart,Nn,dNn)

%This program find hermite representation of the SDE variable w to 
%Second Order as a function of two hermite polynomials.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite 
%polynomials. It is not the density which can however be found by change of variable 
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
%  Z(wnStart:Nn) represents the associated value of standard Gaussian on
%  the grid. w and Z can be associated by matching appropriate CDF points
%  on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and 
% Z(n) and Z(n+1) two points on the grid.
%Prob is the probability mass associated with each point on the grid.

wMean=sum(w(1:Nn).*Prob(1:Nn));  %Calculate the mean. This is in general 
%not on the discrete grid of w. 

%Below calculate the associated value of ZMean that maps wMean to ZMean.
%Also the grid point index nnMean which is associated with the largest
% w(nn) such that w(nn)<=wMean.
for nn=1:Nn-1
    if( (w(nn+1)>wMean) && (w(nn)<=wMean))
        ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
        nnMean=nn;
    end
end

%Below Calculate the vol associated with each point on the grid
for nn=1:Nn
    Sigma0(nn)=(w(nn)-wMean)./(Z(nn)-ZMean);  
 end

%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial. We chose this volatility equal to
%volatility at mean since we are expanding around the mean.
%If mean has zero volatility, we will have to find volatility at mean by
%interpolating the volatility at one higher grid point and one lower grid
%point. Volatility at mean must not be zero (unless we have a delta density).
Sigma00=Sigma0(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma0(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);

%Following is the operation related to equation(4) in post 860.     
SSigma1(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end

%Below is the averaging
Sigma1(wnStart:nnMean-1)=SSigma1(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma1(nnMean+1:Nn)=SSigma1(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma1(nnMean)=(Sigma1(nnMean-1)+Sigma1(nnMean+1))/2.0;

%Now below we represent the SDE variable w in terms of first two hermite
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
   
 w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
    
%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
 plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1(wnStart:Nn),'r')
str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 1st, 2020, 1:06 pm

I know friends would have liked if I ran the above code with a real example (I was thinking and I realized that). Hopefully tomorrow I will be uploading my another program for evolution of density of SDEs with moving and expanding SD fractions framework that I have been talking about earlier in post 859. I have used the above program in post 861 on every time step in that program that I will post tomorrow. So friends would be able to test the above program with any possible SDE they like. Please bear with me for another day. 
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 2nd, 2020, 9:39 am

In this program, I am simulating the SDE given as
dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

I have not directly simulated the SDE but simulated the transformed Bessel process version of the SDE and then changed coordinates to retreive
the SDE in original coordinates.
In this version, the mean correction has been disabled. If you are simulating standard mean-reverting SV type SDEs, please enable the mean 
correction by uncommenting the appropriate line in the body of code.
you can specify any general mu1 and mu2 and beta1 and beta2. The present program is not perfect but quite very good.
if there is only one term in drift with beta=0, there will be minor errors. This case is yet not covered.
function [] = ItoHermiteMethodWilmott08Improved()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*1;
TtM=Tt/1;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

x0=.35;   % starting value of SDE
beta1=0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.850;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.250;   %mean reversion parameter.
theta=.1250;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.750;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
    %Above calculate probability mass in each probability subdivision.
end

 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;


tic

for tt=1:Tt
  
    yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

    theta1=mu1;
    theta2=mu2;
    theta3=.5*sigma0^2*(-gamma);
    omega1=beta1-gamma;
    omega2=beta2-gamma;
    omega3=gamma-1;
    
    GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
    dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
        omega3*theta3*yy(wnStart:Nn).^(omega3-1);
    DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
    dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
        theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
        theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
    QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
    ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
    d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
    ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
    d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);
    
    VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
    d_dGG_VV_dyy(wnStart:Nn)= ...
        sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2) + ...
        sigma0.*omega2.*theta2.*(omega2+gamma-1).*yy(wnStart:Nn).^(omega2+gamma-2) + ...
        sigma0.*omega3.*theta3.*(omega3+gamma-1).*yy(wnStart:Nn).^(omega3+gamma-2);    

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this code block, we use second hermite representation of the density
%which is used to solve for mean reversion.
%It is the value of Sigma1 below whcih is used later in the main SDE one
%time step update formula.

    if(tt>1)
        [wMean,ZMean,Sigma00,Sigma1,Sigma0] = CalculateHermiteRepresentationOfDensityOrderTwo(w,Z,ZProb,wnStart,Nn,dNn);    
   
        %w1New(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
        %    Sigma00*Z(wnStart:Nn)+ ...
        %    Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
        %Uncomment to see a graph comparison of w constructed from hermite
        %representation and the original w.
        %  plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1New(wnStart:Nn),'r')
        %  str=input('Look at the graphs--1');
    else
        Sigma1(wnStart:Nn)=0.0;
    end
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %In this block we determine whether the SDE is mean reverting.
   %This is done, possibly naiively, by taking a derivative of drift at 
   %value and if the derivative is negative, we calculate KappaMR which is 
   %parameter that determines the strength of mean reversion. KappaMR is
   %used in the main SDE one step update formula.
   
    
  
    yyMean=sum(ZProb(wnStart:Nn).*yy(wnStart:Nn));
    dDriftAtMean=mu1*beta1*yyMean.^(beta1-1)+mu2*beta2*yyMean.^(beta2-1);
    KappaMR(wnStart:Nn)=0;
    if(dDriftAtMean<0)
        if(beta1~=0)
            KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)-mu1*yy(wnStart:Nn).^(beta1-1);
        end
        if(beta2~=0)
            KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)-mu2*yy(wnStart:Nn).^(beta2-1);
        end
        KappaMR(wnStart:Nn)=abs(KappaMR(wnStart:Nn));
    end
    if(tt<2)
        KappaMR(wnStart:Nn)=0.0;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the main SDE variable one step update equation. 


 w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
          1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
    CTheta=1.0;
    D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
    D2wl=(35/12*w(NnMidl-2)	-26/3*w(NnMidl-3)+	19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
    wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
    wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
    wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
    wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
    %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%The last equation is not used since point w(NnMidh+3)is given but it is used in 
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.

    D1wh=(-11/6*w(NnMidh+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2)	-26/3*w(NnMidh+3)+	19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
%     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;


    whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
    whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
    whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;

    %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;


    w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
    w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
    w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
    w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
    w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

        Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
        Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        Pn=1.0;%%Sum(ZProb(1:Nn))
        Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);
    
        Y0Pn=Y0*Pn;
        W0=1/(YnPn-Y0Pn);
        YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
        wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
   %     w(wnStart:Nn)=wCorrected(wnStart:Nn);
        
   %I have disabled the mean correction. The above mean correction is only valid for 
   %standard mean reverting stochastic volatility type SDEs. To enable mean
   %correction please uncomment the above last line in the block.
        
        
end

 y_w(1:Nn)=0;
 y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
 Dfy_w(wnStart:Nn)=0;
 Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
    fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
    fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc
str=input('Analytic Values calculated; Press a key to start monte carlo');

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.
 rng(85109634, 'twister')
 paths=200000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=500;
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here are recent dependencies of the program.
function [wMean,ZMean,Sigma00,Sigma1,Sigma0] = CalculateHermiteRepresentationOfDensityOrderTwo(w,Z,Prob,wnStart,Nn,dNn)

%This program find hermite representation of the SDE variable w to 
%Second Order as a function of two hermite polynomials.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite 
%polynomials. It is not the density which can however be found by change of variable 
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
%  Z(wnStart:Nn) represents the associated value of standard Gaussian on
%  the grid. w and Z can be associated by matching appropriate CDF points
%  on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and 
% Z(n) and Z(n+1) two points on the grid.
%Prob is the probability mass associated with each point on the grid.

wMean=sum(w(1:Nn).*Prob(1:Nn))  %Calculate the mean. This is in general 
%not on the discrete grid of w. 

%Below calculate the associated value of ZMean that maps wMean to ZMean.
%Also the grid point index nnMean which is associated with the largest
% w(nn) such that w(nn)<=wMean.
for nn=1:Nn-1
    if( (w(nn+1)>wMean) && (w(nn)<=wMean))
        ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
        nnMean=nn;
    end
end

%Below Calculate the vol associated with each point on the grid
for nn=1:Nn
    Sigma0(nn)=(w(nn)-wMean)./(Z(nn)-ZMean);  
 end

%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial. We chose this volatility equal to
%volatility at mean since we are expanding around the mean.
%If mean has zero volatility, we will have to find volatility at mean by
%interpolating the volatility at one higher grid point and one lower grid
%point. Volatility at mean must not be zero (unless we have a delta density).
Sigma00=Sigma0(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma0(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);

%Following is the operation related to equation(4) in post 860.     
SSigma1(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end

%Below is the averaging
Sigma1(wnStart:nnMean-1)=SSigma1(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma1(nnMean+1:Nn)=SSigma1(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma1(nnMean)=(Sigma1(nnMean-1)+Sigma1(nnMean+1))/2.0;

%Now below we represent the SDE variable w in terms of first two hermite
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
   
 w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
    
%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
% plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1(wnStart:Nn),'r')
%str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2nd recent dependency is given here for graphing
function [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(X,Paths,NoOfBins,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(X(p)<0)
%X(p)=0;
%end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

%IndexMax=NoOfBins+1;
BinSize=(Xmax-Xmin)/NoOfBins;
%IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1

XDensity(1:IndexMax)=0.0;

for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);
if(real(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

Other dependencies have been already been uploaded in the previous code distribution posts. For example you could look at recent post 850 for other monte carlo dependencies.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 5th, 2020, 6:01 am

I have been able to drastically improve the previous version of the Ito-Hermite density evolution algorithm that I posted above a few days ago. I have added the third hermite polynomial correction to the previous program with just second hermite correction and made a few other minor changes and it works in an excellent fashion. The new program is just very close to perfect. I will be posting the program later tonight or tomorrow after checking it for some more time. It now works well for all exponents in drift equation. I am looking forward to sharing the new program by tomorrow by the latest.
For those friends who would want to ask how I calculated KappaMR in the previous program, I will refer them to post 788 in the same thread. Here I used kappaMR not with the bessel coordinates but with the original coordinates. It works for all equations mean-reverting or non mean reverting after we have added the third order skew effect hrermite polynomial. It was not working well for some exponents since for those exponents the effect of skew which was far greater and had not been added. In the updated new program that I will post tomorrow, I am using the same KappaMR with third hermite polynomial.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 6th, 2020, 4:46 am

Here is the new program with third hermite polynomial effect added. All the instructions for program in post 863 follow. This program is quite close to perfect for vols less than 1.0 i.e sigma0<=1.0. I expect there will be further some minor changes in the program after all the calculations are properly done with the "mean reverting exponential" (misnomer?) as I calculated in post 788 (but with original coordinates) and the orthogonal hermite polynomials. I have not made any precise calculations and basically just extended the previous code so I expect there will be further some minor changes in next few days or possibly a week. Here is the new version with updated dependency function that also calculates decomposition of SDE variable up to third orthogonal polynomial.
In the next version, I expect to add the following things:
1. A proper mean correction for (almost) every possible SDE,
2. Try to find a better solution for instability in the mid of the density.
3. Make any changes as theoretical calculations on "mean reverting exponential" and orthogonal hermite polynomials are done and become available.

function [] = ItoHermiteMethodWilmott08Improved02()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*2;
TtM=Tt/2;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.80;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.0;   %mean reversion parameter.
theta=.250;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.95;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
    %Above calculate probability mass in each probability subdivision.
end

 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;


tic

for tt=1:Tt
  
    yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

    theta1=mu1;
    theta2=mu2;
    theta3=.5*sigma0^2*(-gamma);
    omega1=beta1-gamma;
    omega2=beta2-gamma;
    omega3=gamma-1;
    
    GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
    dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
        omega3*theta3*yy(wnStart:Nn).^(omega3-1);
    DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
    dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
        theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
        theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
    QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
    ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
    d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
    ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
    d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);
    
    VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
    d_dGG_VV_dyy(wnStart:Nn)= ...
        sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2) + ...
        sigma0.*omega2.*theta2.*(omega2+gamma-1).*yy(wnStart:Nn).^(omega2+gamma-2) + ...
        sigma0.*omega3.*theta3.*(omega3+gamma-1).*yy(wnStart:Nn).^(omega3+gamma-2);    

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this code block, we use second hermite representation of the density
%which is used to solve for mean reversion.
%It is the value of Sigma1 below whcih is used later in the main SDE one
%time step update formula.

    if(tt>1)
        [wMean,ZMean,Sigma00,Sigma11,Sigma1,Sigma2,Sigma0,w1New,w2New] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);

      %  tt
    else
        Sigma1(wnStart:Nn)=0.0;
        Sigma2(wnStart:Nn)=0.0;
        Sigma11=0.0;
    end
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %In this block we determine whether the SDE is mean reverting.
   %This is done, possibly naiively, by taking a derivative of drift at 
   %value and if the derivative is negative, we calculate KappaMR which is 
   %parameter that determines the strength of mean reversion. KappaMR is
   %used in the main SDE one step update formula.
   
    
  
  %  yyMean=sum(ZProb(wnStart:Nn).*yy(wnStart:Nn));
  %  dDriftAtMean=mu1*beta1*yyMean.^(beta1-1)+mu2*beta2*yyMean.^(beta2-1);
    KappaMR(wnStart:Nn)=0;
    %if(dDriftAtMean<0)
       % if(beta1~=0)
            KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu1*yy(wnStart:Nn).^(beta1-1);
      %  end
      %  if(beta2~=0)
            KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu2*yy(wnStart:Nn).^(beta2-1);
      %  end
   %     SignKappaMR(wnStart:Nn)=sign(mu1*beta1*yy(wnStart:Nn).^(beta1-1)+mu2*beta2*yy(wnStart:Nn).^(beta2-1));
        KappaMR(wnStart:Nn)=abs(KappaMR(wnStart:Nn));
    %end
    if(tt<2)
        KappaMR(wnStart:Nn)=0.0;
        %SignKappaMR(wnStart:Nn)=1.0;
    end
    %KappaMR(wnStart:Nn)
    %SignKappaMR(wnStart:Nn)
    %str=input('Look at KappaMR');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the main SDE variable one step update equation. 


 w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
          sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+...
          .5*SignKappaMR(wnStart:Nn).*(2*KappaMR(wnStart:Nn)).^1.0.*c1(wnStart:Nn).^3.*Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));




%  w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
%           sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
%           sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
%           sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
%           1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*sigma0.^2.*dt.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
    CTheta=1.0;
    D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
    D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
    wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
    wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
    wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
    wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
    %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%The last equation is not used since point w(NnMidh+3)is given but it is used in 
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.

    D1wh=(-11/6*w(NnMidh+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2) -26/3*w(NnMidh+3)+ 19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
%     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;


    whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
    whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
    whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;

    %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;


    w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
    w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
    w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
    w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
    w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

        Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
        Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        Pn=1.0;%%Sum(ZProb(1:Nn))
        Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);
    
        Y0Pn=Y0*Pn;
        W0=1/(YnPn-Y0Pn);
        YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
        wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
    %    w(wnStart:Nn)=wCorrected(wnStart:Nn);
        
   %I have disabled the mean correction. The above mean correction is only valid for 
   %standard mean reverting stochastic volatility type SDEs. To enable mean
   %correction please uncomment the above last line in the block.
        
        
end

 y_w(1:Nn)=0;
 y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
 Dfy_w(wnStart:Nn)=0;
 Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
    fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
    fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc
str=input('Analytic Values calculated; Press a key to start monte carlo');

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.
 rng(85109634, 'twister')
 paths=200000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=500*gamma^2/.25*sigma0;
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is the new dependency function code.
function [wMean,ZMean,Sigma00,Sigma11,Sigma1,Sigma2,Sigma0,w1,w2] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,Prob,wnStart,Nn,dNn,NnMidh,NnMidl)

%This program find hermite representation of the SDE variable w to 
%tjird Order as a function of three hermite polynomials.
%This program is not perfectly precise and has very small errors so you
%will have to be careful using it everywhere.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite 
%polynomials. It is not the density which can however be found by change of variable 
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
%  Z(wnStart:Nn) represents the associated value of standard Gaussian on
%  the grid. w and Z can be associated by matching appropriate CDF points
%  on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and 
% Z(n) and Z(n+1) two points on the grid.
%Prob is the probability mass associated with each point on the grid.

wMean=sum(w(1:Nn).*Prob(1:Nn))  %Calculate the mean. This is in general 
%not on the discrete grid of w. 

%Below calculate the associated value of ZMean that maps wMean to ZMean.
%Also the grid point index nnMean which is associated with the largest
% w(nn) such that w(nn)<=wMean.
ZMean=(w(NnMidh)+w(NnMidl))/2.0;
nnMean=NnMidl;
for nn=1:Nn-1
    if( (w(nn+1)>wMean) && (w(nn)<=wMean))
        ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
        nnMean=nn;
    end
end

%Below Calculate the vol associated with each point on the grid
for nn=1:Nn
    Sigma0(nn)=(w(nn)-wMean)./(Z(nn)-ZMean);  
 end

%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial. We chose this volatility equal to
%volatility at mean since we are expanding around the mean.
%If mean has zero volatility, we will have to find volatility at mean by
%interpolating the volatility at one higher grid point and one lower grid
%point. Volatility at mean must not be zero (unless we have a delta density).
Sigma00=Sigma0(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma0(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);

%Following is the operation related to equation(4) in post 860.     
SSigma1(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end

%Below is the averaging
Sigma1(wnStart:nnMean-1)=SSigma1(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma1(nnMean+1:Nn)=SSigma1(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma1(nnMean)=(Sigma1(nnMean-1)+Sigma1(nnMean+1))/2.0;


Sigma11=Sigma1(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma1(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));



%Now below we represent the SDE variable w in terms of first two hermite
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
   
 w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
 


dSigma1dZ(wnStart)=(-3*Sigma1(wnStart)+4*Sigma1(wnStart+1)-1*Sigma1(wnStart+2))/(2*dNn);
dSigma1dZ(wnStart+1:Nn-1)=(-1*Sigma1(wnStart:Nn-2)+1*Sigma1(wnStart+2:Nn))/(2*dNn);
dSigma1dZ(Nn)=(3*Sigma1(Nn)-4*Sigma1(Nn-1)+1*Sigma1(Nn-2))/(2*dNn);

% Sigma2(NnMidh)=dSigma1dZ(NnMidh);
% for nn=NnMidh+1:Nn
%     Sigma2(nn)=(Sigma2(nn-1)*Z(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end
% Sigma2(NnMidl)=dSigma1dZ(NnMidl);
% for nn=NnMidl-1:-1:wnStart
%     Sigma2(nn)=(Sigma2(nn+1)*Z(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end

SSigma2(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma2(nn)=SSigma2(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma2(nn)=SSigma2(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn;
end

Sigma2(wnStart:nnMean-1)=SSigma2(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma2(nnMean+1:Nn)=SSigma2(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma2(nnMean)=(Sigma2(nnMean-1)+Sigma2(nnMean+1))/2.0;



w2(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma11.*(ZMean.^2-1)+ ...
 3*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) - ...
 3*Sigma2(wnStart:Nn).*Z(wnStart:Nn).*(ZMean.^2-1) + ... 
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma11.*(Z(wnStart:Nn).^2-1)+ ...   
 0*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) + ...
 Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
 



%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
 plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1(wnStart:Nn),'r',Z(wnStart:Nn),w2(wnStart:Nn),'b')
Sigma00
 str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 7th, 2020, 5:27 am

Sorry friends, there was an inadvertent mistake that stopped the program from running properly. I noticed it only today. I am re-posting the corrected code. Here is the corrected code and its function. 
function [] = ItoHermiteMethodWilmott08Improved02Corrected()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*2;
TtM=Tt/2;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.80;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.0;   %mean reversion parameter.
theta=.250;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.95;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
    %Above calculate probability mass in each probability subdivision.
end

 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;


tic

for tt=1:Tt
  
    yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

    theta1=mu1;
    theta2=mu2;
    theta3=.5*sigma0^2*(-gamma);
    omega1=beta1-gamma;
    omega2=beta2-gamma;
    omega3=gamma-1;
    
    GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
    dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
        omega3*theta3*yy(wnStart:Nn).^(omega3-1);
    DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
    dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
        theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
        theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
    QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
    ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
    d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
    ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
    d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);
    
    VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
    d_dGG_VV_dyy(wnStart:Nn)= ...
        sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2) + ...
        sigma0.*omega2.*theta2.*(omega2+gamma-1).*yy(wnStart:Nn).^(omega2+gamma-2) + ...
        sigma0.*omega3.*theta3.*(omega3+gamma-1).*yy(wnStart:Nn).^(omega3+gamma-2);    

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this code block, we use second hermite representation of the density
%which is used to solve for mean reversion.
%It is the value of Sigma1 below whcih is used later in the main SDE one
%time step update formula.

    if(tt>1)
        [wMean,ZMean,Sigma00,Sigma11,Sigma1,Sigma2,Sigma0,w1New,w2New] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);

      %  tt
    else
        Sigma1(wnStart:Nn)=0.0;
        Sigma2(wnStart:Nn)=0.0;
        Sigma11=0.0;
    end
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %In this block we determine whether the SDE is mean reverting.
   %This is done, possibly naiively, by taking a derivative of drift at 
   %value and if the derivative is negative, we calculate KappaMR which is 
   %parameter that determines the strength of mean reversion. KappaMR is
   %used in the main SDE one step update formula.
   
    
    KappaMR(wnStart:Nn)=0;
    KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu1*yy(wnStart:Nn).^(beta1-1);
    KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu2*yy(wnStart:Nn).^(beta2-1);
    KappaMR(wnStart:Nn)=abs(KappaMR(wnStart:Nn));
    %end
    if(tt<2)
        KappaMR(wnStart:Nn)=0.0;
        %SignKappaMR(wnStart:Nn)=1.0;
    end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the main SDE variable one step update equation. 


 w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
          sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+...
          .5.*(2*KappaMR(wnStart:Nn)).^1.0.*c1(wnStart:Nn).^3.*Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));




%  w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
%           sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
%           sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
%           sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
%           1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*sigma0.^2.*dt.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
    CTheta=1.0;
    D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
    D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
    wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
    wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
    wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
    wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
    %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%The last equation is not used since point w(NnMidh+3)is given but it is used in 
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.

    D1wh=(-11/6*w(NnMidh+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2) -26/3*w(NnMidh+3)+ 19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
%     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;


    whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
    whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
    whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;

    %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;


    w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
    w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
    w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
    w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
    w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

        Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
        Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        Pn=1.0;%%Sum(ZProb(1:Nn))
        Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);
    
        Y0Pn=Y0*Pn;
        W0=1/(YnPn-Y0Pn);
        YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
        wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
    %    w(wnStart:Nn)=wCorrected(wnStart:Nn);
        
   %I have disabled the mean correction. The above mean correction is only valid for 
   %standard mean reverting stochastic volatility type SDEs. To enable mean
   %correction please uncomment the above last line in the block.
        
        
end

 y_w(1:Nn)=0;
 y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
 Dfy_w(wnStart:Nn)=0;
 Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
    fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
    fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc
str=input('Analytic Values calculated; Press a key to start monte carlo');

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.
 rng(85109634, 'twister')
 paths=200000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=500*gamma^2/.25*sigma0;
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
function [wMean,ZMean,Sigma00,Sigma11,Sigma1,Sigma2,Sigma0,w1,w2] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,Prob,wnStart,Nn,dNn,NnMidh,NnMidl)

%This program find hermite representation of the SDE variable w to 
%tjird Order as a function of three hermite polynomials.
%This program is not perfectly precise and has very small errors so you
%will have to be careful using it everywhere.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite 
%polynomials. It is not the density which can however be found by change of variable 
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
%  Z(wnStart:Nn) represents the associated value of standard Gaussian on
%  the grid. w and Z can be associated by matching appropriate CDF points
%  on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and 
% Z(n) and Z(n+1) two points on the grid.
%Prob is the probability mass associated with each point on the grid.

wMean=sum(w(1:Nn).*Prob(1:Nn))  %Calculate the mean. This is in general 
%not on the discrete grid of w. 

%Below calculate the associated value of ZMean that maps wMean to ZMean.
%Also the grid point index nnMean which is associated with the largest
% w(nn) such that w(nn)<=wMean.
ZMean=(w(NnMidh)+w(NnMidl))/2.0;
nnMean=NnMidl;
for nn=1:Nn-1
    if( (w(nn+1)>wMean) && (w(nn)<=wMean))
        ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
        nnMean=nn;
    end
end

%Below Calculate the vol associated with each point on the grid
for nn=1:Nn
    Sigma0(nn)=(w(nn)-wMean)./(Z(nn)-ZMean);  
 end

%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial. We chose this volatility equal to
%volatility at mean since we are expanding around the mean.
%If mean has zero volatility, we will have to find volatility at mean by
%interpolating the volatility at one higher grid point and one lower grid
%point. Volatility at mean must not be zero (unless we have a delta density).
Sigma00=Sigma0(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma0(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);

%Following is the operation related to equation(4) in post 860.     
SSigma1(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end

%Below is the averaging
Sigma1(wnStart:nnMean-1)=SSigma1(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma1(nnMean+1:Nn)=SSigma1(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma1(nnMean)=(Sigma1(nnMean-1)+Sigma1(nnMean+1))/2.0;


Sigma11=Sigma1(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma1(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));



%Now below we represent the SDE variable w in terms of first two hermite
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
   
 w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
 


dSigma1dZ(wnStart)=(-3*Sigma1(wnStart)+4*Sigma1(wnStart+1)-1*Sigma1(wnStart+2))/(2*dNn);
dSigma1dZ(wnStart+1:Nn-1)=(-1*Sigma1(wnStart:Nn-2)+1*Sigma1(wnStart+2:Nn))/(2*dNn);
dSigma1dZ(Nn)=(3*Sigma1(Nn)-4*Sigma1(Nn-1)+1*Sigma1(Nn-2))/(2*dNn);

% Sigma2(NnMidh)=dSigma1dZ(NnMidh);
% for nn=NnMidh+1:Nn
%     Sigma2(nn)=(Sigma2(nn-1)*Z(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end
% Sigma2(NnMidl)=dSigma1dZ(NnMidl);
% for nn=NnMidl-1:-1:wnStart
%     Sigma2(nn)=(Sigma2(nn+1)*Z(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end

SSigma2(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma2(nn)=SSigma2(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma2(nn)=SSigma2(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn;
end

Sigma2(wnStart:nnMean-1)=SSigma2(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma2(nnMean+1:Nn)=SSigma2(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma2(nnMean)=(Sigma2(nnMean-1)+Sigma2(nnMean+1))/2.0;



w2(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma11.*(ZMean.^2-1)+ ...
 3*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) - ...
 3*Sigma2(wnStart:Nn).*Z(wnStart:Nn).*(ZMean.^2-1) + ... 
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma11.*(Z(wnStart:Nn).^2-1)+ ...   
 0*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) + ...
 Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
 



%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
% plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1(wnStart:Nn),'r',Z(wnStart:Nn),w2(wnStart:Nn),'b')
%Sigma00
% str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end

 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 7th, 2020, 9:21 am

I did some calculations and also played around with the models and I believe that we have to add the effect of "Mean Reverting exponential" on the existing w(t) curve and it is not to be applied on the noise part(This correction is not about the effect of exponential on the noise integral part of the equation). So if you replace c1(wnStart:Nn).^2 by dt and c1(wnStart:Nn).^3 by dt^1.5 in the main update equation where they are multiplying with the respective orthogonal polynomials, you will get far better results for every volatility level. This is first update and their could be further changes but making this change will hugely improve results at all volatility levels. Just for clarity here c1(.)=sigma(.)*sqrt(dt). we do not want to include sigma(.) in the calculations as this effect is about the curve and not the noise integral part.
Here is what I believe that main update equation should be to third order

w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
          sqrt(2*KappaMR(wnStart:Nn)).*(dt).*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+...
          sqrt(2/3).*(KappaMR(wnStart:Nn)).^1.0.*dt.^1.5.*Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));

I will be posting an update program in a few days but friends can make this correction and get far better results.

Edit: I am attaching the file with above update included.
function [] = ItoHermiteMethodWilmott08Improved02Corrected()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*2;
TtM=Tt/2;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.750;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.50;   %mean reversion parameter.
theta=.150;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.850;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
    %Above calculate probability mass in each probability subdivision.
end

 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;


tic

for tt=1:Tt
  
    yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

    theta1=mu1;
    theta2=mu2;
    theta3=.5*sigma0^2*(-gamma);
    omega1=beta1-gamma;
    omega2=beta2-gamma;
    omega3=gamma-1;
    
    GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
    dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
        omega3*theta3*yy(wnStart:Nn).^(omega3-1);
    DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
    dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
        theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
        theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
    QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
    ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
    d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
    ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
    d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);
    
    VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
    d_dGG_VV_dyy(wnStart:Nn)= ...
        sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2) + ...
        sigma0.*omega2.*theta2.*(omega2+gamma-1).*yy(wnStart:Nn).^(omega2+gamma-2) + ...
        sigma0.*omega3.*theta3.*(omega3+gamma-1).*yy(wnStart:Nn).^(omega3+gamma-2);    

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this code block, we use second hermite representation of the density
%which is used to solve for mean reversion.
%It is the value of Sigma1 below whcih is used later in the main SDE one
%time step update formula.

    if(tt>1)
        [wMean,ZMean,Sigma00,Sigma11,Sigma1,Sigma2,Sigma0,w1New,w2New] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);

      %  tt
    else
        Sigma1(wnStart:Nn)=0.0;
        Sigma2(wnStart:Nn)=0.0;
        Sigma11=0.0;
    end
   
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %In this block we determine whether the SDE is mean reverting.
   %This is done, possibly naiively, by taking a derivative of drift at 
   %value and if the derivative is negative, we calculate KappaMR which is 
   %parameter that determines the strength of mean reversion. KappaMR is
   %used in the main SDE one step update formula.
   
    
    KappaMR(wnStart:Nn)=0;
    KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu1*yy(wnStart:Nn).^(beta1-1);
    KappaMR(wnStart:Nn)=KappaMR(wnStart:Nn)+mu2*yy(wnStart:Nn).^(beta2-1);
    KappaMR(wnStart:Nn)=abs(KappaMR(wnStart:Nn));
    %end
    if(tt<2)
        KappaMR(wnStart:Nn)=0.0;
        %SignKappaMR(wnStart:Nn)=1.0;
    end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the main SDE variable one step update equation. 


 w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
          sqrt(2*KappaMR(wnStart:Nn)).*(dt).*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+...
          sqrt(2/3).*(KappaMR(wnStart:Nn)).^1.0.*dt.^1.5.*Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));




%  w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
%           sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
%           sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
%           sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))- ...
%           1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*sigma0.^2.*dt.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%           %1*sqrt(2*KappaMR(wnStart:Nn)).*c1(wnStart:Nn).^2.*Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
    CTheta=1.0;
    D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
    D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
    wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
    wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
    wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
    wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
    %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%The last equation is not used since point w(NnMidh+3)is given but it is used in 
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.

    D1wh=(-11/6*w(NnMidh+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2) -26/3*w(NnMidh+3)+ 19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/dNn.^2;
    DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
%     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;


    whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
    whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
    whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;

    %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;


    w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
    w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
    w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
    w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
    w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

        Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
        Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
        Pn=1.0;%%Sum(ZProb(1:Nn))
        Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);
    
        Y0Pn=Y0*Pn;
        W0=1/(YnPn-Y0Pn);
        YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
        wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
    %    w(wnStart:Nn)=wCorrected(wnStart:Nn);
        
   %I have disabled the mean correction. The above mean correction is only valid for 
   %standard mean reverting stochastic volatility type SDEs. To enable mean
   %correction please uncomment the above last line in the block.
        
        
end

 y_w(1:Nn)=0;
 y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
 Dfy_w(wnStart:Nn)=0;
 Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
    fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
    fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc
str=input('Analytic Values calculated; Press a key to start monte carlo');

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.
 rng(85109634, 'twister')
 paths=200000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=500*gamma^2/.25*sigma0;
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 
end
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 9th, 2020, 10:59 am

Though we have reasonably good results with our exponential integral as calculated in post 788, I was thinking that I should try and see how the results are with Girsanov integral which is another exponential integral used to convert noise SDE into SDE with drift. Just in case. So I will be doing that experiment and reporting the results here.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 9th, 2020, 6:01 pm

For those friends who are familiar with the program very well, here is the update. I tried two different versions of KappaMR given here.

KappaMR=mu1*yy^(beta1-gamma)+mu2*yy^(beta2-gamma)-.5*gamma*sigma0^2*yy^(gamma-1)
and 
KappaMR=(1-gamma)*(mu1*yy^(beta1-gamma)+mu2*yy^(beta2-gamma)-.5*gamma*sigma0^2*yy^(gamma-1))

You can keep powers on KappaMR, dt and hermite polynomials exactly the same in the last two lines of the density update equation and play around with the coefficients, you will get very good results with one of the above. Especially close to zero, the results were thrilling around x0=.05 to x0=.15;
I wanted to complete and post the programs but my programs were corrupted so I just posted this update.
 
User avatar
Amin
Topic Author
Posts: 2499
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

February 11th, 2020, 11:52 am

For all the friends who read my threrad, I have been doing research independently on monte carlo simulation methods, Ito-Taylor  SDE expansions, and other interesting methods of the SDE densities evolution. Recently I also worked on a technique that is very promising for path integrals simulation and functionals of SDEs. And on the way I also proposed a method for the solution of the ODEs and systems of ODEs. 
I have had reasonable success in my research over past three or four years. But my knowledge about various branches of mathematics and even stochastics is rather limited. I was never trained as a mathematician and my highest degree is a masters in Economics though I learnt a lot of mathematics when I did my electrical engineering degree. I did continue to learn a reasonable lot of mathematics and stochastics when I was working for a Japanese derivatives modeling company. And then through my independent research. So I have a lot of background in stochastics modeling, and also I have done some decent work in filtering, parameter estimation and inference. So I have decent knowledge in stochastics and mathematics but there are really a lot of areas and advanced topics where my knowledge is totally limited.
I wanted to request some seasoned and expert professors in stochastics and mathematics if they would be willing to guide my research in stochastics, possibly discussing ideas and giving advice and suggesting research directions. If you share research interests, I hope that it could become a possibly productive effort to contribute to mathematics and stochastics. I look forward to devoting next few years doing research in stochastics and mathematics. So if you liked my past research efforts and would like to guide me in the future and share the desire to promote interesting research in mathematics and stochastics, I would request you to email me or mention here. If friends can mention/email me about some expert professors who could be willing to guide me and give suggestions, I would love to contact the respected professors on my own. There is absolutely no involvement of money here and simply just some joint and guided effort to do interesting research.
I do want to mention that I would want to continue to share my research with public as it evolves just like I have been doing previously. And I really want to remain based in Pakistan independently.
My email is anan2999(at)yahoo.com
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


Twitter LinkedIn Instagram

JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...


GZIP: On