Page 53 of 56

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

Posted: January 28th, 2019, 2:44 pm
I am testing two methods in this program.1. Integral of the squared vol method.
2. Squared orthogonal increments method.
Second method is far more accurate but suffers from instability in the mid of the SDE density. I have ideas how to fix it and will soon be uploading
a better version of this. I have earlier posted dependencies of this program but will post them again in a few hours.
function [] = ItoTaylorDensityWilmottSimplev12()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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 directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.
%The program has instability in 3-4 subdivisions in the mid of the density
%and I will be relasing a program in another few days taking care of that
%and fixing the problem. You can do it on your own as well.
% I am testing two methods.1. Integral of the squared vol method.
% 2. Squared orthogonal increments method.
%Second method is far more accurate but suffers from instability in the mid
%of the SDE density. I have ideas how to fix it and will soon be uploading
%a better version of this.
dt=.125/4;   % Simulation time interval.
Tt=32;     % Number of simulation levels. Terminal time= Tt*dt;
OrderA=4;  %expansion order for the analytic SDE.
OrderM=4;  %expansion order for the monte carlo simulation.
%Due to variable declaration constrint is that orderA>=OrderM.
%Current max expansion order is 4.
dNn=.1/1;   % Normal density width. would change with number of subdivisions
Nn=80;  % No of normal density subdivisions
%NnMidn=81;%Median cell/subdivision.
NnMidl=40;
NnMidh=41;

NnMid=4.0;
x0=1.0;   % starting value of SDE

beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.550;   % volatility power.
kappa=.80;   %mean reversion parameter.
theta=1.0;   %mean reversion target
%Below you can get rid of kappa and theta and freely choose your own values
%for both drift coefficients.
mu1=0*theta*kappa;   % first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.650;%Volatility value
alpha=1;% x^alpha is being expanded
alpha1=1-gamma;%For the expansion in transformed coordinates.

for nn=1:Nn
x(nn)=x0^(1-gamma)/(1-gamma); %Initial value is assigned to each
%subdivision in transformed bessel
%coordinates.
y(nn)=x0;
w(nn)=x0^(1-gamma)/(1-gamma);
end
%Please read notes at the end of the program for understanding of what is
%Y(l1,l2,l3,l4) and what are l1,l2, l3 and l4. In the ItoTaylor expansion
%l1 is first drift term expansion order index, l2 is for the second drift
%term, l3 is for quadratic variation drift term index and l4 is the
%volatility expansion order index.
%These sigma11, mu11, mu22, sigma22 can be lumped with Y(l1,l2,l3,l4) but I
%kept them separate for a general program where these values could change
%on every step while exponents on x in each term would remain same in
%pre-calcualted Y(l1,l2,l3,l4)
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:Tt+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) = dt^((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));
%Below power uses alpha1 associated with bessel
%coordinates. alpha1 is initial power being expanded. It will be further divided by (1-gamma)
%When used in the body of the program.
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

for nn=1:Nn
Z(nn)=((nn-.5)*dNn-NnMid);
end
%Z
%str=input('Show Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=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
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%The routine below calls the Ito-Taylor Coeffs in transformed coordinates
%and uses alpha1 as the power where alpha1=1-gamma.
Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)

for nn=1:Nn
Z(nn)=((nn-.5)*dNn-NnMid);%*Rho1(nn);
HermiteP(1,nn)=0;
HermiteP(2,nn)=Z(nn);%rho1(nn)
HermiteP(3,nn)=(Z(nn)^2-1);%rho2(nn)
HermiteP(4,nn)=(Z(nn)^3-3*Z(nn));
HermiteP(5,nn)=(Z(nn)^4-6*Z(nn)^2+3);

end
HermiteVar(1)=0;
HermiteVar(2)=1;%rho1(nn)
HermiteVar(3)=2;%rho2(nn)
HermiteVar(4)=6;
HermiteVar(5)=24;
HermiteVar(6)=120;
HermiteVar(7)=720;
HermiteVar(8)=720*7;
HermiteVar(9)=720*7*8;
nStart=1;
wnStart=1;

for tt=1:Tt
for nn=1:Nn
x1(nn)=x(nn);
if(tt==1)
Sigma2dt(nn)=0.0;
Sigma22dt(nn)=0.0;
Sigma2dtPrev(nn)=0.0;
end
Mu0dt(nn)=0.0;

wMu0dt(nn)=0.0;

Sigma2dtPrev(nn)=.5*Sigma2dt(nn)+.5*Sigma22dt(nn);
for k = 0 : OrderA
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<=1)
%Mu0dt is mean associated with integral of vol
%method.
Mu0dt(nn)=Mu0dt(nn) + YqCoeff(l1,l2,l3,l4) .* ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* Ft(l1,l2,l3,l4);

%wMu0dt is mean associated with squared orthogonal
%increments method.

