Friends, very sorry for delay in the new program that works somewhat better close to zero. I am copying the condition I have imposed on coefficients from my matla
here a(nn) is nn'th coefficient of power series that multiplies Z^nn, Please feel free to play around with this growth restriction on absolute value of the coefficients.
I am posting the new program here.
function [] = FPESeriesCoeffsSolutionWilmottOdt2Downloaded()
%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*4; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*5/4;%128*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=.0625/1;
TtM=T/dtM;
dNn=.1/1; % Normal density subdivisions width. would change with number of subdivisions
Nn=46*2; % No of normal density subdivisions
%NnMidl=18;%One half density Subdivision left from mid of normal density(low)
%NnMidh=19;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
%NnT=46;
%NnD=(NnT-Nn)/2;
%NnD=5;
x0=.250; % starting value of SDE
beta1=0.0;
beta2=1.0; % Second drift term power.
gamma=.65;%50; % volatility power.
kappa=1.0;%.950; %mean reversion parameter.
theta=.05;%mean reversion target
sigma0=.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;
w0=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-6.5)*dNn-NnMid);
%Z(1:Nn)=(((1:Nn)+1.5)*dNn-NnMid);
%ZT(1:NnT)=(((1:NnT)-3.5)*dNn-NnMid);
%ZT
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=6;
wnStart=1;
tic
for tt=1:Tt
if(tt==1)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(dt)
c(2:10)=0.0;
w0=c0;
end
%The program intended to find solution of Z-series to 10th power but
%higher order coefficients would make some SDEs unstable so I just
%turned the coefficients c(7:10) and b(7:10) into zero as soon as they
%are calculated and just go with solution to 6th power of Z.
if(tt>1)
%below wMu0dta is drift and dwMu0dtdwa are its derivatives.
%wMu1dt is its first derivative and dwMu1dtdw are derivatives of first
%derivative.
%[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08Bounded(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
%The program CalculateDriftbCoeffs04 takes value of drift and its
%derivatives wrt w and also takes Z-series of w and converts them into
%Z-series of drift. It is used to calculate Z-series of drift and Then
%to calculate Z_series of first derivative of drift wrt w.
[b0,b] = CalculateDriftbCoeffs08(wMu0dta,dwMu0dtdwa,a,gamma);
[b10,b1] = CalculateDriftbCoeffs08(wMu1dt,dwMu1dtdw,a,gamma);
b(ExpnOrder+1:10)=0.0;
b1(ExpnOrder+1:10)=0.0;
%The function below evolves the previous solution in coefficients of Z
%into next time step solution in coefficietns of Z-series.
[c0,c] = EvolveSDESeriesSolution08New(a0,a,b0,b,b10,b1,sigma0,dt,ExpnOrder,gamma);
c(ExpnOrder+1:10)=0.0;
w0=c0;
end
a0=c0;
a(1:ExpnOrder)=c(1:ExpnOrder);
a(ExpnOrder+1:10)=0;
Bound=abs(a(1)/a0)*gamma/(1-gamma);
for nn=1:ExpnOrder-1
%if(abs(a(nn+1))>abs(a(nn))*Bound)
% a(nn+1)=sign(a(nn+1))*abs(a(nn))*Bound;
%end
if(abs(a(nn+1))>abs(a(nn))*2^(-nn+1))
a(nn+1)=sign(a(nn+1))*abs(a(nn))*2^(-nn+1);
end
end
end
c0=a0;
c=a;
wnStart=1;
w(wnStart:Nn)=c0;
for nn=1:ExpnOrder
w(wnStart:Nn)=w(wnStart:Nn)+c(nn)*Z(wnStart:Nn).^nn;
end
c0
c
%str=input('Look at numbers')
yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy0=((1-gamma).*w0).^(1/(1-gamma));
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
Dfyy(Nn)=Dfyy(Nn-1);
Dfyy(1)=Dfyy(2);
pyy(1:Nn)=0;
for nn = wnStart:Nn
pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(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*dt*Tt)%Mean reverting SDE original variable true average
yy0
theta1=1;
rng(29079137, 'twister')
paths=50000;
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:Nn).*ZProb(wnStart:Nn)) %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(1*900*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
.
.
Rest of the functions are old. When you run this program, you will get the output as
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
and the following graph.