Serving the Quantitative Finance Community

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

I am attaching the program for path integrals. Please note that it will work only for CEV type noises. I hope to post another program in a week that will solve for the path integrals of drifting and mean-reverting SDEs.
.
.
function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral01()

%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
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*2;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;

%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end

dtM=dt;
TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46*2;  % No of normal density subdivisions

NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.00;%mean reversion target
sigma0=.650;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
%yy(1:Nn)=x0;
w0=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-6.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

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

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.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

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

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
tt

if(tt==1)
%dt=ds(1);
%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
%       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);

[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);

%  dwMu0dtdw
%  dwMu0dtdwb
%  wMu0dt
%  wMu0dtb
%  str=input('Look at drift');
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;

elseif(tt<=6)

%dt=ds(tt);

%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
%[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
c0=c0+b0;
c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
c(2)=b(2);
c(3)=b(3);
c(4)=b(4);
c(5)=b(5);

w0=c0;

else

%below wMu0dt is drift and dwMu0dtdw are its derivatives.
%  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%size(b)
%ba(1:6)=b(1:6);

[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

%[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);

[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);

a0=c0+b0;

a(1:5)=c(1:5)+b(1:5);

c0=a0;
c(1)=a1;
c(2)=a2;
c(3)=a3;
c(4)=a4;
c(5)=a5;

w0=c0;

end

%Below Calculate the Z-series expansion for SDE variable in original coordinates
yw0=((1-gamma).*c0).^(1/(1-gamma));
dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;

if(tt>1)
fH0prev=fH0;
fHprev=fH;
end

%Below f0 and f define the Z_series in original coordinates
[f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
%fH0 and fH(1:5) define the coefficients of hermite polynomial
%representation of SDE variable in original coordinates.
[fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrder);

%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step. The take double array index in
%simulation level, tt and hermite polynomial count (1:5).
if(tt==1)
cH0=fH0;
cH(tt,1:SeriesOrder)=fH(1:SeriesOrder);

end
if(tt>1)
dH0=fH0-fH0prev;

dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-fHprev(1:SeriesOrder)).*sqrt(abs(sign(fH(1:SeriesOrder)).*fH(1:SeriesOrder).^2-sign(fHprev(1:SeriesOrder)).*fHprev(1:SeriesOrder).^2));

cH0(tt)=dH0;
cH(tt,1:SeriesOrder)=dH(1:SeriesOrder);
end

end
wnStart=1;

yy0=((1-gamma).*w0).^(1/(1-gamma));

w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Below calculate the hermite representation of dt-integral in original
%coordinates
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step and bH0 and bH define the hermite
%representation of the dt-integral after continued squared addition of
%scaled/weighted values of hermite coefficients of orthogonal increments.

bH0=cH0(1)*dt*Ts;
bH(1:SeriesOrder)=cH(1,1:SeriesOrder)*dt*Ts;

for tt=2:Ts
ts=Ts-tt+1;
bH0=bH0+cH0(tt)*dt*ts;
bH(1:SeriesOrder)=sign(bH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(bH(1:SeriesOrder)).*bH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));
end

%We convert hermite represenation given by bH0 and bH(1:5) into
%dt-integral Z-series Coeffs, e0 and e(1:5)

[e0,e] = ConvertHCoeffsToZCoeffs(bH0,bH,SeriesOrder);
yydt(1:Nn)=e0;
for nn=1:ExpnOrder
yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
end

%c0
%c

%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);

pyy(1:Nn)=0;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0

theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;

Random1(1:paths)=0;
for tt=1:TtM

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;

YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average

MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

end
.
.
.
function [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,a,SeriesOrder)

b0=wMu0dt;

b(1)=a(1) *dwMu0dtdw(1);

b(2)=1/2*(2* a(2) *dwMu0dtdw(1)+a(1)^2 *dwMu0dtdw(2));

b(3)=1/6*(6 *a(3) *dwMu0dtdw(1)+6 *a(1)* a(2) *dwMu0dtdw(2)+a(1)^3 *dwMu0dtdw(3));

b(4)=1/24*(24* a(4) *dwMu0dtdw(1)+12* a(2)^2 *dwMu0dtdw(2)+24 *a(1)* a(3) *dwMu0dtdw(2)+ ...
12* a(1)^2* a(2) *dwMu0dtdw(3)+a(1)^4 *dwMu0dtdw(4));

if(SeriesOrder>=5)
b(5)=1/120*(120* a(5) *dwMu0dtdw(1)+120* a(2) *a(3) *dwMu0dtdw(2)+120* a(1)* a(4) *dwMu0dtdw(2)+ ...
60 *a(1)* a(2)^2 *dwMu0dtdw(3)+60* a(1)^2* a(3) *dwMu0dtdw(3)+ ...
20 *a(1)^3* a(2) *dwMu0dtdw(4)+a(1)^5 *dwMu0dtdw(5));
end
if(SeriesOrder>=6)
b(6)=1/720*(720* a(6) *dwMu0dtdw(1)+360* a(3)^2 *dwMu0dtdw(2)+720* a(2)* a(4) *dwMu0dtdw(2)+720* a(1)* a(5) *dwMu0dtdw(2)+ ...
120* a(2)^3 *dwMu0dtdw(3)+720* a(1)* a(2)* a(3) *dwMu0dtdw(3)+360* a(1)^2* a(4) *dwMu0dtdw(3)+ ...
180* a(1)^2* a(2)^2 *dwMu0dtdw(4)+120* a(1)^3* a(3) *dwMu0dtdw(4)+ ...
30* a(1)^4 *a(2) *dwMu0dtdw(5)+a(1)^6 *dwMu0dtdw(6));
end

if(SeriesOrder>=7)
b(7)=1/(720*7)*(5040* a(7)*dwMu0dtdw(1)+5040* a(3)* a(4)*dwMu0dtdw(2)+5040* a(2)* a(5) *dwMu0dtdw(2)+5040 *a(1)* a(6) *dwMu0dtdw(2)+ ...
2520* a(2)^2* a(3) *dwMu0dtdw(3)+2520* a(1)* a(3)^2 *dwMu0dtdw(3)+5040 *a(1)* a(2)* a(4) *dwMu0dtdw(3)+2520* a(1)^2 *a(5) *dwMu0dtdw(3)+ ...
840* a(1)* a(2)^3 *dwMu0dtdw(4)+2520* a(1)^2* a(2)* a(3) *dwMu0dtdw(4)+840* a(1)^3* a(4) *dwMu0dtdw(4)+ ...
420* a(1)^3* a(2)^2 *dwMu0dtdw(5)+210* a(1)^4* a(3) *dwMu0dtdw(5)+ ...
42* a(1)^5* a(2) *dwMu0dtdw(6)+a(1)^7 *dwMu0dtdw(7));