wMu0dt(nn)=wMu0dt(nn) + YqCoeff(l1,l2,l3,l4) .* ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
((1-gamma)*w(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* Ft(l1,l2,l3,l4);

%((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) is
%term for the original untransformed SDE variable.
%This is basically ODE for each subdivision where
%initial value of ODE is reset at each time level
%after addition of diffusion increment fro previous
%time levelthat is different for each subdivision.
end

if(l4>=2)

%%Sigma2dt integrates over time. It is squared vol
%%associated with bessel rpocess and used for
%%integral of vol method.
Sigma2dt(nn)=Sigma2dt(nn) + YqCoeff(l1,l2,l3,l4).^2 .* ...
(mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)).^2.* ...
((1-gamma)*x(nn)).^(2*Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
HermiteVar(l4) .* Ft(l1,l2,l3,l4).^2;

end
end
end
end
end
dw(nn)=sigma0*sqrt(dt)*Z(nn);

%Below are intermediate calculations for integral of vol method.
if(tt==1)
x2(nn)=x0^(1-gamma)/(1-gamma)+ Mu0dt(nn) + sqrt(Sigma2dt(nn))*Z(nn);

end
if(tt>1)
x2(nn)=x1(nn)+ Mu0dt(nn) + (sqrt(Sigma2dt(nn))-sqrt(Sigma2dtPrev(nn)))*Z(nn);
end
%Below I tried some sort of predictor corrector but it does not seem to
%have much benefit.

Mu2dt(nn)=0.0;
%Sigma22dt(nn)=0.0;
for k = 0 : OrderA
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<=1)

Mu2dt(nn)=Mu2dt(nn) + YqCoeff(l1,l2,l3,l4) .* ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
((1-gamma)*x2(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* Ft(l1,l2,l3,l4);

end

if(l4>=2)

%%Sigma2dt integrates over time
Sigma22dt(nn)=Sigma22dt(nn) + YqCoeff(l1,l2,l3,l4).^2 .* ...
(mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)).^2.* ...
((1-gamma)*x2(nn)).^(2*Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
HermiteVar(l4) .* Ft(l1,l2,l3,l4).^2;

end
end
end
end
end
%Below I have advanced the value of bessel cooridinates variable in
%squared vol method.
if(tt==1)
x(nn)=x0^(1-gamma)/(1-gamma)+ .5*(Mu0dt(nn)+Mu2dt(nn)) + sqrt(Sigma2dt(nn))*Z(nn);
end
if(tt>1)
x(nn)=x1(nn)+.5* (Mu2dt(nn)+Mu0dt(nn)) + (sqrt(.5*Sigma2dt(nn)+.5*Sigma22dt(nn))-sqrt(Sigma2dtPrev(nn)))*Z(nn);
end

end
EnterLoop1=1;
EnterLoop2=1;
for nn=Nn-1:-1:nStart
%First loop makes sure that values are monotonic and the values
%very close to zero that become unbounded are turned to zero.
if((x(nn)>x(nn+1))&&(EnterLoop1==1))
x(1:nn)=0;
nStart=nn+1;
EnterLoop1=0;
end
%Second loop makes sure that values smaller than zero that would
%become complex are turned to zero.
if((x(nn)<0)&&(EnterLoop2==1))
x(1:nn)=0.0;
nStart=nn+1;
EnterLoop2=0;
end
end
%Now calculations start for squared orthogonal increments method.
if(tt>0)
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.

for nn=wnStart:Nn
if(w(nn)>=wMeanPrev)
sign_w(nn)=1;
else
sign_w(nn)=-1;
end

if(dw(nn)>=0.0)
sign_dw(nn)=1;
else
sign_dw(nn)=-1;
end
if((w(nn)-wMeanPrev)+dw(nn)>=0)
sign_Sum(nn)=+1;
else
sign_Sum(nn)=-1;
end

end

for nn=wnStart:Nn
%This is the squared addition for t>0 time steps.
w(nn)=wMeanPrev+wMu0dt(nn)+sign_Sum(nn).*sqrt(abs(sign_w(nn)*  ...
(w(nn)-wMeanPrev).^2+sign_dw(nn).*dw(nn).^2));

end

EnterLoop1=1;
EnterLoop2=1;
for nn=Nn-1:-1:wnStart
%Makes sure values decrease monotonically close to zero.
if((w(nn)>w(nn+1))&&(EnterLoop1==1))
w(1:nn)=0;
wnStart=nn+1;
EnterLoop1=0;
end
%Makes sure values are not smaller than zero.
if((w(nn)<0)&&(EnterLoop2==1))
w(1:nn)=0.0;
wnStart=nn+1;
EnterLoop2=0;
end
end

end

end
%x(x<0)=0.0;
for nn=nStart+1:Nn
y_x(nn) = ((1-gamma)*x(nn)).^(1/(1-gamma)); %Convert squared vol values
%from bessel to original coordinates
y_w(nn) = ((1-gamma)*w(nn)).^(1/(1-gamma)); %Convert squared orthogonal
%increment method values from bessel to original coordinates
end
for nn=2:Nn-1
Dfx(nn) = (x(nn + 1) - x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));

Dfy_x(nn) = (y_x(nn + 1) - y_x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end

fw(1:Nn)=0;
fx(1:Nn)=0;
fy_x(1:Nn)=0;
for nn = 2:Nn-1
fx(nn) = normpdf(Z(nn),0, 1)/abs(Dfx(nn)); %Bessel process density
fw(nn) = normpdf(Z(nn),0, 1)/abs(Dfw(nn)); %Bessel process density

fy_x(nn) = normpdf(Z(nn),0, 1)/abs(Dfy_x(nn));%Origianl coordinates density
fy_w(nn) = normpdf(Z(nn),0, 1)/abs(Dfy_w(nn));%Origianl coordinates density
end
str=input('Analytic Values Calculated. Press a key to start Monte Carlo');
%Calculate the density with monte carlo below.
rng(85109634, 'twister')
paths=100000;
X(1:paths)=x0^(1-gamma)/(1-gamma); %Bessel process Monte carlo
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:Tt
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;
%HermiteP1(6,1:paths)=Random1(1:paths).^5-10*Random1(1:paths).^3+15*Random1(1:paths);
%HermiteP1(7,1:paths)=Random1(1:paths).^6-15*Random1(1:paths).^4+45*Random1(1:paths).^2-15;
XX0=X;
YY0=YY;
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);

%Above is bessel process monte carlo equation.

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

X(X<0)=0;
sum(X(:))/paths   %monte carlo average for bessel coordinates
sum(x(nStart:Nn-1).*ZProb(nStart:Nn-1))  % Ito-Taylor density average for bessel coordinates
sum(w(wnStart:Nn-1).*ZProb(wnStart:Nn-1))% Ito-Taylor density average for
%with squared orthogoanl increments in bessel coordinates
YY(YY<0)=0;
sum(YY(:))/paths %origianl coordinates monte carlo average.
sum(y_x(nStart:Nn-1).*ZProb(nStart:Nn-1)) %Original process average from coordinates
sum(y_w(wnStart:Nn-1).*ZProb(wnStart:Nn-1))   %derived from bessel process
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=30;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff);
plot(w(wnStart:Nn-1),fw(wnStart:Nn-1),'b',x(nStart:Nn-1),fx(nStart:Nn-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
str=input('Look at Bessel graphs comparison between monte carlo density and analytic density.Blue graph is from squared orthogonal incrments method and red is from squared volatility, green is monte carlo.');

[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b',y_x(nStart+1:Nn-1),fy_x(nStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('Look at Original variable graphs comparison between monte carlo and analytic density');
str=input('Blue graph is from squared orthogonal incrments method and red is from squared volatility, green is monte carlo.');
end

% In the expansion of X^alpha
% With SDE
% dX(t)=mu1 * X(t)^beta1 dt + mu2* X(t)^beta2 dt + sigma * X(t)^gamma dz(t)
% The first order expansion of X^alpha with the above SDE is
% X(t+1)^alpha= X(t)^alpha +alpha* X(t)^(alpha-1) * mu1 * X(t)^beta1 *dt
% +alpha* X(t)^(alpha-1) * mu2 * X(t)^beta2 *dt
% + .5* alpha *(alpha-1)* X(t)^(alpha-2) * sigma^2 * X(t)^(2*gamma) *dt
% + alpha * X(t)^(alpha-1) * sigma * X(t)^gamma *dz(t)
% In my program, for the above four term expansion, I use a four dimensional
%array notation that shows the order of differentiation on each of these
%four types of terms. So if l1 represents drift one, l2 represents drift
%two and l3 represents quadratic variation term and l4 represents
%volatility term, we can use the notation Y(l1,l2,l3,l4) for the
%coefficients of expansion of the above term. So in the above first order expansion
% Y(1,0,0,0)= alpha * mu1
% Y(0,1,0,0)=alpha *mu2
% Y(0,0,1,0)= .5* alpha *(alpha-1)* sigma^2
% Y(0,0,0,1)= alpha * sigma
%
% For higher orders of expansion, we will have to carefully compute the
%coefficients separately based on path dependence and then lump them
%together in Y according to array indices. Once we have done that,
%we can write one term of the general expansion as
% Y(l1,l2,l3,l4) * x^(alpha+l1 * beta1 +l2* beta2 +l3* 2* gamma +l4 * gamma
%  - l1 -l2 - 2*l3 -l4) * HermiteP(l4) * dt^(l1 + l2 +l3 + .5*l4);
%
% In my program, since matlab does not take indices starting from zero,
%I have slightly adjusted the above notation.
%In the density simulation program, we pre-compute all coefficients and values
%needed in simulation but we have to calculate the part of above
%expression where power of x is evaluated.
% If we have four terms in Ito-hermite SDE expansion at first order,
%In every next level, each of the four terms would take four descendent
%terms and there will be sixteen terms on second level, and on third level
%there will be 64 descendent terms and 256 descendent terms on fourth level.
%However we can reduce the number of terms in Ito-hermite expansion to
%number of terms in polynomial expansion. For example Ito-hermite expansion
%of SDE with four terms on first order can be reduced to the number of
%terms in expansion of a polynomial with four terms on each level.
%To do this, we will have to pre-calculate the Coefficients carefully and
%the lump them together as in a polynomial like expansion.
% Again the above array Y(l1,l2,l3,l4) is relatively sparse and
%does not take all possible indices. Just like in polynomial expansion,
%on a certain expansion order, it takes only indices so that sum of indices
%goes to expansion order. For example on second order, only those indices
%could be found such that l1+l2+l3+l4=2. And similarly for third ito-hermite
%expansion level only those indices would be populated such that l1+l2+l3+l4=3.
% In order to do that I have done simple loops in a fashion shown below
%
% for k = 0 : (Order4)
%     for m = 0:k
%         l4 = k - m;
%         for n = 0 : m
%             l3 = m - n;
%             for j = 0:n
%                 l2 = n - j;
%                 l1 = j;
% First loop variable is for expansion order. And then rest of the loop
%variables take values such that their sum equals the expansion order.
%When we would add a fifth term in the SDE expansion for example,
%we will have to add another loop variable with the same logic that
%sum of loop variables equal to expansion order. In my matlab program,
%I have given offset to l4, l2,l3, l1 since array indices cannot start at zero.

%
%
% % In the expansion of X^alpha
% % With SDE
% % dX(t)=mu1 * X(t)^beta1 dt + mu2* X(t)^beta2 dt + sigma * X(t)^gamma dz(t)
% % The first order expansion of X^alpha with the above SDE is
% % X(t+1)^alpha= X(t)^alpha +alpha* X(t)^(alpha-1) * mu1 * X(t)^beta1 *dt
% % +alpha* X(t)^(alpha-1) * mu2 * X(t)^beta2 *dt
% % + .5* alpha *(alpha-1)* X(t)^(alpha-2) * sigma^2 * X(t)^(2*gamma) *dt
% % + alpha * X(t)^(alpha-1) * sigma * X(t)^gamma *dz(t)
% % In my program, for the above four term expansion, I use a four dimensional
% %array notation that shows the order of differentiation on each of these
% %four types of terms. So if l1 represents drift one, l2 represents drift
% %two and l3 represents quadratic variation term and l4 represents
% %volatility term, we can use the notation Y(l1,l2,l3,l4) for the
% %coefficients of expansion of the above term. So in the above first order expansion
% % Y(1,0,0,0)= alpha * mu1
% % Y(0,1,0,0)=alpha *mu2
% % Y(0,0,1,0)= .5* alpha *(alpha-1)* sigma^2
% % Y(0,0,0,1)= alpha * sigma
% %
% % For higher orders of expansion, we will have to carefully compute the
% %coefficients separately based on path dependence and then lump them
% %together in Y according to array indices. Once we have done that,
% %we can write one term of the general expansion as
% % Y(l1,l2,l3,l4) * x^(alpha+l1 * beta1 +l2* beta2 +l3* 2* gamma +l4 * gamma
% %  - l1 -l2 - 2*l3 -l4) * HermiteP(l4) * dt^(l1 + l2 +l3 + .5*l4);
% %
% % In my program, since matlab does not take indices starting from zero,
% %I have slightly adjusted the above notation.
% %In the density simulation program, we pre-compute all coefficients and values
% %needed in simulation but we have to calculate the part of above
% %expression where power of x is evaluated.
% % If we have four terms in Ito-hermite SDE expansion at first order,
% %In every next level, each of the four terms would take four descendent
% %terms and there will be sixteen terms on second level, and on third level
% %there will be 64 descendent terms and 256 descendent terms on fourth level.
% %However we can reduce the number of terms in Ito-hermite expansion to
% %number of terms in polynomial expansion. For example Ito-hermite expansion
% %of SDE with four terms on first order can be reduced to the number of
% %terms in expansion of a polynomial with four terms on each level.
% %To do this, we will have to pre-calculate the Coefficients carefully and
% %the lump them together as in a polynomial like expansion.
% % Again the above array Y(l1,l2,l3,l4) is relatively sparse and
% %does not take all possible indices. Just like in polynomial expansion,
% %on a certain expansion order, it takes only indices so that sum of indices
% %goes to expansion order. For example on second order, only those indices
% %could be found such that l1+l2+l3+l4=2. And similarly for third ito-hermite
% %expansion level only those indices would be populated such that l1+l2+l3+l4=3.
% % In order to do that I have done simple loops in a fashion shown below
% %
% % for k = 0 : (Order4)
% %     for m = 0:k
% %         l4 = k - m;
% %         for n = 0 : m
% %             l3 = m - n;
% %             for j = 0:n
% %                 l2 = n - j;
% %                 l1 = j;
% % First loop variable is for expansion order. And then rest of the loop
% %variables take values such that their sum equals the expansion order.
% %When we would add a fifth term in the SDE expansion for example,
% %we will have to add another loop variable with the same logic that
% %sum of loop variables equal to expansion order. In my matlab program,
% %I have given offset to l4, l2,l3, l1 since array indices cannot start at zero. 

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

Posted: January 28th, 2019, 2:49 pm
Here are the dependencies of the above program.
function [Y] = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma)

%In the coefficient calculation program which calculates Y(l1,l2,l3,l4),
%I have used four levels of looping each for relevant expansion order.
%The first loop takes four values and second loop takes 16 values and
%third loop takes 64 values and so on. And then each coefficient
%term can be individually calculated while carefully accounting
%for path dependence.
%So for example in a nested loop structure
%m1= 1:mDim
% 	m2=1:mDim
%		m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;

%in the above looping loop takes values from one to four with one
%indicating the first drift term, two indicating the second drift term
%and three indicating quadratic variation term and
%four indicating the volatility term. And with this looping structure
%we can So in the above looping m1=1 would mean that all terms are
%descendent of first drift term and m2=4 would mean that all terms are
%descendent of first drift term on first expansion order and descendent
%of volatility term on the second order and so we can keep track of path
%dependence perfectly.
%And on each level, we individually calculate the descendent terms. While
%keeping track of path dependence and calculating the coefficients with
%careful path dependence consideration, we update the appropriate element
%in our polynomial like expansion coefficient array

%explaining the part of code
%m1= 1:mDim
%	m2=1:mDim
%		m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;
%Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

%Here l(1) denotes l1 but written as l(1) so it can be conveniently
%updated with the loop variable when the loop variable takes value one
%indicating first drift term . And l(2) could be conveniently updated when
%the loop variable takes value equal to two indicating second
%drift term and so on.
%Here is the part of code snippet for that

%for m1=1:mDim
%    l(1)=1;
%    l(2)=1;
%    l(3)=1;
%    l(4)=1;
%    l(m1)=l(m1)+1;
%CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
%    CoeffDX2 = CoeffDX1 - 1;
%    ArrIndex0=m1;
%    ArrIndex=(m1-1)*mDim;
%    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%    Y2(1+ArrIndex)=Coeff1st;
%    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(2+ArrIndex)=Coeff1st;
%    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(3+ArrIndex)=Coeff2nd;
%    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
%    Y2(4+ArrIndex)=Coeff1st;
%    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));

%The first four lines update the array indices according to the parent term.
%And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.
%ArrIndex0=m1; calculates the array index of the parent term
%And ArrIndex=(m1-1)*mDim; calculates the array index of the descendent terms
%And coefficient of the drift and volatility descendent terms is
%calculated by multiplying the coefficient of the parent term by
%Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%And coefficient of the quadratic variation descendent terms is
%calculated by multiplying the coefficient of the parent term by
%Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%And then each of the four descendent terms are updated with Coeff1st
%if they are drift or volatility descendent terms or Coeff2nd if
%they are quadratic variation descendent terms.
%Here Y1 indicates the temporary coefficient array with parent terms on
%first level. And Y2 denotes temporary coefficient array with parent terms
%on second level and Y3 indicates temporary coefficient array with parent terms on third level.

[IntegralCoeff,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs();

n1(1)=2;
n1(2)=2;
n1(3)=2;
n1(4)=3;
%n1(1), n1(2), n1(3) are drift and quadratic variation term variables
%and take a value equal to 2 indicating a dt integral.
%n1(4) is volatility term variable and indicates a dz-integral by taking a
%value of 3.

mDim=4; % four descendent terms in each expansion
Y(1:5,1:5,1:5,1:5)=0;
Y1(1:mDim)=0;
Y2(1:mDim*mDim)=0;
Y3(1:mDim*mDim*mDim)=0;

%First Ito-hermite expansion level starts here. No loops but four
%descendent terms.
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
CoeffDX1 = alpha;
CoeffDX2 = CoeffDX1 - 1;
Coeff1st=CoeffDX1;
Coeff2nd=.5*CoeffDX1*CoeffDX2;

Y1(1)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(2)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(3)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,2);
Y1(4)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,3);

%Second Ito-hermite expansion level starts. It has a loop over four parent
%terms and there are four descendent terms for each parent term.
%The coefficient terms are then lumped in a polynomial-like expansion
%array of coefficients.
for m1=1:mDim
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
l(m1)=l(m1)+1;
CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
- (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
CoeffDX2 = CoeffDX1 - 1;
ArrIndex0=m1;
ArrIndex=(m1-1)*mDim;
Coeff1st=Y1(ArrIndex0)*CoeffDX1;
Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
Y2(1+ArrIndex)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
Y2(2+ArrIndex)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
Y2(3+ArrIndex)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
Y2(4+ArrIndex)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));
%Third Ito-hermite expansion level starts and it is a nested loop with
%a total of sixteen parents and each parent takes four descendent
%terms.
for m2=1:mDim
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
l(m1)=l(m1)+1;
l(m2)=l(m2)+1;
CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
- (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
CoeffDX2=CoeffDX1-1;
ArrIndex0=(m1-1)*mDim+m2;
ArrIndex=((m1-1)*mDim+(m2-1))*mDim;
Coeff1st=Y2(ArrIndex0)*CoeffDX1;
Coeff2nd=Y2(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
Y3(1+ArrIndex)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
Y3(2+ArrIndex)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
Y3(3+ArrIndex)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,2,n1(m2),n1(m1));
Y3(4+ArrIndex)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,3,n1(m2),n1(m1));
%fourht Ito-hermite expansion level starts and it is a triply-nested loop with
%a total of sixteen parents and each parent takes four descendent
%terms. We then lump the terms in a relatively sparse polynomial
%like expansion coefficient array that has smaller number of
%non-zero terms.

