Friends, here is the program with ODE evolution of densities of SDEs in original coordinates. I have done first 15 steps simulation in bessel coordinates since second derivative would be too small for stable evolution in original coordinates. Zero boundary is still not covered and friends would have to work on that a bit carefully. For densities that go into zero, program is still not appropriate in this basic form.
Again in this program, I have used the ODE for evolution of constant CDF points on the grid just as the ODE is derived in previous few posts.
Here is the program.
function [] = FPERevisitedTransProb08ABwNew04B()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)
%I have not directly simulated the SDE but simulated the transformed
%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/16/2/2/2; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*2*2*2*2; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
OrderA=4; %
OrderM=4; %
dtM=dt*4*2*4;
TtM=Tt/4/2/4;
dNn=.2/1; % Normal density subdivisions width. would change with number of subdivisions
Nn=50; % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
x0=1.0; % starting value of SDE
beta1=0.0;
beta2=1.0; % Second drift term power.
gamma=.95;%50; % volatility power.
kappa=4.0;%.950; %mean reversion parameter.
theta=.250;%mean reversion target
sigma0=1.50;%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;
w(1:Nn)=x0^(1-gamma)/(1-gamma);
x(1:Nn)=x0;
%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
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);
%Above calculate probability mass in each probability subdivision.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wnStart=1;%
tic
Zt1(wnStart:Nn)=0.0;
Zt2(wnStart:Nn)=0.0;
for tt=1:Tt
x(isnan(x)==1)=0.00;
[xMu0dt,c1] = CalculateDriftAndVolA4Original(x,wnStart,Nn,YCoeff0,Fp,gamma,dt);
%Loop below tackles first ten time steps in Bessel coordinates
if(tt<=10)
w(wnStart:Nn)=x(wnStart:Nn).^(1-gamma)/(1-gamma);
[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;
[dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);
Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dwdZ(wnStart:Nn)).^2+sigma0^2*dt)).*Z(wnStart:Nn);
%x2(wnStart:Nn)=((1-gamma)*(Zt2(wnStart:Nn)-0*Zt1(wnStart:Nn)+wMid)).^(1.0/(1-gamma));
end
[dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,x,Z);
% [xMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),x(NnMidl-3),x(NnMidl-2),x(NnMidl-1),x(NnMidl),x(NnMidh),x(NnMidh+1),x(NnMidh+2),x(NnMidh+3));
%
% Zt1(wnStart:Nn)=x(wnStart:Nn)-xMid;
%
% [dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,x,Z);
%
% C0(wnStart:Nn)=Zt1(wnStart:Nn)-dxdZ(wnStart:Nn).*Z(wnStart:Nn);%-.5*d2xdZ2A(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%
% %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dxdZ(wnStart:Nn)).^2+sigma0^2*x(wnStart:Nn).^(2*gamma).*dt)).*Z(wnStart:Nn);%+.5.*(sqrt((d2xdZ2A(wnStart:Nn).^2)+(sigma0^2.*(gamma).*x(wnStart:Nn).^(2*gamma-1)*dt).^2)).*(Z(wnStart:Nn).^2-1);
%
%
% %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dxdZ(wnStart:Nn).*x(wnStart:Nn).^gamma).^2+c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
% [dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
% x2=Zt2(wnStart:Nn)+xMid;
% tt
if(tt>10)
x(wnStart:Nn)= +x(wnStart:Nn)-sigma0.^2*gamma*x(wnStart:Nn).^(2*gamma-1)*dt ...
+xMu0dt(wnStart:Nn) ...
+.5*sigma0^2*x(wnStart:Nn).^(2*gamma).*(dxdZ(wnStart:Nn)).^(-2).*(d2xdZ2A(wnStart:Nn)).*dt ...
+.5*sigma0^2*x(wnStart:Nn).^(2*gamma).*Z(wnStart:Nn).*(dxdZ(wnStart:Nn)).^(-1)*dt;
else
%Bessel coordinates evolution associated with first ten steps.
w(wnStart:Nn)= wMid+Zt2(wnStart:Nn)+wMu0dt(wnStart:Nn);
x(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
%x(wnStart:Nn)= +x2(wnStart:Nn)-sigma0.^2*gamma*x(wnStart:Nn).^(2*gamma-1)*dt+ ...
% xMu0dt(wnStart:Nn);%+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);%+ ...
end
%yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
w(wnStart:Nn)=x(wnStart:Nn).^(1-gamma)/(1-gamma);
% [wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
%
% w1(wnStart:Nn-1)=w(wnStart:Nn-1);
% w1(Nn)=x(Nn);
% w2(wnStart:Nn-1)=w(wnStart+1:Nn);
% w2(Nn)=wE;
% w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
% % Change 3:7/25/2020: I have improved zero correction in above.
% w(w<0)=0.0;
for nn=wnStart:Nn
if(x(nn)<=0)
wnStart=nn+1;
end
end
end
%yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy(wnStart:Nn)=x(wnStart:Nn);
Dfyy(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
pyy(1:Nn)=0;
for nn = wnStart:Nn-1
pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
end
toc
ItoHermiteMean=sum(yy(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %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*dt*Tt)%Mean reverting SDE original variable true average
theta1=1;
rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0; %Original process monte carlo.
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;
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);
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+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
MaxCutOff=30;
NoOfBins=round(300*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.');
end
.
.
.
true Mean only applicable to standard SV mean reverting type models otherwise disregard
true Mean only applicble to standard SV mean reverting type models otherwise disregard