end

if(SeriesOrder>=8)
b(8)=1/(720*7*8)*(40320* a(8) *dwMu0dtdw(1)+20160* a(4)^2 *dwMu0dtdw(2)+40320* a(3)* a(5) *dwMu0dtdw(2)+40320* a(2) *a(6) *dwMu0dtdw(2)+ ...
40320* a(1) *a(7) *dwMu0dtdw(2)+ ...
20160* a(2)* a(3)^2 *dwMu0dtdw(3)+20160* a(2)^2* a(4) *dwMu0dtdw(3)+40320* a(1)* a(3)* a(4) *dwMu0dtdw(3)+40320* a(1)* a(2)* a(5) *dwMu0dtdw(3)+ ...
20160* a(1)^2* a(6) *dwMu0dtdw(3)+ ...
1680* a(2)^4 *dwMu0dtdw(4)+20160* a(1)* a(2)^2* a(3) *dwMu0dtdw(4)+10080* a(1)^2* a(3)^2 *dwMu0dtdw(4)+20160* a(1)^2* a(2) *a(4) *dwMu0dtdw(4)+ ...
6720* a(1)^3* a(5) *dwMu0dtdw(4)+ ...
3360* a(1)^2* a(2)^3 *dwMu0dtdw(5)+6720* a(1)^3 *a(2)* a(3) *dwMu0dtdw(5)+1680* a(1)^4* a(4) *dwMu0dtdw(5)+ ...
840* a(1)^4* a(2)^2 *dwMu0dtdw(6)+336* a(1)^5* a(3) *dwMu0dtdw(6)+ ...
56* a(1)^6* a(2) *dwMu0dtdw(7)+a(1)^8 *dwMu0dtdw(8));
end
%b(9:10)=0.0;

%Bound=abs(b(1)/b0)*gamma/(1-gamma);
%Bound=1/abs(b0/b(1))*gamma;%/(1-gamma);

%for nn=1:9
%
%    if(abs(b(nn+1))>abs(b(nn)*Bound))
%  %      b(nn+1)=sign(b(nn+1))*abs(b(nn))*Bound;
%    end
%end

%[b0,b] = ApplyBoundsOnSeriesCoeffs(b0,b);

%
%
%
% b(9)=1/(720*7*8*9)*(362880* a(9) *dwMu0dtdw(1)+ ...
%     362880* a(4)* a(5) *dwMu0dtdw(2)+ 362880* a(3)* a(6) *dwMu0dtdw(2) +362880* a(2)* a(7) *dwMu0dtdw(2) +362880* a(1) *a(8) *dwMu0dtdw(2) + ...
%     60480* a(3)^3 *dwMu0dtdw(3)+362880* a(2)* a(3)* a(4) *dwMu0dtdw(3)+181440* a(1)* a(4)^2 *dwMu0dtdw(3)+181440* a(2)^2* a(5) *dwMu0dtdw(3) + ...
%     362880* a(1)* a(3)* a(5) *dwMu0dtdw(3)+362880 *a(1)* a(2)* a(6) *dwMu0dtdw(3) +181440 *a(1)^2* a(7) *dwMu0dtdw(3)+ ...
%     60480* a(2)^3* a(3) *dwMu0dtdw(4)+181440* a(1)* a(2)* a(3)^2 *dwMu0dtdw(4)+181440* a(1)* a(2)^2 *a(4) *dwMu0dtdw(4)+ ...
%     181440* a(1)^2* a(3)* a(4) *dwMu0dtdw(4)+181440* a(1)^2* a(2)* a(5) *dwMu0dtdw(4)+60480* a(1)^3* a(6) *dwMu0dtdw(4)+ ...
%     15120* a(1)* a(2)^4 *dwMu0dtdw(5)+90720* a(1)^2* a(2)^2* a(3) *dwMu0dtdw(5)+30240* a(1)^3* a(3)^2 *dwMu0dtdw(5)+ ...
%     60480 *a(1)^3* a(2)* a(4) *dwMu0dtdw(5)+15120* a(1)^4* a(5) *dwMu0dtdw(5)+ ...
%     10080* a(1)^3* a(2)^3 *dwMu0dtdw(6)+15120* a(1)^4 *a(2)* a(3) *dwMu0dtdw(6)+3024* a(1)^5* a(4) *dwMu0dtdw(6)+ ...
%     1512* a(1)^5* a(2)^2 *dwMu0dtdw(7)+504* a(1)^6* a(3) *dwMu0dtdw(7)+ ...
%     72* a(1)^7 *a(2) *dwMu0dtdw(8)+a(1)^9 *dwMu0dtdw(9));
%
%
%
%
% b(10)=1/(720*7*8*9*10)*(3628800* a(10)* dwMu0dtdw(1)+ ...
%     1814400 * a(5)^2 *dwMu0dtdw(2)+3628800 *a(4)* a(6)* dwMu0dtdw(2)+3628800* a(3)* a(7)* dwMu0dtdw(2)+3628800 *a(2)* a(8)* dwMu0dtdw(2)+ ...
%     3628800* a(1)* a(9)* dwMu0dtdw(2)+ ...
%     1814400* a(3)^2* a(4)* dwMu0dtdw(3)+1814400* a(2)* a(4)^2* dwMu0dtdw(3)+3628800* a(2)* a(3)* a(5)* dwMu0dtdw(3)+3628800* a(1)* a(4)* a(5) *dwMu0dtdw(3)+ ...
%     1814400* a(2)^2* a(6)* dwMu0dtdw(3)+3628800* a(1)* a(3)* a(6)* dwMu0dtdw(3)+3628800 *a(1)* a(2)* a(7)* dwMu0dtdw(3)+1814400* a(1)^2* a(8)* dwMu0dtdw(3)+ ...
%     907200* a(2)^2 *a(3)^2* dwMu0dtdw(4)+604800* a(1)* a(3)^3* dwMu0dtdw(4)+604800* a(2)^3* a(4)* dwMu0dtdw(4)+3628800* a(1)* a(2)* a(3)* a(4)* dwMu0dtdw(4)+ ...
%     907200* a(1)^2* a(4)^2* dwMu0dtdw(4)+1814400* a(1)* a(2)^2* a(5)* dwMu0dtdw(4)+1814400* a(1)^2* a(3)* a(5)* dwMu0dtdw(4)+ ...
%     1814400* a(1)^2* a(2)* a(6)* dwMu0dtdw(4)+ 604800* a(1)^3* a(7) *dwMu0dtdw(4)+ ...
%     30240 *a(2)^5* dwMu0dtdw(5)+604800* a(1)* a(2)^3* a(3)* dwMu0dtdw(5)+907200 *a(1)^2* a(2)* a(3)^2* dwMu0dtdw(5)+907200* a(1)^2* a(2)^2* a(4) *dwMu0dtdw(5)+ ...
%     604800* a(1)^3* a(3)* a(4)* dwMu0dtdw(5)+604800* a(1)^3* a(2)* a(5)* dwMu0dtdw(5)+151200* a(1)^4* a(6)* dwMu0dtdw(5)+ ...
%     75600* a(1)^2* a(2)^4* dwMu0dtdw(6)+302400 *a(1)^3* a(2)^2* a(3)* dwMu0dtdw(6)+75600* a(1)^4 *a(3)^2* dwMu0dtdw(6)+151200* a(1)^4* a(2)* a(4) *dwMu0dtdw(6)+ ...
%     30240* a(1)^5* a(5)* dwMu0dtdw(6)+ ...
%     25200* a(1)^4* a(2)^3* dwMu0dtdw(7)+30240* a(1)^5* a(2)* a(3)* dwMu0dtdw(7)+5040* a(1)^6* a(4)* dwMu0dtdw(7)+ ...
%     2520* a(1)^6* a(2)^2* dwMu0dtdw(8)+720* a(1)^7* a(3)* dwMu0dtdw(8)+ ...
%     90 *a(1)^8* a(2)* dwMu0dtdw(9)+ ...
%     a(1)^10* dwMu0dtdw(10));