for m3=1:mDim
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
l(m1)=l(m1)+1;
l(m2)=l(m2)+1;
l(m3)=l(m3)+1;
CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
- (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
CoeffDX2=CoeffDX1-1;
ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
%ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
Coeff1st=Y3(ArrIndex0)*CoeffDX1;
Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%Y4(1+ArrIndex)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(2+ArrIndex)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(3+ArrIndex)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(4+ArrIndex)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m3),n1(m2),n1(m1));
end
end
end

end


and
function [IntegralCoeff0,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs()

IntegralCoeff(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdt(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdz(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1,1,1,2)=1;
IntegralCoeff0(1,1,1,3)=1;

IntegralCoeff0(1,1,2,2)=1/2;
IntegralCoeff0(1,1,3,2)=1-1/sqrt(3);
IntegralCoeff0(1,1,2,3)=1/sqrt(3);
IntegralCoeff0(1,1,3,3)=1/2;

IntegralCoeff0(1,2,2,2)=1/6;
IntegralCoeff0(1,3,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5));
IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(1-1/sqrt(5));
IntegralCoeff0(1,3,3,2)=1/2*(1-sqrt(2)/2);
IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5);
IntegralCoeff0(1,3,2,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2;
IntegralCoeff0(1,2,3,3)=1/sqrt(3)/sqrt(2)*1/(2);
IntegralCoeff0(1,3,3,3)=1/6;

IntegralCoeff0(2,2,2,2)=1/24;
IntegralCoeff0(2,2,3,2)=1/2*1/sqrt(5)*1/3*(1-1/sqrt(7));
IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,2,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,3,2,2)=1/2*(1-sqrt(2)/2)*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,2,3,2)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,3,3,2)=1/6*(1-sqrt(3)/sqrt(5));