end


.
.
.
function [aH0,aH] = ConvertZCoeffsToHCoeffs(a0,a,SeriesOrder)

if(SeriesOrder==4)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3);
aH(2)=a(2)+6*a(4);
aH(3)=a(3);
aH(4)=a(4);
end

if(SeriesOrder==5)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4);
aH(3)=a(3)+10*a(5);
aH(4)=a(4);
aH(5)=a(5);
end

if(SeriesOrder==6)
aH0=a0+a(2)+3*a(4)+15*a(6);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4)+45*a(6);
aH(3)=a(3)+10*a(5);
aH(4)=a(4)+15*a(6);
aH(5)=a(5);
aH(6)=a(6);
end

end


.
.
.
function [a0,a] = ConvertHCoeffsToZCoeffs(aH0,aH,SeriesOrder)

if(SeriesOrder==4)
a0=aH0-aH(2)+3*aH(4);
a(1)=aH(1)-3*aH(3);
a(2)=aH(2)-6*aH(4);
a(3)=aH(3);
a(4)=aH(4);
end

if(SeriesOrder==5)
a0=aH0-aH(2)+3*aH(4);
a(1)=aH(1)-3*aH(3)+15*aH(5);
a(2)=aH(2)-6*aH(4);
a(3)=aH(3)-10*aH(5);
a(4)=aH(4);
a(5)=aH(5);
end

if(SeriesOrder==6)
a0=aH0-aH(2)+3*aH(4)-15*aH(6);
a(1)=aH(1)-3*aH(3)+15*aH(5);
a(2)=aH(2)-6*aH(4)+45*aH(6);
a(3)=aH(3)-10*aH(5);
a(4)=aH(4)-15*aH(6);
a(5)=aH(5);
a(6)=aH(6);
end

end


.
.
.
Here is the graph you should get for the path integral when you run the program. Parameters (changeable) that are hard-coded in the program are given on the graph title.
.
.

Here is the output you should see at the end of the program.

Elapsed time is 0.135680 seconds.

ItoHermiteMean =

1.000709620804804

true Mean only applicable to standard SV mean reverting type models otherwise disregard

TrueMean =

1

Original process average from monte carlo

MCMean =

1.001138749614760 - 0.000000001921218i

Original process average from our simulation

ItoHermiteMean =

1.000709620804804

ItoHermiteMeanPI =

2.000219732720310

MCMeanPI =

1.999724593240335 + 0.000000000230177i

true Mean only applicble to standard SV mean reverting type models otherwise disregard

TrueMean =

1

IndexMax =

732

red line is density of SDE from Ito-Hermite method, green is monte carlo.

IndexMax =

367

red line is density of SDE from Ito-Hermite method, green is monte carlo.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I am quickly posting the program for drifting SDEs. In this program, drift is calculated in original coordinates and then orthogonal increments are calculated after accounting for drift i.e. after adding drift(in a particular interval) linearly to the Z-series representation of SDE at the start of the interval. And in the end copies of drift are linearly added after multiplying them time left to terminal time. Orthogonal increments are added in a squared fashion like we did for CEV noises. Sometimes I noticed that there was some slack in representation of drift at the extreme tails but for small to medium drift the program seems to work reasonably out to a few years (Later I want to work on improving the Z-series/hermite representation in original coordinates). The current program will still not work for mean-reverting SDEs.
In this program you cannot specify mean reverting parameters kappa and theta and you have to directly specify the value of one drift coefficient mu1 and its power beta1.
Here is the new program.
.
.
function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral02()

%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
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*2;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;

%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end

dtM=dt;
TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=47*2;  % No of normal density subdivisions

NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=.650;
beta2=0.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.
kappa=0.250;%.950;   %mean reversion parameter.
theta=1.00;%mean reversion target
sigma0=.450;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=.075;   %first drift coefficient.
mu2=0;%-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
%yy(1:Nn)=x0;
w0=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-7.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

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

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.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

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

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
tt

if(tt==1)
%dt=ds(1);
%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
%       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);

[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);

%  dwMu0dtdw
%  dwMu0dtdwb
%  wMu0dt
%  wMu0dtb
%  str=input('Look at drift');
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;
b0=wMu0dt;
b(1:5)=0;

elseif(tt<=3)

%dt=ds(tt);

%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
%[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
c0=c0+b0;
c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
c(2)=b(2);
c(3)=b(3);
c(4)=b(4);
c(5)=b(5);

w0=c0;

else

%below wMu0dt is drift and dwMu0dtdw are its derivatives.
%  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%size(b)
%ba(1:6)=b(1:6);

[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

%[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);

[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);

a0=c0+b0;

a(1:5)=c(1:5)+b(1:5);

c0=a0;
c(1)=a1;
c(2)=a2;
c(3)=a3;
c(4)=a4;
c(5)=a5;

w0=c0;

end

%Below Calculate the Z-series expansion for SDE variable in original coordinates
yw0=((1-gamma).*c0).^(1/(1-gamma));
dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;

if(tt>1)
fH0prev=fH0;%+bH0;
fHprev=fH;%+bH;
end

%Below f0 and f define the Z_series in original coordinates
[f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
%fH0 and fH(1:5) define the coefficients of hermite polynomial
%representation of SDE variable in original coordinates.
[fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrder);

[yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(yw0,YCoeff0,Fp,ds(tt),ExpnOrder)
[by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,f,SeriesOrder);
[byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrder);

by
by0
byH0
byH

%str=input('Look at original drift');

if(tt>1)
fH0prev=fH0prev+byH0;
fHprev=fHprev+byH;
end

%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step. The take double array index in
%simulation level, tt and hermite polynomial count (1:5).
if(tt==1)
cH0=fH0;
cH(tt,1:SeriesOrder)=fH(1:SeriesOrder);

bt0=by0;
bt(tt,1:SeriesOrder)=by(1:SeriesOrder);
end
if(tt>1)
dH0=fH0-fH0prev;

dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-fHprev(1:SeriesOrder)).*sqrt(abs(sign(fH(1:SeriesOrder)).*fH(1:SeriesOrder).^2-sign(fHprev(1:SeriesOrder)).*fHprev(1:SeriesOrder).^2));

cH0(tt)=dH0;
cH(tt,1:SeriesOrder)=dH(1:SeriesOrder);

bt0(tt)=by0;
bt(tt,1:SeriesOrder)=by(1:SeriesOrder);
end

end
wnStart=1;

yy0=((1-gamma).*w0).^(1/(1-gamma));

w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Below calculate the hermite representation of dt-integral in original
%coordinates
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step and bH0 and bH define the hermite
%representation of the dt-integral after continued squared addition of
%scaled/weighted values of hermite coefficients of orthogonal increments.

eH0=cH0(1)*dt*Ts;
eH(1:SeriesOrder)=cH(1,1:SeriesOrder)*dt*Ts;

for tt=2:Ts
ts=Ts-tt+1;
eH0=eH0+cH0(tt)*dt*ts;
eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));
end

%We convert hermite represenation given by bH0 and bH(1:5) into
%dt-integral Z-series Coeffs, e0 and e(1:5)

[e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrder);%

e0=e0+bt0(1)*dt*Ts;
e(1:SeriesOrder)=e(1:SeriesOrder)+bt(1:SeriesOrder)*dt*Ts;
for tt=2:Ts
ts=Ts-tt+1;
e0=e0+bt0(tt)*dt*ts;
e(1:SeriesOrder)=e(1:SeriesOrder)+bt(tt,1:SeriesOrder)*dt*ts;
%eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));
end

yydt(1:Nn)=e0;
for nn=1:ExpnOrder
yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
end

%c0
%c

%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);

pyy(1:Nn)=0;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0

theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;

Random1(1:paths)=0;
for tt=1:TtM

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;

YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average

MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

end
.
.
Here is the graph you should see when you run this program.
.
.

.
.
.
Here is the output you should see when you run the above program.
ItoHermiteMean =

1.153993084863467

true Mean only applicable to standard SV mean reverting type models otherwise disregard

TrueMean =

1

Original process average from monte carlo

MCMean =

1.154116780238912

Original process average from our simulation

ItoHermiteMean =

1.153993084863467

ItoHermiteMeanPI =

2.156851330823339

MCMeanPI =

2.151806582653016

true Mean only applicble to standard SV mean reverting type models otherwise disregard

TrueMean =

1

IndexMax =

378

red line is density of SDE from Ito-Hermite method, green is monte carlo.

IndexMax =

190