IntegralCoeff0(2,2,2,3)=1/6*1/sqrt(7);
IntegralCoeff0(2,2,3,3)=1/2*1/sqrt(5)*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(2,3,2,3)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,2,2,3)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,3,2,3)=1/2*(1-sqrt(2)/2)*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(2,3,3,3)=1/sqrt(3)*1/sqrt(2)*1/sqrt(4)*1/sqrt(3)*(1/sqrt(5));
IntegralCoeff0(3,2,3,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(3,3,3,3)=1/24;

%Can also be calculated with algorithm below.
%here IntegralCoeffdt indicates the coefficients of a dt-integral.
%This dt-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dt with X(t) dynamics given by the SDE
%here IntegralCoeffdz indicates the coefficients of a dz-integral.
%This dz-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dz(t) with X(t) dynamics given by the SDE

%IntegralCoeff is is associated with expansion of X(t)^alpha.
%IntegralCoeff below is not returned and IntegralCoeff0 manually calculated
%above is returned but both are the same.

l0(1:2)=1;

for m4=1:2
l0(1)=1;
l0(2)=1;

%IntegralCoeff4(m4,1,1,1)=1;
%IntegralCoeff4(m4,1,1,1)=1;
%1 is added to m4 since index 1 stands for zero, 2 for one and three
%for two.
IntegralCoeff(1,1,1,m4+1)=1;
l0(m4)=l0(m4)+1;
IntegralCoeffdt(1,1,1,m4+1)=IntegralCoeff(1,1,1,m4+1)* ...
1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
IntegralCoeffdz(1,1,1,m4+1)= IntegralCoeff(1,1,1,m4+1)* ...
1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
for m3=1:2
l0(1)=1;
l0(2)=1;
l0(m4)=l0(m4)+1;

if(m3==1)
IntegralCoeff(1,1,m4+1,m3+1)=1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
end
if(m3==2)
IntegralCoeff(1,1,m4+1,m3+1)= 1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
l0(m3)=l0(m3)+1;
%IntegralCoeff(1,1,m4+1,m3+1)=IntegralCoeff4(m4,m3,1,1);
IntegralCoeffdt(1,1,m4+1,m3+1)=IntegralCoeff(1,1,m4+1,m3+1)* ...
1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
IntegralCoeffdz(1,1,m4+1,m3+1)= IntegralCoeff(1,1,m4+1,m3+1)* ...
1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
for m2=1:2
l0(1)=1;
l0(2)=1;
l0(m4)=l0(m4)+1;
l0(m3)=l0(m3)+1;

if(m2==1)
IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m4+1,m3+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
end
if(m2==2)
IntegralCoeff(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,1,m4+1,m3+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
l0(m2)=l0(m2)+1;
%IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff4(m4,m3,m2,1);
IntegralCoeffdt(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
IntegralCoeffdz(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);

for m1=1:2
l0(1)=1;
l0(2)=1;
l0(m4)=l0(m4)+1;
l0(m3)=l0(m3)+1;
l0(m2)=l0(m2)+1;
if(m1==1)
IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
end
if(m1==2)
IntegralCoeff(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
l0(m1)=l0(m1)+1;
%IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff4(m4,m3,m2,m1);
IntegralCoeffdt(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
IntegralCoeffdz(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
end
end
end

end
And the graphing of simulation file is here.
function [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti(X,Paths,BinSize,MaxCutOff )
%Processes monte carlo paths to return a series Xdensity as a function of IndexOut. IndexMax is the maximum value of index.
%

Xmin=0;
Xmax=0;
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
%if(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=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



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

Posted: January 28th, 2019, 3:39 pm
I posted this on internet about my latest research. I will continue to post my research and programs on parameter estimation of SDEs and state-space processes and related filtering on Wilmott forum. And then further continue with  posting research on AI + Bayesian methods on this forum. Here is my post on linkedin about my research.

Introduction to Newly Developed Research at Infiniti Derivatives Technologies.

Fast and Accurate Monte Carlo Simulation.

Due to recently developed research at Infiniti Derivatives Technologies, we can do accurate monte carlo simulations of almost any stochastic differential equation. This includes arbitrary stochastic volatility SDEs and their related greeks SDEs. We challenge you to throw examples of difficult SDEs that we would solve with our newly developed technology.

Discovery of a new SDE Density Simulation Framework.

We have developed a new SDE density simulation framework that is very fast and reasonably robust. Here are some examples of the applications of this technology
• Fast Calibration of derivatives models and calculation of their greeks.
• New innovative models. For example you could now have a SV SABR model in which volatility is governed by a CEV process of arbitrary quotient. We could for example help you simulate, calibrate and risk manage a SV Libor Market Model framework in which SABR vol is given by a CEV process.
• Formation of alpha generation strategies with options. Our density simulation framework can allow you to calculate the right distribution of the underlying asset and you could then trade options according to mispricing with respect to implied vol in the market.
• Modeling of mean-reverting densities for statistical arbitrage in pairs trading.
•  The method is ideal for calibration of volatility surfaces and can very easily be enhanced to calibrate and simulate the local volatility and LSV models.
""
""
Filtering and Parameter Estimation:

Our research on SDEs has ramifications for non- linear Bayesian filtering and parameters estimation in the area of SDEs and state-space models. Using this new research, we can do on-line parameter optimization and tuning for the parameters of the SDE or the state space process. It can also be possible to use slowly varying parameters framework for filtering of the SDE or the state space process. I have also done research work that obviates the need for expensive integrations in non-linear Bayesian filtering and the whole work can be done in transition probabilities framework. This makes the process both faster and more accurate. With this new research, we can do non-linear filtering for non-linear state-space models where efficient filtering was earlier not possible. These are interesting developments in filtering and parameter estimation using my research on stochastic processes. This can have huge applications within and out of the context of quantitative finance. Here are some applications in quantitative finance:
• We can find the true distribution of financial data generation process using filtering and parameter estimation techniques. We can also find out how volatility mean reverts after using filtering and calculating the data distribution of data generating process on several scales. This will give us true data generating process and tell us about the true mean reversion in the data as opposed to blindly following the implied distributions backed out from the option pricing data. We can for example find out what is the true data generating process for sixty dimensional LIBOR rate process used in LIBOR market model and what is the structure of their daily volatility and how this structure scales to weekly, monthly and quarterly volatility and how future volatility changes. This can be a very exciting application of our parameter estimation algorithm and will help outline the true LIBOR volatilities which can then be compared to LIBORs implied volatility for statistical arbitrage purposes.
• filtering of financial data can be used in the context of quantitative trading.
• We can use filtering in the context of pairs trading where the mean reverting parameter can be calculated online and can also be modeled as a slowly varying parameter. Its intensity could tell about the mean reversion speed of the pair in pairs trading.
I wish to work with your company as a consulting partner to implement some of these new technologies. If you have any questions or you would like me to work for you on these projects for your company, please email me at anan2999@yahoo.com

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

Posted: January 28th, 2019, 8:32 pm
I will continue to post my research and programs on parameter estimation of SDEs and state-space processes and related filtering on Wilmott forum. And then further continue with  posting research on AI + Bayesian methods on this forum.
thanks
I for one always look forward to your posts

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

Posted: January 29th, 2019, 2:38 pm
Thank you Peter for your kind words. Yes, I hope to continue putting interesting stuff on Wilmott here. I wish everybody likes this knowledge sharing. Thank you again.

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

Posted: February 7th, 2019, 9:09 pm
OK friends, I am excited as there is another major breakthrough in the Ito-Taylor density simulation method. I had done a paper last year in which I had claimed that transition probabilities between different CDF subdivisions during the evolution of an SDE can be mapped and are the same as the transition probabilities between different CDF subdivisions of a brownian motion with the same volatility. Using the ideas from that paper, I have worked out a new method in which works well even in original coordinates and can also be used for evolution of dt-integrals and dz-integrals of the SDE. This is a major breakthrough. We noticed that in Ito Taylor method density is mapped on a normal density at each step. I mapped the simulation density before simulation step on a corresponding Brownian motion and then simulated the brownian motion with the same volatility as the SDE and found out what is the transition value of gaussian Z that is required to move from nth CDF interval at time t-1 to nth CDF interval at time t on the brownian motion density. I used the same transition value of gaussian Z in hermite polynomials for the simulation of  SDE and it can be done in the original coordinates as well. I am uploading a very simple program that shows how it is done. I have just simulated CEV noise using the density but it can also be very easily used with SDEs with linear and non-linear drifts. Here is the simple program. I will come up with a very detailed program very soon. Here is the simple program that simulates an SDE by mapping it on Brownian motion and using transition Z from the brownian motion transition values.
function [] = ItoTaylorDensityRecombining()
%In this quick program, I am simulating the driftless noise SDE
%dX(t)=sigma * X(t)^gamma * dz(t)
dt=.125/4;   % Simulation time interval.
Tt=16;       % Number of simulation levels. Terminal time= Tt*dt;
dNn=.1/1;    % Normal density width. would change with number of subdivisions
Nn=80;       % No of normal density subdivisions
NnMidl=40;
NnMidh=41;
NnMid=4.0;
x0=1.0;      % starting value of SDE
gamma=.50;   %volatility power exponent
sigma0=1.0;   %Volatility value

y(1:Nn)=x0;  %SDE variable
for nn=1:Nn
Z(nn)=((nn-.5)*dNn-NnMid);
end
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=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);
end

for tt=1:Tt
for nn=1:Nn
x1(nn)=x0+sigma0*sqrt((tt-1)*dt)*Z(nn);%Simulate brownian motion
%to the last date
end
for nn=1:Nn
x2(nn)=x0+sigma0*sqrt(tt*dt)*Z(nn);%Simulate brownian motion to
%the next date
end
for nn=1:Nn
Zt(nn,nn)=(x2(nn)-x1(nn))/sigma0/sqrt(dt);%Get the transition value
%of Z that is required to go from nth CDF subdivision of brownain
%motion at time t-1 to nth CDF subdivision of brownian motion
%at time t.
%use the same transition value of Z to generate hermite
%polynomials.
HermiteP(1,nn)=0;
HermiteP(2,nn)=Zt(nn,nn);%rho1(nn)
HermiteP(3,nn)=(Zt(nn,nn)^2-1);%rho2(nn)
HermiteP(4,nn)=(Zt(nn,nn)^3-3*Zt(nn,nn));
HermiteP(5,nn)=(Zt(nn,nn)^4-6*Zt(nn,nn)^2+3);
end
for nn=1:Nn
y(nn)=y(nn)+sigma0*sqrt(dt)*y(nn)^gamma.*HermiteP(2,nn) ...
+sigma0^2*gamma*y(nn)^(2*gamma-1)*dt.*HermiteP(3,nn)/2 ...
+sigma0^3*gamma*(2*gamma-1)*y(nn).^(3*gamma-2)*dt^1.5.*(HermiteP(4,nn))/6 ...
+sigma0^4*gamma*(2*gamma-1)*(3*gamma-2).*y(nn).^(4*gamma-3)*dt^2.0.*(HermiteP(5,nn))/24;
end
end

for nn=2:Nn-1
Dfy(nn) = (y(nn + 1) - y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
end
fy(1:Nn)=0;
for nn = 2:Nn-1
fy(nn) = normpdf(Z(nn),0, 1)/abs(Dfy(nn));%Origianl coordinates density
end

str=input('Analytic Values Calculated. Press a key to start Monte Carlo');
rng(85109634, 'twister')
paths=100000;
Y(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:Tt
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;
Y0=Y;

Y(1:paths)=Y(1:paths)+sigma0*sqrt(dt).*Y0(1:paths).^gamma.*HermiteP1(2,1:paths) ...
+sigma0^2*gamma*Y0(1:paths).^(2*gamma-1).*dt.*HermiteP1(3,1:paths)/2.0 ...
+sigma0^3*gamma*(2*gamma-1)*Y0(1:paths).^(3*gamma-2).*dt^1.5.*HermiteP1(4,1:paths)/6.0 ...
+sigma0^4*gamma*(2*gamma-1)*(3*gamma-2).*Y0(1:paths).^(4*gamma-3)*dt^2.0.*HermiteP1(5,1:paths)/24;
end
sum(y(2:Nn-1).*ZProb(2:Nn-1))
Y(Y<0)=0;
sum(Y(:))/paths
BinSize=.0275;
MaxCutOff=30;
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(Y,paths,BinSize,MaxCutOff);
plot(y(2:Nn-1),fy(2:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

end 
Though the above program works with driftless CEV SDEs, I will be uploading a detailed program in a day or two that would simulate every SDE and its dt-integrals and its dz-integrals. This has become possible with the new discovery.

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

Posted: February 22nd, 2019, 6:36 pm
Sorry that I have been away working on my market trading project. But for next few days, I will work on my Ito-Taylor SDE simulation methods research. Here is the agenda for next few days. OK, one more thing that we have been working with Ito-Taylor density simulation part and we have neglected the Ito-Taylor monte carlo simulation. In the next post, I will be putting here a program that will take an arbitrary dimensional correlated system of SDEs of the kinds that appear in mathematical finance and we will be able to do monte carlo simulation up to 2nd or 4th order accuracy as suggested by the user. Then I will upload another example code that will simulate correlated three factor SV LIBOR market model SDEs to 2nd or 4th order accuracy as suggested by the user. So there are going to be interesting examples. On the Ito-Taylor density simulation side, I have also almost completed work on squared orthogonal increments density simulation scheme to take care of the instability and I will also be posting that here before I go underground again to start working on market trading. Hope to have an exciting week ahead.

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

Posted: July 2nd, 2019, 4:26 pm
Very sorry for remaining away for such a long time even before my research was taken to conclusion. I have been working with market trading based on analytics with limit order book dynamics but two weeks ago I returned to the problem of simulation of stochastic differential equations and their stochastic integrals. I would like to take the research to conclusion since I spent more than two years working on the problem. And now I believe that I am extremely close to the solution of the problem and expect that I would be able to present a completely worked out matlab program within a week to ten days. OK during my research I came across some interesting equations that I would like to share since there is a possibility they could be solved by analytical means by experts. One benchmark of my research has been stochastic volatility type equations take two non-linear terms in their drift. These represent enough generality and are also difficult enough. Ok, let us consider the following SDE.

[$]dX_s= \mu_1 X_s^{\beta_1} ds + \mu_2 X_s^{\beta_2} ds + \sigma X_s^{\gamma} dz(s)[$]   Eq(1)

I have been working with bessel coordinates where the SDE can be transformed into bessel variable [$]Y_s=\frac{{X_s}^{(1-\gamma)}}{(1-\gamma)}[$]  as (after the application of Ito formula)

[$]d[\frac{{X_s}^{(1-\gamma)}}{(1-\gamma)}] = \mu_1 X_s^{(\beta_1-\gamma)} ds + \mu_2 X_s^{(\beta_2-\gamma)} ds + .5 \sigma^2 (-\gamma) X_s^{(\gamma-1)} ds + \sigma dz(s)[$]  Eq(2)

I had been trying to use Feynman exponential on this SDE. When we write the Fokker planck of the SDE, there is killing term in the FP PDE that correponds to Feynman exponential. Since there are no random numbers in my deterministic density evolution, it has to learn from the Fokker planck in addition to following the SDE. Drift and diffusion terms would remain the same but Feynman exponential slightly shifts the underlying distribution and the density of the mean reverting SDEs (at 100% vols) start to go off after six or nine months. I believe that this is due to not accounting the Feynman exponential. While I was working on this, I found something interesting.
If we pick an exponential intelligently and multiply it with Bessel form of the SDE and apply Ito formula, we can cancel the drift terms in the Bessel form with the exponential multiplying the Gaussian directly. I could easily solve such a form in my framework but there could possibly be analytic solutions of such forms as well. Here is my choice of exponential and I apply it to bessel transformation and use Ito here
[$]d[\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])[$]
[$] \frac{X_s^{(1-\gamma)}}{(1-\gamma)}][$]  Eq(3)
In the above stochastic product to be solved by Ito, the terms in the exponential have been so chosen that after the application of Ito dt-terms in the first part of Ito expansion that appear from the differential of the integral multiplied by bessel variable would cancel the drift terms in the second part from differential of the bessel variable. These drift terms in second part are from the bessel SDE form given in equation two. The above differential is expanded by Ito as
[$]d[\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv]) \frac{X_s^{(1-\gamma)}}{(1-\gamma)}][$]
[$]=\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])[$]
[$]* (-\mu_1(1-\gamma)X_s^{\beta_1-1}ds-\mu_2(1-\gamma)X_s^{\beta_2-1}ds-.5\sigma^2(-\gamma)(1-\gamma)X_s^{2\gamma-2}ds) * \frac{X_s^{(1-\gamma)}}{(1-\gamma)}[$]
[$]+\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])[$]
[$]*( \mu_1 X_s^{(\beta_1-\gamma)} ds + \mu_2 X_s^{(\beta_2-\gamma)} ds + .5 \sigma^2 (-\gamma) X_s^{(\gamma-1)} ds + \sigma dz(s))[$]

After algebra, the above simplifies to
[$]d[\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv]) \frac{X_s^{(1-\gamma)}}{(1-\gamma)}][$]
[$]=\exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])\sigma dz(s))[$]
After taking integrals on both sides, we get
[$]\exp(\int_0^t[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv]) \frac{X_t^{(1-\gamma)}}{(1-\gamma)}[$]
[$]=\frac{X_0^{(1-\gamma)}}{(1-\gamma)} + \int_0^t \exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])\sigma dz(s))[$]
From where we get the solution as
[$]\frac{X_t^{(1-\gamma)}}{(1-\gamma)}=\frac{X_0^{(1-\gamma)}}{(1-\gamma)}*\exp(\int_0^t[+\mu_1(1-\gamma)X_v^{\beta_1-1}dv+\mu_2(1-\gamma)X_v^{\beta_2-1}dv+.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv]) [$]
[$]+\exp(\int_0^t[+\mu_1(1-\gamma)X_v^{\beta_1-1}dv+\mu_2(1-\gamma)X_v^{\beta_2-1}dv+.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])[$]
[$]*\int_0^t \exp(\int_0^s[-\mu_1(1-\gamma)X_v^{\beta_1-1}dv-\mu_2(1-\gamma)X_v^{\beta_2-1}dv-.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}dv])\sigma dz(s))[$]
writing the integral in a simplified form, the above expression becomes
[$]\frac{X_t^{(1-\gamma)}}{(1-\gamma)}=\frac{X_0^{(1-\gamma)}}{(1-\gamma)}*\exp(\int_0^t K(X_v)dv) [$]
[$]+ \exp(\int_0^t K(X_v)dv) \int_0^t \exp( \int_0^s -K(X_v)dv ) \sigma dz(s) [$]
where [$]K(X_v)=+\mu_1(1-\gamma)X_v^{\beta_1-1}+\mu_2(1-\gamma)X_v^{\beta_2-1}+.5\sigma^2(-\gamma)(1-\gamma)X_v^{2\gamma-2}[$]
I can solve the integrals in my framework but the question is if there is a clever analytical method to solve these integrals.
I will be coming back with more posts in next few days. I hope to conclusively solve the problem now.

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