red line is density of SDE from Ito-Hermite method, green is monte carlo.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I am writing this post to make a few requests. Though it is still April, we are seeing reasonably high temperatures in Lahore. Summers are very difficult in this city. High temperature easily reaches 115 Fahrenheit on some days and there is a span of two months when night-time temperature remains above 90 F even in the early morning before dawn. City is a concrete jungle and traps heat and I would be sweating even in the early morning. For a large number of people this is not a problem since they spend most of their time in air-conditioned environment and travel in air-conditioned cars. But when I try to use air-conditioners, I am greeted with huge amount of mind control gases that make it totally impossible for me to spend time there. Last year, they allowed me to use AC in our drawing room where I would sleep and work but only for 7-10 days and then started using mind control gases so I had to stop using that AC.
For this reason, I have spent many years (>10 years) when I would sleep in the open balcony on first floor in the summers with a simple pedestal fan. And some times they would somehow put mind control gases in the fan and it would release mind control gases with the air and it would be extremely difficult to sleep even in the balcony and I would have to stop the fan and I would be completely wet with sweating. Then I started to use small chargeable fans at night that were difficult to temper with and I had enough air. But then my brother got married and had two children and we decided to give all of the upper portion to my brother's family (year 2014/2015)  and I had to stop sleeping in the balcony on first floor. But I still like to sleep in the semi-open covered garage when weather is somewhat pleasant outside. And that is how I had dengue last year. Summers are still not reaching their peak and it is still reasonably pleasant in the garage at night so I have been sleeping here for past two weeks. It is still not dengue season (that comes after monsoon ) so I am not very worried about dengue. So I want to request friends to ask and pressurize mind control agents to let me sleep in the air-conditioning without mind control gases.
Mind control agents have been increasing mind control for past few days. They drugged bottled water at several places before I would reach those shops and buy water. I came to know since water of both Nestle and Coca-Cola brands is generally good and is not drugged at manufacturing source, so it must have been drugged ahead of my arrival at those shops. In past 3-4 days , I was hit with drugged water several times even though Nestle and Coca-Cola are not drugging their water.
And really most painful thing is that they continue to force mind control chemical into my mouth during my sleep again and again. I have mentioned it several times before but they never end doing it. Just two days ago when (in the night) I would repeatedly walk into the lounge from garage( where I would be sleeping ) and brush my teeth/mouth during the night again and again, they placed some device inside the water supply of wash basin in the lounge and I noticed (after returning from outside) that water of the wash basin was ionized and was charging my mouth. So they started adding mind control chemicals to wash basin's water because I would repeatedly brush my mouth there in the midst of the night. This thing happened just two or three days ago. They are totally determined to charging my mouth and forcing chemicals into my mouth during my sleep.
I want to request friends to please force the mind control agents to be better and end their mind control.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, after writing the previous post, I went out to have food in the morning. I bought four bottles of 1.5 L Nestle water that turned out to be drugged since I bought them from a market close to my home.  I returned three bottles and bought water from elsewhere. After having food from a Gourmet bakery in Muslim town, I returned home. And on the way home, mind control agents gave me extremely strong anxiety attacks by forcing waves and charges through my eyes. I had received antipsychotic injection on Sunday and I am very vulnerable to anxiety attacks like this in a period of 3-10 days after the injections. These injections make me frail and mind control agents can easily play a lot of their tricks on me after that. I was supposed to take my mother to my uncle's house after returning from outside but I was feeling so bad that I requested her to excuse me and my sister had to take my mother to out uncle's house. I did not do any work and only lied on bed and also slept somewhat after that. And of course, mind control agents continued to force mind control chemicals into my mouth all through my sleep and I had to wash my mouth again and again. Even in the evening I was feeling very frail and at some point even had slight difficulty breathing and I was not having a full breath. When I realized that I was feeling short of breath, I decided to write this post.
I want to request friends that I am already 49 years old and very overweight due to antipsychotics that are forced on me. It is very difficult for me to live with antipsychotics at this age. I have been on antipsychotics for most of the time after start of year 1999 despite that I am perfectly mentally healthy. I want to request friends to please force American army and mind control agents to not force antipsychotics on me and let me live without mind control.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, below are a few pairs of graphs. The first graph in each pair shows a comparison of the SDE variable in original coordinates converted from Bessel coordinates on the Z-grid(Green graph), and analytically converted from eight cumulants(or expansion till seventh power of Z -- blue graph), and analytically converted from six cumulants (or expansion till fifth power of Z --black graph). You can see that Z-series representation of the SDE variable in original coordinates is less accurate as compared to original coordinates converted from Bessel coordinates on the Z-grid. You can see the graphs below to get an idea how accurate is Z-series representation from eight cumulants(in Bessel coordinates converted to original coordinates through series conversion) as compared to Z-series representation from six cumulants and both of the above compared with almost true representation(green graph). We like Z-series representation in original coordinates since we can do a large number of analytics with this representation.
Below are the graph comparisons. First graph in each pair shows the comparison of Z-series representation as a function of Z while second graph in each pair shows the comparison of density derived form each of the above three cases.

You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

I just un-liked more than ten posts on linkedin that were shown to be liked by me in my activity page. Some of them really deserved to be liked but I had not originally liked them myself so I un-liked them all. I had just liked two posts in past several months. Some of the content they had liked from my account was also religious stuff that I do not like.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, there were some minor errors in my program to calculate path integrals of drifting SDEs and the errors were causing  a decrease in accuracy. I will come back with an updated version in a few days.
I am busy with work on path integrals of mean-reverting SDEs. I did not make any progress for past three days and was ready to give up but today I realized there were multiple mistakes in my program and I also had some new ideas so I am optimistic again and hopeful that I would be able to do it in next few days.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I have been able to solve for dt-integrals of classical mean-reversion type SDEs. I hope to present the matlab program to solve for these integrals on Monday. I will also post an updated program for calculation of dt-integrals of drifting SDEs as the earlier posted program had some minor errors.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

I will be posting some graphs about my latest research later today or possibly tomorrow. To correct the previous post, I have been able to solve for the path integrals of SDEs where V0=theta i.e where the whole body of SDE is not drifting. When the whole body of SDE is drifting towards a mean-reversion target that is different from initial value of the SDE, I still have not been able to solve for path integrals. But I am quite confident that I would also be able to solve for these SDEs hopefully very soon.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

I am posting some graphs for path integrals of some mean-reverting SDEs in which $V_0=\theta$ i.e. SDEs where mean reversion target is the same as initial value of the SDE.
SDE parameters on title of the graphs.
.
.

You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I am posting some graphs of dt-integrals of drifting SDEs that take one drift term. The coefficient of drift term can be both positive or negative. Here is the relevant SDE.

$dX(t) \,= \, \mu_1 {X(t)}^{\beta_1} \, dt + \, {X(t)}^{\gamma} \, dZ(t) \,$

Most of the graphs are out to eight years at moderate vols and reasonable drifts. I am also going to post the matlab program for calculation of path integrals drifting SDEs in next post (that was used to create these graphs. It is still a preliminary program and I want a finessed version with comments later). The program for calculation of dt-integrals of mean-reverting SDEs will still take another day or two.

Here are the graphs of dt-integrals of drifting SDEs. The parameters are on the title of the graph.
.
.

You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Here is the matlab program that was used to create these graphs in previous post.
.
.
.
function [] = FPESeriesCoeffs_Lite_PIntegral08AB_Drifting()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%or skype ahsan.amin2999

%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*8;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;

%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end

dtM=dt;
TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=47*2;  % No of normal density subdivisions

NnMid=4.0;

x0=1;   % starting value of SDE

%beta1=0.0;
beta2=0.0;   % Second drift term power.
beta1=0.9950;
%beta2=0.0;   % Second drift term power.
gamma=.6;%50;   % volatility power.
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.250;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=-.150;   %first drift coefficient.
mu2=0;%-1*kappa;    % Second drift coefficient.
%mu1=theta*kappa;
%mu2=-1*kappa;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
%yy(1:Nn)=x0;
w0=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-7.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

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

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.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

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

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
tt

if(tt==1)
%dt=ds(1);
%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
%       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);

[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);

%  dwMu0dtdw
%  dwMu0dtdwb
%  wMu0dt
%  wMu0dtb
%  str=input('Look at drift');
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;
b0=wMu0dt;
b(1:5)=0;

elseif(tt<=3)

%dt=ds(tt);

%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
%[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
c0=c0+b0;
c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
c(2)=b(2);
c(3)=b(3);
c(4)=b(4);
c(5)=b(5);

w0=c0;

else

%below wMu0dt is drift and dwMu0dtdw are its derivatives.
%  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%size(b)
%ba(1:6)=b(1:6);

[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

%[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);

[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);

a0=c0+b0;

a(1:5)=c(1:5)+b(1:5);

c0=a0;
c(1)=a1;
c(2)=a2;
c(3)=a3;
c(4)=a4;
c(5)=a5;

w0=c0;

end

%Below Calculate the Z-series expansion for SDE variable in original coordinates
yw0=((1-gamma).*c0).^(1/(1-gamma));
dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;

if(tt==1)
bSumH0=0;%+bH0;
bSumH(1:SeriesOrder+2)=0;%+bH;
f0prev=x0;
fprev(1:SeriesOrder)=0;
end

if(tt>1)
fH0prev=fH0;%+bH0;
fHprev=fH;%+bH;
f0prev=f0;
fprev=f;
end

%Below f0 and f define the Z_series in original coordinates
[f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
f(6:8)=0;

SeriesOrdery=SeriesOrder;

%fH0 and fH(1:5) define the coefficients of hermite polynomial
%representation of SDE variable in original coordinates.
[fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrdery);

[yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(f0prev,YCoeff0,Fp,ds(tt),7)
[by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,fprev,SeriesOrdery);
[byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrdery);

%by
%by0
%byH0
%byH

fH(SeriesOrdery+1:7)=0;
byH(SeriesOrdery+1:7)=0;

%str=input('Look at original drift');

% if(tt>1)
%    fH0prev=fH0prev+byH0;
%    fHprev=fHprev+byH;
% end

if(tt>=1)

bSumH0prev=bSumH0;
bSumHprev=bSumH;
%fH0prev=fH0prev+byH0;
%fHprev=fHprev;%-exp(-kappa*dt).*bSumH;
bSumH0=bSumH0+byH0;
bSumH=bSumH+byH;

end

%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step. The take double array index in
%simulation level, tt and hermite polynomial count (1:5).
if(tt==1)
cH0(tt)=fH0;
cH(tt,1:SeriesOrdery)=fH(1:SeriesOrdery)-byH(1:SeriesOrdery);

bt0=by0;
bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end
if(tt>1)
%dH0=fH0-fH0prev;
dH0=fH0-bSumH0-fH0prev+bSumH0prev;

%dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-fHprev(1:SeriesOrdery)).*sqrt(abs(sign(fH(1:SeriesOrdery)).*fH(1:SeriesOrdery).^2-sign(fHprev(1:SeriesOrdery)).*fHprev(1:SeriesOrdery).^2));
dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));

cH0(tt)=dH0;
cH(tt,1:SeriesOrdery)=dH(1:SeriesOrdery);

bt0(tt)=by0;
bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end

end
wnStart=1;

yy0=((1-gamma).*w0).^(1/(1-gamma));

w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Below calculate the hermite representation of dt-integral in original
%coordinates
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step and bH0 and bH define the hermite
%representation of the dt-integral after continued squared addition of
%scaled/weighted values of hermite coefficients of orthogonal increments.

eH0=cH0(1)*dt*Ts;
%eH(1:SeriesOrdery)=cH(1,1:SeriesOrdery)*dt*Ts;

for tt=2:Ts
%   ts=Ts-tt+1;
%   eH0=eH0+cH0(tt)*dt*ts;
%   eH(1:SeriesOrdery)=sign(eH(1:SeriesOrdery)+dt*(ts).*cH(tt,1:SeriesOrdery)).*sqrt(abs(sign(eH(1:SeriesOrdery)).*eH(1:SeriesOrdery).^2+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2));
end

%eH2(1:SeriesOrdery)=sign(cH(1,1:SeriesOrdery)*dt*Ts).*(cH(1,1:SeriesOrdery)*dt*Ts).^2;

eH2(1:SeriesOrdery)=0;
eH0=0;

for tt=2:Ts
%   ts=Ts-tt+1;
%   eH0=eH0+cH0(tt)*dt*ts;
%   eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2;
%  for ss=1:tt-1
%      tss=tt-ss;
%      %tss=1;
%      %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;
%      eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*2.*tss.*cH(ss,1:SeriesOrdery).^2;
%  end
end
Coeff0(1:Ts)=0;
Coeff(1:Ts)=0;
Coeff1(1:Ts)=0;
Coeff2(1:Ts)=0;
%  for tt=1:Ts
%      Coeff2(tt)=Coeff2(tt)+0.*(Ts-tt+1).*dt^2;
%  end
for tt=1:Ts
for ss=tt:Ts
for rr=ss:Ts
if(rr==ss)

Coeff2(tt)=Coeff2(tt)+dt^2;
else

Coeff2(tt)=Coeff2(tt)+2.*dt^2;
end

end
end
end

for ss=1:Ts
for tt=ss:Ts

Coeff1(ss)=Coeff1(ss)+dt;
end
end

for tt=1:Ts
ts=Ts-tt+1;
eH0=eH0+cH0(tt)*Coeff1(tt);
eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(cH(tt,1:SeriesOrdery)).*Coeff2(tt).*cH(tt,1:SeriesOrdery).^2;
end

eH(1:SeriesOrdery)=sign(eH2(1:SeriesOrdery)).*sqrt(abs(eH2(1:SeriesOrdery)));

%We convert hermite represenation given by bH0 and bH(1:5) into
%dt-integral Z-series Coeffs, e0 and e(1:5)
%eH(SeriesOrdery)=eH(SeriesOrdery)/2;
[e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrdery);%

for ss=1:Ts
%for ss=1:tt
for tt=ss:Ts

e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(ss,1:SeriesOrdery)*dt;
e0=e0+bt0(ss)*dt;

end
end

yydt(1:Nn)=e0;
for nn=1:SeriesOrdery
yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
end

%c0
%c

%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);

pyy(1:Nn)=0;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0

theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;

Random1(1:paths)=0;
for tt=1:TtM

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;

%YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
YYdt(1:paths)=YYdt(1:paths)+YY(1:paths)*dtM;

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average

MaxCutOff=100;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

end
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I am posting some graphs for path integrals of mean-reverting SDEs where long term target theta is equal to the initial value of the SDE. The analytic calculation of dt-itnegrals matches with monte carlo versions quite well out to four or five years (for kappa 1.5 or 2.0) but accuracy starts to deteriorate after that as terminal time increases beyond five years. I am trying to find out what is the reason behind it.
Here are the graphs. In the next post, I will attach the program used to make these graphs.