Posted: July 12th, 2019, 11:38 am
Here is the new improved program for simulation of densities using Ito-Taylor density Simulation method I have worked on. It works quite well. But it is still not perfect for large terminal time with very high volatilities in the SDE. I have worked with squared orthogonal increments method and removed instabilities in the earlier versions and made the algorithm more accurate. Here is the program. It computes the density analytically and then makes a comparison with monte carlo density.

function [] = ItoTaylorDensityMRRecombining07()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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 directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.
dt=.125/8;   % Simulation time interval.%For diffusions close to zero
%decrease dt for accuracy.
Tt=8*4;     % Number of simulation levels. Terminal time= Tt*dt; //.125/8*8*8=1 year;
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=40;  % No of normal density subdivisions
NnMidl=20;%One half left from mid of normal density
NnMidh=21;%One half right from the mid of normal density
NnMid=4.0;
x0=1.0;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.550;   % volatility power.
kappa=.750;   %mean reversion parameter.
theta=1.0;   %mean reversion target
%Below you can get rid of kappa and theta and freely choose your own values
%for both drift coefficients.
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.
y(1:Nn)=x0;
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-.5)*dNn-NnMid);

ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=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;
dtsqrt=sqrt(dt);
dtonehalf=dt.^1.5;

wnStart=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);
VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
%        d2GGVV(nn)=omega1*theta1*sigma0.*(omega1+gamma-1).*yy(nn).^(omega1+gamma-2)+ ...
%            omega2*theta2*sigma0.*(omega2+gamma-1).*yy(nn).^(omega2+gamma-2)+ ...
%            omega3*theta3*sigma0.*(omega3+gamma-1).*yy(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;

dw(wnStart:Nn)=sigma0*Z(wnStart:Nn).*dtsqrt+ ...
dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).*dtonehalf;%+ ...
%0*d2GGVV(nn).*VV(nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(nn).^2-1);

if(tt>0)
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
wPrev=w;
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
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)).*(dw(wnStart:Nn).^2)));

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;

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
%    EnterLoop1=1;
%    EnterLoop2=1;
%    for nn=Nn-1:-1:wnStart
%        %Makes sure values decrease monotonically close to zero.
%        if((w(nn)>w(nn+1))&&(EnterLoop1==1))
%            w(1:nn)=0;
%            wnStart=nn+1;
%            EnterLoop1=0;
%        end
%        %Makes sure values are not smaller than zero.
%        if((w(nn)<0)&&(EnterLoop2==1))
%            w(1:nn)=0.0;
%            wnStart=nn+1;
%            EnterLoop2=0;
%        end
%    end
end

end

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(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));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
end

toc
str=input('Analytic Values calcuated; 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:Tt+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) = dt^((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));
end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Calculate the density with monte carlo below.
rng(85109634, 'twister')
paths=100000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:Tt
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;
YY0=YY;
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)

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;
sum(YY(:))/paths %origianl coordinates monte carlo average.
sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075/2;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=30;
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('Blue graph is from squared orthogonal incrments method and red is from squared volatility, green is monte carlo.');

end

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

Posted: July 13th, 2019, 8:49 am
As the friends would notice that the above program is quite accurate for reasonable volatilities but for very high volatilities, the results start to deteriorate after a year and more. On the other hand, well written PDEs give precise results even for extremely high volatilities. Therefore I decided to add PDEs to my existing research repertoire and substitute numerical solution of Fokker Planck where it is easily available and continue with the analytics. Once a precisely calculated density on a grid is available, it does not matter for further analytics if we obtained the grid from Ito-Taylor density generation method or from numerical solution of a PDE. I will continue to try to make the results of my Ito-Taylor density method sharply precise for extreme volatilities but in the mean time would also would do analytics with numerical solution of forward Fokker-Planck PDEs. Here is a post I made on Linkedin that elaborates on this idea.