.
.
.

You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Here is the program used to make graphs in previous post. It needs to be improved and I suspect it has some minor errors and I will come back with a new version after a few days.
.
.
.
function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral08CC()

%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
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*3;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;

%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end

dtM=dt;
TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=47*2;  % No of normal density subdivisions

NnMid=4.0;

x0=.10;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
%beta1=0.950;
%beta2=0.0;   % Second drift term power.
gamma=.85;%50;   % volatility power.
kappa=2.00;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.750;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=.50;   %first drift coefficient.
mu2=0;%-1*kappa;    % Second drift coefficient.
mu1=theta*kappa;
mu2=-1*kappa;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
%yy(1:Nn)=x0;
w0=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-7.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

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

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.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

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

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
tt

if(tt==1)
%dt=ds(1);
%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
%       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);

[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);

%  dwMu0dtdw
%  dwMu0dtdwb
%  wMu0dt
%  wMu0dtb
%  str=input('Look at drift');
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;
b0=wMu0dt;
b(1:5)=0;

elseif(tt<=3)

%dt=ds(tt);

%[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
%[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
c0=c0+b0;
c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
c(2)=b(2);
c(3)=b(3);
c(4)=b(4);
c(5)=b(5);

w0=c0;

else

%below wMu0dt is drift and dwMu0dtdw are its derivatives.
%  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%size(b)
%ba(1:6)=b(1:6);

[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);

%[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);

[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);

a0=c0+b0;

a(1:5)=c(1:5)+b(1:5);

c0=a0;
c(1)=a1;
c(2)=a2;
c(3)=a3;
c(4)=a4;
c(5)=a5;

w0=c0;

end

%Below Calculate the Z-series expansion for SDE variable in original coordinates
yw0=((1-gamma).*c0).^(1/(1-gamma));
dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;

if(tt==1)
bSumH0=0;%+bH0;
bSumH(1:SeriesOrder+2)=0;%+bH;
fH0prev=x0;%+bH0;
fHprev(1:SeriesOrder)=0;
f0prev=x0;
fprev(1:SeriesOrder)=0;
end

if(tt>1)
fH0prev=fH0;%+bH0;
fHprev=fH;%+bH;
f0prev=f0;%+bH0;
fprev=f;%+bH;
end

%Below f0 and f define the Z_series in original coordinates
[f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
f(6:8)=0;

SeriesOrdery=SeriesOrder;

%fH0 and fH(1:5) define the coefficients of hermite polynomial
%representation of SDE variable in original coordinates.
[fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrdery);

% [yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(f0prev,YCoeff0,Fp,ds(tt),7)
%  [by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,fprev,SeriesOrdery);
% [byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrdery);

% fH0
% byH0
% fH
%
% byH
% fH./byH
%
% 1/(exp(kappa*dt)-1)
% 1/(1-exp(-kappa*dt))
% byH./fH
% (1-exp(-kappa*dt))
%  (fHprev(1:SeriesOrder)+byH(1:SeriesOrder))./fHprev(1:SeriesOrder)
%  exp(-kappa*dt)
%  %fH+byH
%  %(fH+byH)./fH
%
% by0
% byH0
% byH

fH(SeriesOrdery+1:7)=0;
% byH(SeriesOrdery+1:7)=0;

% str=input('Look at original drift');

% if(tt>1)
%    fH0prev=fH0prev+byH0;
%    fHprev=fHprev+byH;
% end

if(tt>=1)

bSumH0prev=0;
bSumHprev(1:SeriesOrdery)=0;
%fH0prev=fH0prev+byH0;
%fHprev=fHprev;%-exp(-kappa*dt).*bSumH;
bSumH0=0;
bSumH(1:SeriesOrdery)=0;

end

%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step. The take double array index in
%simulation level, tt and hermite polynomial count (1:5).
if(tt==1)
cH0(tt)=fH0;
cH(tt,1:SeriesOrdery)=fH(1:SeriesOrdery);%-byH(1:SeriesOrdery);

% bt0=by0;
% bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end
if(tt>1)
%dH0=fH0-fH0prev;
dH0=fH0-bSumH0-fH0prev+bSumH0prev;

%dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-fHprev(1:SeriesOrdery)).*sqrt(abs(sign(fH(1:SeriesOrdery)).*fH(1:SeriesOrdery).^2-sign(fHprev(1:SeriesOrdery)).*fHprev(1:SeriesOrdery).^2));
dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(exp(0*kappa*dt).*fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(exp(-kappa*dt).*fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));
%dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-exp(-kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(exp(1*kappa*tt*dt).*fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*exp(2*(tt-1)*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));

cH0(tt)=dH0;
cH(tt,1:SeriesOrdery)=dH(1:SeriesOrdery);

% bt0(tt)=by0;
% bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end

end
wnStart=1;

yy0=((1-gamma).*w0).^(1/(1-gamma));

w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Below calculate the hermite representation of dt-integral in original
%coordinates
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
%increments at each simulation step and bH0 and bH define the hermite
%representation of the dt-integral after continued squared addition of
%scaled/weighted values of hermite coefficients of orthogonal increments.

eH0=cH0(1)*dt*Ts;
%eH(1:SeriesOrdery)=cH(1,1:SeriesOrdery)*dt*Ts;

for tt=2:Ts
%   ts=Ts-tt+1;
%   eH0=eH0+cH0(tt)*dt*ts;
%   eH(1:SeriesOrdery)=sign(eH(1:SeriesOrdery)+dt*(ts).*cH(tt,1:SeriesOrdery)).*sqrt(abs(sign(eH(1:SeriesOrdery)).*eH(1:SeriesOrdery).^2+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2));
end

%eH2(1:SeriesOrdery)=sign(cH(1,1:SeriesOrdery)*dt*Ts).*(cH(1,1:SeriesOrdery)*dt*Ts).^2;

eH2(1:SeriesOrdery)=0;
eH0=0;

for tt=2:Ts
%   ts=Ts-tt+1;
%   eH0=eH0+cH0(tt)*dt*ts;
%   eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2;
%  for ss=1:tt-1
%      tss=tt-ss;
%      %tss=1;
%      %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;
%      eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*2.*tss.*cH(ss,1:SeriesOrdery).^2;
%  end
end
Coeff0(1:Ts)=0;
Coeff(1:Ts)=0;
Coeff1(1:Ts)=0;
Coeff2(1:Ts)=0;
for tt=1:Ts
Coeff2(tt)=Coeff2(tt)+0.*(Ts-tt+1).*dt^2;
end
for tt=1:Ts
for ss=tt:Ts
for rr=ss:Ts
%Coeff2(tt)=Coeff2(tt)+0*.75* exp(-.5*kappa*(dt*((rr)-ss))).*dt^2;
%Coeff2(tt)=Coeff2(tt)+2/.77* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
%coef=max(.8-.8/sqrt(1)/kappa,0);
%coef=max(.75-.8/sqrt(2)/kappa,0);
%coef=max(.75-.8/sqrt(4)/sqrt(kappa)+2/3*sqrt(kappa)*sqrt(Ts*dt)-sqrt(kappa)*sqrt(.5),0);
%Coeff2(tt)=Coeff2(tt)+(2+coef)* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
%Coeff2(tt)=Coeff2(tt)+(2+coef)* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
if(rr==ss)

Coeff2(tt)=Coeff2(tt)+(exp(-0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;
else
%Coeff2(tt)=Coeff2(tt)+sqrt(3)/sqrt(kappa)*(exp(-2*kappa*dt*(rr-ss))).*dt^2;
Coeff2(tt)=Coeff2(tt)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt))).*dt^2;
end
%  rr
%  ss
%  (exp(-1*kappa*dt*(rr-ss)))
%  (exp(-2*kappa*dt*(rr-ss)))
%  str=input('Look at numbers');
%Coeff2(tt)=Coeff2(tt)+(2)* (exp(-kappa*dt*(rr-ss)))*(1+.125*kappa*dt*(ss)).*dt^2;
%Coeff2(tt)=Coeff2(tt)+(2)* (exp(-kappa*dt*(rr-ss)))*(1/.8-.125/4*kappa*dt*(tt)/2).*dt^2;
%Coeff2(tt)=Coeff2(tt)+1.*(exp(-kappa*dt*(tt/2))).*dt^2;

%Coeff2(tt)=Coeff2(tt)+1/(2).* dt^2;
end
end
end

for tt=1:Ts
ts=Ts-tt+1;
eH0=eH0+cH0(tt)*ts*dt;
eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(cH(tt,1:SeriesOrdery)).*Coeff2(tt).*cH(tt,1:SeriesOrdery).^2;
end

%   for ss=1:Ts-1
%         %tss=tt-ss;
%         tss=Ts-ss+1;
%         %tss=1;
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;
%         eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts+ss)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;
%  end

eH(1:SeriesOrdery)=sign(eH2(1:SeriesOrdery)).*sqrt(abs(eH2(1:SeriesOrdery)));

%We convert hermite represenation given by bH0 and bH(1:5) into
%dt-integral Z-series Coeffs, e0 and e(1:5)
%eH(SeriesOrdery)=eH(SeriesOrdery)/2;
[e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrdery);%

for ss=1:Ts
%for ss=1:tt
for tt=ss:Ts

%eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;
%eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;
%e(1:SeriesOrdery)=e(1:SeriesOrdery)+1/sqrt(2^0)*exp(-kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;
%        e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-1*kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;
%        e0=e0+exp(-kappa*dt*(tt-ss)).*bt0(ss)*dt;
%        e(1:SeriesOrdery)=0;
%        e0=0;

end
end

% for tt=1:Ts
%      for ss=tt:Ts
%
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;
%         %e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;
%         e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-kappa*dt*(ss-tt)).*bt(tt,1:SeriesOrdery)*dt;
%         %e0=e0+exp(-kappa*dt*(tt-ss)).*bt0(ss)*dt;
%         e0=e0+exp(-kappa*dt*(ss-tt)).*bt0(tt)*dt;
%      end
% end

%   e0=e0+bt0(1)*dt*Ts;
%   e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(1:SeriesOrdery)*dt*Ts;
%   for tt=2:Ts-1
%      ts=Ts-tt;
%      e0=e0+bt0(tt)*dt*ts;
%      e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(tt,1:SeriesOrdery)*dt*ts.*exp(-kappa*ts*dt);
%      %eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));
%   end