Introducing New Technologies in Mathematical Finance Discipline. A Totally New Paradigm.

I want to mention some interesting research achievements at Infiniti Derivatives Technologies in the area of stochastics leading to huge improvements in mathematical methods of derivatives valuation, hedging, trading and risk management.

1.      Valuation of Path Dependent Stochastic Integrals on a PDE Grid.
It has been common historical wisdom that path dependent stochastic integrals cannot be calculated on a PDE grid. In our new research, we debunk these past notions and show that there are methods of introducing path dependence on a PDE grid. We have introduced new technologies that can calculate any stochastic path integral on a one or more dimensional PDE grid in fraction of a second after the PDE has been solved. Of course calculating path integrals on higher dimensional PDEs would require larger time but the time required to solve stochastic path integrals would strictly be far lesser than required to solve the PDE itself. Some of these stochastic integrals could be for example conditional value of path integral of squared volatility on the PDE grid for valuation of derivatives with stochastic volatility. Once the PDE has been solved conditional backward stochastic integrals like exponential of path integral of interest rates representing conditional bond prices or discount factor could also be calculated on the same PDE grid. From the PDE of interest rates, you could calculate the backward evolution of difficult stochastic path integrals depending on interest rates of whose PDE cannot be written easily. For Asian options, you could find the path integral of underlying as a forward stochastic integral on the Fokker planck equation solved by a one dimensional PDE or on a two dimensional PDE with stochastic volatility. The conditional values of these stochastic path integrals are available anywhere on the PDE grid along the backward or forward evolution of the PDE. PDE methods can only calculate backward evolution of payoffs whose explicit PDE has been written with arbitrary final data but we can calculate backward dynamic programming valuation of far more complex stochastic integrals dependent on PDE underlying of whose explicit PDE is not available but conditional of the value of the underlying as calculated in the forward PDE. We do not use any simulation or random numbers and only rely on fast analytics to achieve the above mentioned goals.

2.      Calculation of Optimal Quantizers for Monte Carlo Simulations.
Once a Fokker planck PDE has been solved, we can calculate optimal quantizers in one or two dimensions in fraction of a second. You could use these optimal quantizer schemes for forward evolution of stochastic integrals or backward dynamic programming in a monte carlo setting. As opposed to existing optimal quantizer methods that are based on monte carlo discretization schemes coupled with optimization techniques, we take a Fokker planck PDE solution as a starting point which is more cheaply available and simply use very fast proprietary analytics to generate optimal quantizer grid that can later be used for extremely fast monte carlo simulations and backward dynamic programming.

If you find any of this interesting and would like to implement this solution at your institution, you can email me at anan2999@yahoo.com for further information and discussions. I want to mention again that none of these new technologies are in public domain.

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

Posted: August 1st, 2019, 8:19 am
Over the weekend, I will be writing two posts. First would explain the algorithm I used to remove instability from the middle subdivisions in the density generation method and give a perfectly smooth output. This is important since many friends would like to alter the program in their own fashion and would like to know how to remove the instability later as the current method to remove instability is tailored to the grid I used and may fail on other grids. Once the basics of the method are known, friends would be able to fix any grid with this method to remove instability.
Second post will give another algorithm that will be used if you have a probability distribution that is very close to true distribution but the mean is slightly off. The algorithm would alter the probability distribution in a way that new probability distribution is a valid distribution that it integrates to one and its mean integrates to true mean. I will be doing both these things over the weekend.

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

Posted: August 1st, 2019, 8:36 am
I have decided to write a formal paper on the work I have done on this new density generation method for the diffusions/SDEs. The paper will also show how to calculate stochastic integrals of the SDEs on the same grid used to calculate the density in the first place. The paper will also show how to calculate the moments of these stochastic integrals on the same grid. For example these mth moments can be written as
[$]E[ [ \int_0^T f(x) dt ]^m] [$]  and [$]E[ [ \int_0^T g(x) dz ]^m] [$]  where we are given the SDE of x(t) and we will first evolve this SDE using our method and then calculate the stochastic integrals and then calculate moments of the same stochastic integrals. This will be a very interesting paper and I hope that it will take no more than a month to complete it and put it on SSRN. I might also add a section on backward dynamic programming for American options with non-linear payoffs.
I will be posting many more programs here on Wilmott in coming weeks.
Though the method is not perfect with mean reverting equations with very high volatility and high mean reversions, I do find solace in the fact that it works very well for a lot of other simpler SDEs. When we think of monte carlo simulations, we know that monte carlo has been there for more than half a century but we have been able to take care of nonlinear mean reversion SDE equations with monte carlo only very recently. So hopefully, other intelligent people will make more interesting contributions to solve non-linear mean reverting SDEs in a few years.

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

Posted: August 4th, 2019, 6:36 pm
In this new version of the program, I have also simulated the dt-integral given as [$]\int_0^t x(s)^{mp} ds[$]
for this I used the traditional formula given as
[$]\int_{t0}^{t1} x(s)^{mp} ds[$]
[$]=t1* x(t1)^{mp}- t0* x(t0)^{mp}[$]
[$]- mp * x(t0)^{(mp-1)}*(mu1*x(t0)^{beta1} + mu2*x(t0)^{beta2}) *(t1^2-t0^2)/2[$]
[$]-mp * x(t0)^{(mp-1)}*(sigma0 *x(t0)^{gamma})* N * (t1^{1.5}-t0^{1.5})/sqrt(3)[$]
[$] -.5*mp*(mp-1)*sigma0^2*x(t0)^{(mp+2*gamma-2)}*(t1^2-t0^2)/2.0[$]

where SDE for the evolution of x(t) is given as
[$]dx(t)=mu1* x(t)^{beta1} dt + mu2* x(t)^{beta2} dt +sigma0* x(t)^{gamma} dz(t)[$]

I have used the same dt-integral simulation formula in monte carlo and my transition probabilities framework. I have used original coordinates both in monte carlo and transition probabilities framework for calculation of dt-integrals.
Reference for the method for path-integrals/dt-integral I used in my program is my paper on SSRN titled "Solution of Stochastic Volatility Models Using Variance Transition Probabilities and Path Integrals"
https://ssrn.com/abstract=2149231 or
http://dx.doi.org/10.2139/ssrn.2149231
For diffusions going into zero, you may see instabilities, I have not corrected for that yet in this initial version.

function [] = ItoTaylorDensityMRRecombining14A()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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 directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.

%I have also simulated the dt-integral given as
% \int_0^t x(s)^mp ds  %%%here mp is power on x(s)
%for this I used the traditional formula given as
%\int_t0^t1 x(s)^mp ds=t1* x(t1)^mp- t0* x(t0)^mp - ...
% mp * x(t0)^(mp-1)*(mu1*x(t0)^beta1 + mu2*x(t0)^beta2) *(t1^2-t0^2)/2 -...
% mp * x(t0)^(mp-1)*(sigma0 *x(t0)^gamma)* N * (t1^1.5-t0^1.5)/sqrt(3) -...
%.5*mp*(mp-1)*sigma0^2*x(t0)^(mp+2*gamma-2)*(t1^2-t0^2)/2.0;

%I have used the same simulation formula in monte carlo and my
%transition probabilities framework. I have used original coordinates
%both in monte carlo and transition probabilities framework for calculation
% of dt-integrals.
% REference for the method for path-integrals/dt-integral I used in my
%program is my paper on
% SSRN titled "Solution of Stochastic Volatility Models Using Variance
%Transition Probabilities and Path Integrals"
%https://ssrn.com/abstract=2149231 or
%http://dx.doi.org/10.2139/ssrn.2149231
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.
dt=.125/8;   % Simulation time interval.%For diffusions close to zero
%decrease dt for accuracy.
Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/8*8*8=1 year;
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=40;  % No of normal density subdivisions
NnMidl=20;%One half left from mid of normal density
NnMidh=21;%One half right from the mid of normal density
NnMid=4.0;
x0=1.0;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.850;   % volatility power.
kappa=.950;   %mean reversion parameter.
theta=1.0;   %mean reversion target
%Below you can get rid of kappa and theta and freely choose your own values
%for both drift coefficients.
mu1=0*theta*kappa;   % first drift coefficient.
mu2=-0*kappa;    % Second drift coefficient.
sigma0=.450;%Volatility value
mp=.5; %power on x(t) in dt-integral. I have tested powers .5,1, and 2.0
%may fail for extrme high or low powers.
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
y(1:Nn)=x0;
w(1:Nn)=x0^(1-gamma)/(1-gamma);
x2dt0(1:Nn)=0.0;
x2dt(1:Nn,1:Tt)=0.0;
Z(1:Nn)=(((1:Nn)-.5)*dNn-NnMid);

ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=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;
dtsqrt=sqrt(dt);
dtonehalf=dt.^1.5;

wnStart=1;
tic
Zt(1:Nn,1:Nn)=0.0;

for tt=1:Tt

t0=(tt-1)*dt;
t1=tt*dt;
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);
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;%

dw(wnStart:Nn)=(sigma0.*Z(wnStart:Nn).^1.*dtsqrt+ ...
dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).^1.*dtonehalf + ... %;%+ ...
d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(wnStart:Nn).^2-1));
%I have added an extr term in above dw

w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.

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)).*((dw(wnStart:Nn)).^2)));

%Calculation of dt-integral starts
yyPrev(wnStart:Nn)=yy(wnStart:Nn);
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%Below we calculate transition probabilities between the process in
%original coordinates based on brownian motion transition probabilities.
%Below b1 is brownian motion subdivisions at time t0 and b2 are brownian
%motion subdivisions at time t1=t0+dt;
b1(1:Nn)=x0+sigma0*sqrt((tt-1)*dt)*Z(1:Nn);%Simulate brownian motion
%to the last date
b2(1:Nn)=x0+sigma0*sqrt(tt*dt)*Z(1:Nn);%Simulate brownian motion to
%recent date
for mm=1:Nn
Zt(mm,1:Nn)=(b2(1:Nn)-b1(mm))/sigma0/sqrt(dt);%Get the transition value
end
%Zt is the normal random number value which is used to mth subdivision
%of b1 at t0 to nth subdivision of b2 at t1. This will remain same for
%all the diffusions with volatility sigma0. And we later use the same
% normal random number value for our calculations of evolution of
%dt-integrals.

%dZ2(2:Nn-1)=(b2(3:Nn)-b2(1:Nn-2))/2.0;
%dZ2(1)=(b2(2)-b2(1));
%dZ2(Nn)=(b2(Nn)-b2(Nn-1));

dZ2(1:Nn)=sigma0*sqrt(t1)*dNn;

%Above dZ2 are subdivision widths. they are also equal to
%sigma0*sqrt(t1)*dNn. They are used to calculate integrated transition
%probability mass that arrives into a certain subdivision of b2 at time
%t1

Mm=Nn;
for nn=1:Nn
pt(nn,1:Mm)=exp(-.5* ((b2(1:Mm)-b1(nn))/(sigma0*sqrt(dt))).^2)./sqrt(2*pi*sigma0^2*dt).* ...
(dZ2(1:Mm));
%above pt is transition probability
end

%for nn=1:Nn
%Sumpt(nn)=sum(pt(nn,:));
%end
%Sumpt
if(tt==1)
pt0(1:Nn)=ZProb(1:Nn);
end

for mm1=1:Nn
for mm0=1:Nn

if(tt==1)
x2dt(mm1,tt)=x2dt(mm1,tt)+x2dt0(mm0).*pt0(mm1)+ ...
(yy(mm1).^mp*t1-yyPrev(mm0).^mp*t0- ...
(yy(mm1)^mp-yyPrev(mm0).^mp)*t0).*pt(mm0,mm1).*ZProb(mm0);
end
if(tt>1)

x2dt(mm1,tt)=x2dt(mm1,tt)+x2dt(mm0,tt-1).*pt(mm0,mm1)+ ...
(yy(mm1).^mp*t1-yyPrev(mm0).^mp*t0- ...
mp*(mu1*yyPrev(mm0).^(mp+beta1-1)+mu2*yyPrev(mm0).^(mp+beta2-1)).*(t1^2-t0^2)/2 - ...
.5*mp*(mp-1).*sigma0^2*yyPrev(mm0).^(mp+2*gamma-2).*(t1^2-t0^2)/2 - ...
mp*sigma0*yyPrev(mm0).^(mp+gamma-1) .*Zt(mm0,mm1).*(t1^1.5-t0^1.5)/sqrt(3)) .* ...
pt(mm0,mm1).*ZProb(mm0);

%This calculates the dt-integral using transition probabilities
%framework
end
end
end

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;

%
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;
%

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
end

xx2dt(1:Nn)=0.0;
xx2dt(1:Nn)=x2dt(1:Nn,Tt);%./ZProb(1:Nn);
xx2dt(1:Nn)=xx2dt(1:Nn)./ZProb(1:Nn); %We have to remove probability from
%the variable x2dt. In calculations, we had x2dt equal to
% dt-integral *probability of being in the particular subdivision.
y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfx2dt(1:Nn)=0;
for nn=2:Nn-1
Dfx2dt(nn) = (xx2dt(nn + 1) - xx2dt(nn - 1))/(Z(nn + 1) - Z(nn - 1));
end

fx2dt(1:Nn)=0;
for nn = 2:Nn-1
fx2dt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfx2dt(nn));%Origianl coordinates density
end
fx2dt
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
end

toc
str=input('Analytic Values calcuated; 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:Tt+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) = dt^((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));
end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,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;
Random1(1:paths)=0;
for tt=1:Tt
t1=tt*dt;
t0=(tt-1)*dt;
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;
YY0=YY;
YYPrev=YY;
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)

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
YY2dt(1:paths)=YY2dt(1:paths)+t1*YY(1:paths).^mp-t0*YYPrev(1:paths).^mp- ...
mp.*(mu1*YYPrev(1:paths).^(mp+beta1-1)+mu2*YYPrev(1:paths).^(mp+beta2-1)).*(t1^2-t0^2)/2.0 - ...
.5*mp*(mp-1).*sigma0^2*YYPrev(1:paths).^(mp+2*gamma-2).*(t1^2-t0^2)/2 - ...
mp*sigma0*YYPrev(1:paths).^(mp+gamma-1) .*HermiteP1(2,1:paths).*(t1^1.5-t0^1.5)/sqrt(3);
end
YY(YY<0)=0;

sum(YY2dt(:))/paths   %dt-integral average from monte carlo
sum(xx2dt(2:Nn-1).*ZProb(2:Nn-1)) %dt-integral from transition
%probabilities framework.

sum(YY(:))/paths %origianl coordinates monte carlo average.
sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075*2;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=30;
[Y2dtDensity,IndexOutY2dt,IndexMaxY2dt] = MakeDensityFromSimulation_Infiniti(YY2dt,paths,BinSize,MaxCutOff);
plot(xx2dt(2:Nn-1),fx2dt(2:Nn-1),'r',IndexOutY2dt(1:IndexMaxY2dt),Y2dtDensity(1:IndexMaxY2dt),'g');
str=input('red line is dt-integral from transition probability framework , green is monte carlo.');

BinSize=.0075/1;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=30;
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('red line is density of SDE from squared orthogonal incrments method, green is monte carlo.');

end

In a day or two, I will also post algorithm for removing instability from my programs as I promised.
Please note that transition probabilities dt-integral/Ito-integral density calculation program posted only works when diffusion remains strictly positive. This is a programming problem and a limitation only in the introductory version and in later versions, I will fix this problem and remove the restriction.  The part of the program that calculates the density of SDE works in almost all cases even when the density hits zero, it is only the new transition probabilities Ito-Integral part that is is currently limited to positive domain.  The algorithm to calculate the density of SDEs is 99.999% perfect when the diffusion takes just one term in the drift. Only when a second term is introduced in the drift and there are non-linear dynamics that the program has problems for highly non-linear equations.

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

Posted: August 23rd, 2019, 12:16 pm
I have been able to alter the new method to find solution to densities of mean reverting SDEs with very high volatility and very high mean reversions to perfect precision. I will be posting the matlab program over the weekend after properly fromatting the program and making notes/comments. I Was able to add second order effect from the existing volatility. After addition of second order volatility effect, we get perfect fit to shape of the monte carlo overlaid density and the mean is remarkably exact to within a few basis points from the true mean. I will be posting the program over the weekend. This is a very remarkable discovery and I am sure friends would love the new program.

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

Posted: August 25th, 2019, 5:29 pm
Sorry Double Post. Deleted.