%e(5)=1/3*e(5);

yydt(1:Nn)=e0;
for nn=1:SeriesOrdery
yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
end

%c0
%c

%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);

pyy(1:Nn)=0;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0

theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;

Random1(1:paths)=0;
for tt=1:TtM

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;

YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average

MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,kappa=%.3f,theta=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

end
.
.
.
Here is the output you should see when you run this program.

ItoHermiteMean =

0.100021004761355

true Mean only applicable to standard SV mean reverting type models otherwise disregard

TrueMean =

0.100000000000000

Original process average from monte carlo

MCMean =

0.099800196176705

Original process average from our simulation

ItoHermiteMean =

0.100021004761355

ItoHermiteMeanPI =

0.299902402491998

MCMeanPI =

0.300090652148871

true Mean only applicble to standard SV mean reverting type models otherwise disregard

TrueMean =

0.100000000000000

IndexMax =

1145

red line is density of SDE from Ito-Hermite method, green is monte carlo.

IndexMax =

573

red line is density of SDE from Ito-Hermite method, green is monte carlo.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

Amin
Topic Author
Posts: 1838
Joined: July 14th, 2002, 3:00 am

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

Friends, I have been able to solve for the densities of CEV noises by addition of variances in bessel coordinates. The key is to keep drift and variance values distinct and apart from each other. We have to calculate drift at each point in time as a Z-series and find its sum across time. We also have to add volatility to previous volatility in a squared fashion (through hermite polynomials) and keep track of squared volatility/variance. Adding the sum of drift up to a certain time with square root of summed variance till that time helps us construct the Z-series of SDE variable at that time. Please note that we are dealing with CEV noises in bessel coordinates where their transformed SDE takes a drift term. This technique works well to simulate the CEV noises merely by adding variances.
When I tried this earlier last week, I was not keeping separate track of sum of drift and sum of squared volatility that seems necessary to properly evolve the SDE. Simulating CEV noise densities by adding variances makes their simulation very simple.
Though extension of this method to SDEs with one drift term is obvious, for mean-reverting SDEs, the distinction between drift and variance gets blurred so more work will be needed to simulate them by adding variances.
I will post my simple matlab program for simulation of CEV noises by adding variances later today or tomorrow.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal