Serving the Quantitative Finance Community

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

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

November 20th, 2021, 1:39 pm

Friends, I tried the new improved formula for simulation of constant CDF lines and densities of SDEs. It gives both higher accuracy everywhere and more stability closer to zero. But there is not a highly marked difference from the previous formula. We can easily increase the step size by two to four as compared to step size that we were earlier using with first order version but it is still difficult to increase the step size by a factor of ten or more. After the weekend, I will be posting a new program for simulation of one dimensional SDEs with the new advanced formula.
About the power series formula, its performance with first order was not very good and I want to try to see how it does with improved order. I have several different and meaningful ideas about power series formula and I hope that I would be able to make something interesting out of it but its results would also come out next week.
.
For most of the easy SDE diffusions where we earlier needed a step size of 1/512 or 1/256, we can now easily go by a step size of 1/100 and still not lose accuracy(This is with initial shorter stepping at start of diffusion simulation at delta origin). This is a good improvement but certainly not dramatic improvement. And closer to zero, you still need a relatively larger step for accuracy.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 21st, 2021, 4:25 pm

Friends, I was able to get the series recursive solution to work. I have retained only six powers in the series since higher order beyond six was more off with several SDEs and would need second  order time correction. My programs are still very initial and I have not done any higher 2nd order correction yet that I hope to do in next few days. Using this initial basic solution, I have drawn some graphs with various values of SDE parameters out to one year. Though the results are very encouraging, I believe there is huge room for improvement especially with second order correction that I hope to implement very soon. I really hope that final series solution will be vastly improved as compared to this version.
Here are the graphs from series solution. I will post the program used to create these graphs in the next post.
.
.
Image

Image

Image

Image

Image

Image

Image

Image
Last edited by Amin on November 21st, 2021, 4:36 pm, edited 2 times in total.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 21st, 2021, 4:34 pm

Here are the new files used to calculate the series solution. Please note that this is very basic initial solution and I intend to do a lot of effort to improve it to second order accuracy. I have not run the SDE solutions with this program beyond one year but I hope that with improved solution, we would easily be able to run it out to ten years.
.
function [] = FPESeriesCoeffsSolutionWilmott()

%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;%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=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % 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=1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.5;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=2.0;%.950;   %mean reversion parameter.
theta=.25;%mean reversion target
sigma0=.80;%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)-5.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);

%ZProbT(1)=normcdf(.5*ZT(1)+.5*ZT(2),0,1)-normcdf(.5*ZT(1)+.5*ZT(2)-dNn,0,1);
%ZProbT(NnT)=normcdf(.5*ZT(NnT)+.5*ZT(NnT-1)+dNn,0,1)-normcdf(.5*ZT(NnT)+.5*ZT(NnT-1),0,1);
%ZProbT(2:NnT-1)=normcdf(.5*ZT(2:NnT-1)+.5*ZT(3:NnT),0,1)-normcdf(.5*ZT(2:NnT-1)+.5*ZT(1:NnT-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

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

ExpnOrder=10;
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)

    %The function below calculates drift and its derivatives at w(Z=0);    
    [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
    
    a(7:10)=0.0;
    %The function below converts derivatives of drift into coefficients of
    %Z-series.
    [b0,b] = CalculateDriftbCoeffs04(wMu0dt,dwMu0dtdw,a);
%    b0
    b(7: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] = EvolveSDESeriesSolution04(a0,a,b0,b,sigma0,dt);
    c(7:10)=0.0;
    w0=c0;
    
%     tt
%     a0
%     a
%     b0
%     b
%     c0
%     c
%     tt
    %str=input('Look at solution');
    
    end

    %The newly calculated solution is assigned to previous solution
    %coefficients so that we could calculate further future solution with 
    %these coefficients.
    a0=c0;
    a(1:6)=c(1:6);
    
    

end    
wnStart=1;
    w(wnStart:Nn)=c0+c(1)*Z(wnStart:Nn)+c(2)*Z(wnStart:Nn).^2+c(3)*Z(wnStart:Nn).^3+ ...
        c(4)*Z(wnStart:Nn).^4+c(5)*Z(wnStart:Nn).^5+c(6)*Z(wnStart:Nn).^6+ ...
        c(7)*Z(wnStart:Nn).^7+c(8)*Z(wnStart:Nn).^8+c(9)*Z(wnStart:Nn).^9+ ...
        c(10)*Z(wnStart:Nn).^10;
    


%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(wnStart:Nn) =(c(1)+2*c(2)*Z(wnStart:Nn)+3*c(3)*Z(wnStart:Nn).^2+ ...
%         4*c(4)*Z(wnStart:Nn).^3+5*c(5)*Z(wnStart:Nn).^4+6*c(6)*Z(wnStart:Nn).^5+ ...
%         7*c(7)*Z(wnStart:Nn).^6+8*c(8)*Z(wnStart:Nn).^7+9*c(9)*Z(wnStart:Nn).^8+ ...
%         10*c(10)*Z(wnStart:Nn).^9)*dNn;

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
yy0


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(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
.
.
function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)

Fp2=Fp1(1,1,2,1)/(1-gamma);

%ExpnOrder=10;
NoOfTerms=19;
NoOfTerms=9;%excluding dt^3 terms

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;


YqCoeff(1)=YqCoeff0(1,1,2,1)*dt;
YqCoeff(2)=YqCoeff0(1,2,1,1)*dt;
YqCoeff(3)=YqCoeff0(2,1,1,1)*dt;

YqCoeff(4)=YqCoeff0(1,1,3,1)*dt^2;
YqCoeff(5)=YqCoeff0(1,2,2,1)*dt^2;
YqCoeff(6)=YqCoeff0(2,1,2,1)*dt^2;
YqCoeff(7)=YqCoeff0(1,3,1,1)*dt^2;
YqCoeff(8)=YqCoeff0(2,2,1,1)*dt^2;
YqCoeff(9)=YqCoeff0(3,1,1,1)*dt^2;

YqCoeff(10)=YqCoeff0(1,1,4,1)*dt^3;
YqCoeff(11)=YqCoeff0(1,2,3,1)*dt^3;
YqCoeff(12)=YqCoeff0(2,1,3,1)*dt^3;
YqCoeff(13)=YqCoeff0(1,3,2,1)*dt^3;
YqCoeff(14)=YqCoeff0(2,2,2,1)*dt^3;
YqCoeff(15)=YqCoeff0(3,1,2,1)*dt^3;
YqCoeff(16)=YqCoeff0(1,4,1,1)*dt^3;
YqCoeff(17)=YqCoeff0(2,3,1,1)*dt^3;
YqCoeff(18)=YqCoeff0(3,2,1,1)*dt^3;
YqCoeff(19)=YqCoeff0(4,1,1,1)*dt^3;

Fp(1)=Fp1(1,1,2,1);
Fp(2)=Fp1(1,2,1,1);
Fp(3)=Fp1(2,1,1,1);

Fp(4)=Fp1(1,1,3,1);
Fp(5)=Fp1(1,2,2,1);
Fp(6)=Fp1(2,1,2,1);
Fp(7)=Fp1(1,3,1,1);
Fp(8)=Fp1(2,2,1,1);
Fp(9)=Fp1(3,1,1,1);

Fp(10)=Fp1(1,1,4,1);
Fp(11)=Fp1(1,2,3,1);
Fp(12)=Fp1(2,1,3,1);
Fp(13)=Fp1(1,3,2,1);
Fp(14)=Fp1(2,2,2,1);
Fp(15)=Fp1(3,1,2,1);
Fp(16)=Fp1(1,4,1,1);
Fp(17)=Fp1(2,3,1,1);
Fp(18)=Fp1(3,2,1,1);
Fp(19)=Fp1(4,1,1,1);

Fp2(1:19)=Fp(1:19)/(1-gamma);


wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)*(1-gamma)/((1-gamma)*w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp2(mm)-(nn-1))*(1-gamma)/((1-gamma)*w0);
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
        
    
end



end

.
.
function [b0,b] = CalculateDriftbCoeffs04(wMu0dt,dwMu0dtdw,a)



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


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


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


b(7:10)=0.0;

% 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));
% 
% 
% 
% 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));
% 
% 
% 
% 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 [c0,c] = EvolveSDESeriesSolution04(a0,a,b0,b,sigma0,dt)

c0=0.0;
c(1:10)=0.0;

c0=1/a(1)^2*(a0*a(1)^2+b0*a(1)^2+ ...
    (.5*sigma0^2*dt*0)+ ...
    (.5*sigma0^2*dt*2*a(2)));


% c(1)=1/a(1).^2*((+a(1).^3  + 4* a0* a(1)* a(2) ) + ...
%     (+4 *a(1)* a(2)* b0  + a(1).^2 *b(1) )- ...
%     (+4 *a(1)* a(2) *c0 )+ ...
%     (.5*sigma0^2*dt*a(1))+ ...
%     (.5*sigma0^2*dt*6*a(3)));


A1=(+a(1).^3  + 4* a0* a(1)* a(2));
B1=+4 *a(1)* a(2)* b0  + a(1).^2 *b(1);
C1=+4 *a(1)* a(2) *c0;
D1=.5*sigma0^2*dt*a(1);
%D1=sigma0*sqrt(dt)*a(1)^2;
E1=.5*sigma0^2*dt*6*a(3);


c(1)=1/a(1).^2*(A1+B1-C1+D1+E1);
%c(1)=1/a(1).^2*(A1-C1+sqrt((A1)^2+D1^2)+E1);



c(2)=1/a(1)^2*( (+5 *a(1)^2* a(2)  + 4* a0* a(2)^2  + 6* a0* a(1)* a(3) )+ ...
    +(+4* a(2)^2* b0  + 6 *a(1)* a(3)* b0  + 4* a(1)* a(2)* b(1)  + a(1)^2* b(2) )- ...
    (+4* a(2)^2* c0 +6* a(1)* a(3)* c0 +4* a(1)* a(2)* c(1) )+ ...
    (.5*sigma0^2*dt*2*a(2))+ ...
    (.5*sigma0^2*dt*12*a(4)));
c(3)=1/a(1)^2*((+8* a(1)* a(2)^2  + 7* a(1)^2* a(3)  + 12* a0* a(2)* a(3)  + 8* a0* a(1)* a(4) )+...
    (+12* a(2)* a(3)* b0  + 8* a(1)* a(4)* b0  + 4* a(2)^2* b(1)  + 6* a(1)* a(3)* b(1)  + ...
    4* a(1)* a(2)* b(2)  + a(1)^2* b(3)  )-...
    (+12* a(2)* a(3)* c0 +8* a(1)* a(4)* c0 +4* a(2)^2* c(1) +6* a(1)* a(3)* c(1) +4* a(1)* a(2)* c(2) )+ ...
    (.5*sigma0^2*dt*3* a(3) )+ ...
    (.5*sigma0^2*dt*20* a(5) ));

c(4)=1/a(1)^2*((+4* a(2)^3  + 22* a(1)* a(2)* a(3)  + 9* a0* a(3)^2  + 9* a(1)^2* a(4)  + ... 
    16* a0* a(2)* a(4)  + 10* a0* a(1)* a(5) )+...
    (+9* a(3)^2* b0  + 16* a(2)* a(4)* b0  + 10* a(1)* a(5)* b0  + ...
    12* a(2)* a(3)* b(1)  + 8* a(1)* a(4)* b(1)  + 4* a(2)^2* b(2)  + 6* a(1)* a(3)* b(2)  + ...
    4* a(1)* a(2)* b(3)  + a(1)^2* b(4) )-...
    (+9* a(3)^2* c0 +16* a(2)* a(4)* c0 +10 *a(1)* a(5)* c0 +12* a(2)* a(3)* c(1) +8* a(1)* a(4)* c(1) +4* a(2)^2* c(2) +6* a(1)* a(3)* c(2) + ...
    4* a(1)* a(2)* c(3) )+ ...
    (.5*sigma0^2*dt*4* a(4) )+ ...
    (.5*sigma0^2*dt*30* a(6) ));
    

c(5)=1/a(1)^2*((+16* a(2)^2* a(3)  + 15* a(1)* a(3)^2  + 28* a(1)* a(2)* a(4)  + ...
    24* a0* a(3)* a(4)  + 11* a(1)^2* a(5)  + 20* a0* a(2)* a(5)  + 12* a0* a(1)* a(6) )+ ...
    (+24* a(3)* a(4)* b0  + 20* a(2)* a(5)* b0  + 12* a(1)* a(6)* b0  + ...
    9* a(3)^2* b(1)  + 16* a(2)* a(4)* b(1)  + 10* a(1)* a(5)* b(1)  + ...
    12* a(2)* a(3)* b(2)  + 8* a(1)* a(4)* b(2)  + 4* a(2)^2* b(3)  + 6* a(1)* a(3)* b(3)  + ...
    4* a(1)* a(2)* b(4)  + a(1)^2* b(5) )-...
    (+24* a(3)* a(4)* c0 +20* a(2)* a(5)* c0 +12* a(1)* a(6)* c0 +9* a(3)^2* c(1) +16* a(2)* a(4)* c(1) +10* a(1)* a(5)* c(1) +12* a(2)* a(3)* c(2) + ...
    8* a(1)* a(4)* c(2) +4* a(2)^2* c(3) +6* a(1)* a(3)* c(3) +4* a(1)* a(2)* c(4) )+...
    (.5*sigma0^2*dt*5* a(5) )+ ...
    (.5*sigma0^2*dt*42* a(7) ));

c(6)=1/a(1)^2*((+21* a(2)* a(3)^2  + 20* a(2)^2* a(4)  + 38* a(1)* a(3)* a(4)  + 16* a0* a(4)^2  + ...
    34* a(1)* a(2)* a(5)  + 30* a0* a(3)* a(5)  + 13* a(1)^2* a(6)  + ...
    24* a0* a(2)* a(6)  + 14* a0* a(1)* a(7) )+...
    (+16* a(4)^2* b0  + 30* a(3)* a(5)* b0  + 24* a(2)* a(6)* b0  + ...
    14* a(1)* a(7)* b0  + 24* a(3)* a(4)* b(1)  + 20* a(2)* a(5)* b(1)  + ...
    12* a(1)* a(6)* b(1)  + 9* a(3)^2* b(2)  + 16* a(2)* a(4)* b(2)  + ...
    10* a(1)* a(5)* b(2)  + 12* a(2)* a(3)* b(3)  + 8* a(1)* a(4)* b(3)  + 4* a(2)^2* b(4)  + ...
    6* a(1)* a(3)* b(4)  + 4* a(1)* a(2)* b(5)  + a(1)^2* b(6) )- ...
    (+16* a(4)^2* c0 +30* a(3)* a(5)* c0 +24* a(2)* a(6)* c0 +14* a(1)* a(7)* c0 +24 *a(3)* a(4)* c(1) +20* a(2)* a(5)* c(1) +12* a(1)* a(6)* c(1) + ...
    9* a(3)^2* c(2) +16 *a(2)* a(4)* c(2) +10* a(1)* a(5)* c(2) +12* a(2)* a(3)* c(3) +8* a(1)* a(4)* c(3) +4* a(2)^2* c(4) +6* a(1)* a(3)* c(4) + ...
    4* a(1)* a(2)* c(5) )+ ...
    (.5*sigma0^2*dt*6* a(6) )+ ...
    (.5*sigma0^2*dt*56* a(8) ));

c(7:10)=0.0;    

% c(7)=1/a(1)^2*((+9* a(3)^3  + 52* a(2)* a(3)* a(4)  + 24* a(1)* a(4)^2  + 24* a(2)^2* a(5)  +... 
%    46* a(1)* a(3)* a(5)  + 40* a0* a(4)* a(5)  + 40* a(1)* a(2)* a(6)  + ...
%    36* a0* a(3)* a(6)  + 15* a(1)^2* a(7)  + 28* a0* a(2)* a(7)  + 16* a0* a(1)* a(8) )+...
%    (+40* a(4)* a(5)* b0  + 36* a(3)* a(6)* b0  + 28* a(2)* a(7)* b0  + ...
%    16* a(1)* a(8)* b0  + 16* a(4)^2* b(1)  + 30* a(3)* a(5)* b(1)  + ...
%    24* a(2)* a(6)* b(1)  + 14* a(1)* a(7)* b(1)  + 24* a(3)* a(4)* b(2)  + ...
%    20* a(2)* a(5)* b(2)  + 12* a(1)* a(6)* b(2)  + 9* a(3)^2* b(3)  + ...
%    16* a(2)* a(4)* b(3)  + 10* a(1)* a(5)* b(3)  + 12* a(2)* a(3)* b(4)  + ...
%    8* a(1)* a(4)* b(4)  + 4* a(2)^2* b(5)  + 6* a(1)* a(3)* b(5)  + 4* a(1)* a(2)* b(6)  + ... 
%    a(1)^2* b(7) )-...     
%    (+40* a(4)* a(5)* c0 +36* a(3)* a(6)* c0 +28* a(2)* a(7)* c0 +16* a(1)* a(8)* c0 +16* a(4)^2* c(1) +30* a(3)* a(5)* c(1) + ...
%     24* a(2)* a(6)* c(1) +14* a(1)* a(7)* c(1) +24* a(3)* a(4)* c(2) +20* a(2)* a(5)* c(2) +12* a(1)* a(6)* c(2) +9* a(3)^2* c(3) + ...
%     16* a(2)* a(4)* c(3) +10* a(1)* a(5)* c(3) +12* a(2)* a(3)* c(4) +8* a(1)* a(4)* c(4) +4* a(2)^2* c(5) +6* a(1)* a(3)* c(5) + ...
%     4* a(1)* a(2)* c(6) )+ ...
%     (.5*sigma0^2*dt*7* a(7) )+ ...
%     (.5*sigma0^2*dt* 72* a(9) ));
%     
% c(8)=1/a(1)^2*((33* a(3)^2* a(4)  + 32* a(2)* a(4)^2  + 62* a(2)* a(3)* a(5)  + 58 *a(1) *a(4)* a(5)  +... 
%    25* a0* a(5)^2  + 28* a(2)^2* a(6)  + 54 *a(1)* a(3)* a(6)  + ...
%    48* a0* a(4)* a(6)  + 46* a(1)* a(2)* a(7)  + 42* a0* a(3)* a(7)  + ...
%    17* a(1)^2* a(8)  + 32* a0* a(2)* a(8)  + 18* a0* a(1)* a(9) )+ ...
%    (+25* a(5)^2* b0  + 48* a(4)* a(6)* b0  + 42* a(3)* a(7)* b0  + ...
%    32* a(2)* a(8)* b0  + 18* a(1)* a(9)* b0  + 40* a(4)* a(5)* b(1)  + ...
%    36* a(3)* a(6)* b(1)  + 28* a(2)* a(7)* b(1)  + 16* a(1)* a(8)* b(1)  + ...
%    16* a(4)^2* b(2)  + 30* a(3)* a(5)* b(2)  + 24* a(2)* a(6)* b(2)  + ...
%    14* a(1)* a(7) *b(2)  + 24* a(3)* a(4)* b(3)  + 20* a(2)* a(5)* b(3)  + ...
%    12* a(1)* a(6)* b(3)  + 9* a(3)^2* b(4)  + 16* a(2)* a(4)* b(4)  + ...
%    10* a(1)* a(5)* b(4)  + 12* a(2)* a(3)* b(5)  + 8* a(1)* a(4)* b(5)  + 4* a(2)^2* b(6)  + ...
%    6* a(1)* a(3)* b(6)  + 4* a(1)* a(2)* b(7)  + a(1)^2* b(8) )-...
%    (+25* a(5)^2* c0 +48* a(4)* a(6)* c0 +42* a(3)* a(7)* c0 +32* a(2)* a(8)* c0 +18* a(1)* a(9)* c0 +40* a(4)* a(5)* c(1) +36* a(3)* a(6)* c(1) + ...
%     28* a(2)* a(7)* c(1) +16* a(1)* a(8)* c(1) +16* a(4)^2* c(2) +30* a(3)* a(5)* c(2) +24* a(2)* a(6)* c(2) +14* a(1)* a(7)* c(2) +24* a(3)* a(4)* c(3) +...
%     20* a(2)* a(5)* c(3) +12* a(1)* a(6)* c(3) +9* a(3)^2* c(4) +16* a(2)* a(4)* c(4) +10* a(1)* a(5)* c(4) +12* a(2)* a(3)* c(5) +8* a(1)* a(4)* c(5) + ...
%     4* a(2)^2* c(6) +6* a(1)* a(3)* c(6) +4* a(1)* a(2)* c(7) )+...
%     (.5*sigma0^2*dt*8* a(8) )+ ...
%     (.5*sigma0^2*dt* 90* a(10) ));
%     
% 
% c(9)=1/a(1)^2*((+20* a0* a(1)* a(10)  + 40* a(3)* a(4)^2  + 39* a(3)^2* a(5)  + ... 
%    76* a(2)* a(4)* a(5)  + 35* a(1)* a(5)^2  + 72* a(2)* a(3)* a(6)  + ...
%    68* a(1)* a(4)* a(6)  + 60* a0* a(5)* a(6)  + 32* a(2)^2* a(7)  + ...
%    62* a(1)* a(3)* a(7)  + 56* a0* a(4)* a(7)  + 52* a(1)* a(2)* a(8)  + ...
%    48* a0* a(3)* a(8)  + 19* a(1)^2* a(9)  + 36* a0* a(2)* a(9) )+ ...
%    (+20* a(1)* a(10)* b0  + 60* a(5)* a(6)* b0  + 56* a(4)* a(7)* b0  + ...
%     48* a(3)* a(8)* b0  + 36* a(2)* a(9)* b0  + 25* a(5)^2* b(1)  + ...
%     48* a(4)* a(6)* b(1)  + 42* a(3)* a(7)* b(1)  + 32* a(2)* a(8)* b(1)  + ...
%     18* a(1)* a(9)* b(1)  + 40* a(4)* a(5) *b(2)  + 36* a(3)* a(6)* b(2)  + ...
%     28* a(2)* a(7)* b(2)  + 16* a(1)* a(8)* b(2)  + 16* a(4)^2* b(3)  + ...
%     30* a(3)* a(5)* b(3)  + 24* a(2)* a(6)* b(3)  + 14* a(1)* a(7)* b(3)  + ...
%     24* a(3)* a(4)* b(4)  + 20* a(2)* a(5)* b(4)  + 12* a(1)* a(6)* b(4)  + ...
%     9* a(3)^2* b(5)  + 16* a(2)* a(4)* b(5)  + 10* a(1)* a(5)* b(5)  + ...
%     12* a(2)* a(3)* b(6)  + 8* a(1)* a(4)* b(6)  + 4* a(2)^2* b(7)  + 6* a(1)* a(3)* b(7)  + ...
%     4* a(1)* a(2)* b(8)  + a(1)^2* b(9) )-...
%     (+20* a(1)* a(10)* c0 +60* a(5)* a(6)* c0 +56* a(4)* a(7)* c0 +48* a(3)* a(8)* c0 +36* a(2)* a(9)* c0 +25* a(5)^2* c(1) +48* a(4)* a(6)* c(1) + ...
%     42* a(3)* a(7)* c(1) +32* a(2)* a(8)* c(1) +18* a(1)* a(9)* c(1) +40* a(4)* a(5)* c(2) +36* a(3)* a(6)* c(2) +28* a(2)* a(7)* c(2) +16* a(1)* a(8)* c(2) + ...
%     16* a(4)^2* c(3) +30* a(3)* a(5)* c(3) +24* a(2)* a(6)* c(3) +148* a(1)* a(7)* c(3) +24* a(3)* a(4)* c(4) +20* a(2)* a(5)* c(4) +12* a(1)* a(6)* c(4) + ...
%     9* a(3)^2* c(5) +16* a(2)* a(4)* c(5) +10* a(1)* a(5)* c(5) +12* a(2)* a(3)* c(6) +8* a(1)* a(4)* c(6) +4* a(2)^2* c(7) +6* a(1)* a(3)* c(7) + ...
%     4* a(1)* a(2)* c(8) )+ ...
%     (.5*sigma0^2*dt*9* a(9) )+ ...
%     (.5*sigma0^2*dt* 0));
% 
% c(10)=1/a(1)^2*((+21* a(1)^2* a(10)  + 40* a0* a(10)* a(2)  + 16* a(4)^3  + ... 
%     94* a(3)* a(4)* a(5)  + 45* a(2)* a(5)^2  + 45* a(3)^2* a(6)  + ...
%     88* a(2)* a(4)* a(6)  + 82* a(1)* a(5)* a(6)  + 36* a0* a(6)^2  + ...
%     82* a(2)* a(3)* a(7)  + 78* a(1)* a(4)* a(7)  + 70* a0* a(5)* a(7)  + ...
%     36* a(2)^2* a(8)  + 70* a(1)* a(3)* a(8)  + 64* a0* a(4)* a(8)  + ...
%     58* a(1)* a(2)* a(9)  + 54* a0* a(3)* a(9) )+ ...     
%     (+40* a(10)* a(2)* b0  + 36* a(6)^2* b0  + 70* a(5)* a(7)* b0  + ...
%     64* a(4)* a(8)* b0  + 54* a(3)* a(9)* b0  + 20* a(1)* a(10)* b(1)  + ...
%     60* a(5)* a(6)* b(1)  + 56* a(4)* a(7)* b(1)  + 48* a(3)* a(8)* b(1)  + ...
%     36* a(2)* a(9)* b(1)  + a(1)^2* b(10)  + 25* a(5)^2* b(2)  + ...
%     48* a(4)* a(6)* b(2)  + 42* a(3)* a(7)* b(2)  + 32* a(2)* a(8)* b(2)  +... 
%     18* a(1)* a(9)* b(2)  + 40* a(4)* a(5)* b(3)  + 36* a(3)* a(6)* b(3)  + ...
%     28* a(2)* a(7)* b(3)  + 16* a(1)* a(8)* b(3)  + 16* a(4)^2* b(4)  + ...
%     30* a(3)* a(5)* b(4)  + 24* a(2)* a(6)* b(4)  + 14* a(1)* a(7)* b(4)  + ...
%     24* a(3)* a(4)* b(5)  + 20* a(2)* a(5)* b(5)  + 12* a(1)* a(6)* b(5)  + ...
%     9* a(3)^2* b(6)  + 16* a(2)* a(4)* b(6)  + 10 *a(1) *a(5) *b(6)  + ...
%     12 *a(2) *a(3) *b(7)  + 8 *a(1) *a(4) *b(7)  + 4 *a(2)^2 *b(8)  + ...
%     6 *a(1)* a(3)* b(8)  + 4* a(1)* a(2)* b(9) )-...
%     (+40* a(10)* a(2)* c0 +36* a(6)^2* c0 +70* a(5)* a(7)* c0 +64* a(4)* a(8)* c0 +54* a(3)* a(9)* c0 +20* a(1)* a(10)* c(1) +60* a(5)* a(6)* c(1) + ...
%     56* a(4)* a(7)* c(1) +48* a(3)* a(8)* c(1) +36* a(2)* a(9)* c(1) +25* a(5)^2* c(2) +48* a(4)* a(6)* c(2) +42* a(3)* a(7)* c(2) + ...
%     32* a(2)* a(8)* c(2) +18* a(1)* a(9)* c(2) +40* a(4)* a(5)* c(3) +36* a(3)* a(6)* c(3) +28* a(2)* a(7)* c(3) +16* a(1)* a(8)* c(3) +16* a(4)^2* c(4) + ...
%     30* a(3)* a(5)* c(4) +24* a(2)* a(6)* c(4) +14* a(1)* a(7)* c(4) +24* a(3)* a(4)* c(5) +20* a(2)* a(5)* c(5) +12* a(1)* a(6)* c(5) +9* a(3)^2* c(6) + ...
%     16* a(2)* a(4)* c(6) +10* a(1)* a(5)* c(6) +12* a(2)* a(3)* c(7) +8* a(1)* a(4)* c(7) +4* a(2)^2* c(8) +6* a(1)* a(3)* c(8) +4* a(1)* a(2)* c(9) )+ ...
%     (.5*sigma0^2*dt*10*a(10) )+ ...
%     (.5*sigma0^2*dt* 0 ));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55


end

.
.
.
function [wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)

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

Fp2=Fp1/(1-gamma);

wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
    YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
    YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
    (YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
    YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
    YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
    YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
    YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
    YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2 + ...
     (YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
     YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
     YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
     YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
     YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
     YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
     YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
     YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
     YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
     YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3+ ...
     (YqCoeff0(1,1,5,1).*yy(wnStart:Nn).^Fp1(1,1,5,1)+ ...
      YqCoeff0(1,2,4,1).*yy(wnStart:Nn).^Fp1(1,2,4,1)+ ...
      YqCoeff0(2,1,4,1).*yy(wnStart:Nn).^Fp1(2,1,4,1)+ ...
      YqCoeff0(1,3,3,1).*yy(wnStart:Nn).^Fp1(1,3,3,1)+ ...
       YqCoeff0(2,2,3,1).*yy(wnStart:Nn).^Fp1(2,2,3,1)+ ...
       YqCoeff0(3,1,3,1).*yy(wnStart:Nn).^Fp1(3,1,3,1)+ ...
       YqCoeff0(1,4,2,1).*yy(wnStart:Nn).^Fp1(1,4,2,1)+ ...
       YqCoeff0(2,3,2,1).*yy(wnStart:Nn).^Fp1(2,3,2,1)+ ...
       YqCoeff0(3,2,2,1).*yy(wnStart:Nn).^Fp1(3,2,2,1)+ ...
       YqCoeff0(4,1,2,1).*yy(wnStart:Nn).^Fp1(4,1,2,1)+ ...
       YqCoeff0(1,5,1,1).*yy(wnStart:Nn).^Fp1(1,5,1,1)+ ...
       YqCoeff0(2,4,1,1).*yy(wnStart:Nn).^Fp1(2,4,1,1)+ ...
       YqCoeff0(3,3,1,1).*yy(wnStart:Nn).^Fp1(3,3,1,1)+  ...
       YqCoeff0(4,2,1,1).*yy(wnStart:Nn).^Fp1(4,2,1,1)+ ...
       YqCoeff0(5,1,1,1).*yy(wnStart:Nn).^Fp1(5,1,1,1))*dt^4;

   
   dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,1,2,1)-1)+ ...
    YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,2,1,1)-1)+ ...
    YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(2,1,1,1)-1))*dt ;
   
   
%wMu0dtOdt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)*0+ ...
%    YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
%    YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt ;

c1(wnStart:Nn)=((YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt))+ ...
     (YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
     YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
     (YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
     YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
     YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5+ ...
     (YqCoeff0(1,1,4,2).*yy(wnStart:Nn).^Fp1(1,1,4,2)+YqCoeff0(1,2,3,2).*yy(wnStart:Nn).^Fp1(1,2,3,2)+ ...
     YqCoeff0(2,1,3,2).*yy(wnStart:Nn).^Fp1(2,1,3,2)+YqCoeff0(1,3,2,2).*yy(wnStart:Nn).^Fp1(1,3,2,2)+ ...
     YqCoeff0(2,2,2,2).*yy(wnStart:Nn).^Fp1(2,2,2,2)+ YqCoeff0(3,1,2,2).*yy(wnStart:Nn).^Fp1(3,1,2,2)+ ...
     YqCoeff0(1,4,1,2).*yy(wnStart:Nn).^Fp1(1,4,1,2)+YqCoeff0(2,3,1,2).*yy(wnStart:Nn).^Fp1(2,3,1,2)+ ...
     YqCoeff0(3,2,1,2).*yy(wnStart:Nn).^Fp1(3,2,1,2)+YqCoeff0(4,1,1,2).*yy(wnStart:Nn).^Fp1(4,1,1,2)).*dt^3.5);



c2(wnStart:Nn)=(YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
     YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2;%+ ...
     (YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
     YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
     YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;
 
 
% (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+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths)
 
 
 
end
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 22nd, 2021, 1:39 pm

Friends, I made some improvements in the program. I just made the drift second order accurate while work still remains to be done to add further terms that I will do in next few days. But the program was working great out to five years so I decided to distribute a new interim version.
For volatility exponents close to .5, some condition has to be applied how fast the coefficients can be allowed to increase. I hope it will increase substantially when I increase order. But the new program is still very good.
Here are some graphs for densities of SDEs out to five years. I am going to post the exact program used to create these graphs in a minute in next post.
.
.
Image

Image

Image

Image

Image

Image

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

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

November 22nd, 2021, 1:44 pm

Here is the new interim program for recursive calculation of series solution to densities of SDEs.
.
.
function [] = FPESeriesCoeffsSolutionWilmott()

%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*5;%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=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46;  % 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=1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=1.0;%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)-3.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);

%ZProbT(1)=normcdf(.5*ZT(1)+.5*ZT(2),0,1)-normcdf(.5*ZT(1)+.5*ZT(2)-dNn,0,1);
%ZProbT(NnT)=normcdf(.5*ZT(NnT)+.5*ZT(NnT-1)+dNn,0,1)-normcdf(.5*ZT(NnT)+.5*ZT(NnT-1),0,1);
%ZProbT(2:NnT-1)=normcdf(.5*ZT(2:NnT-1)+.5*ZT(3:NnT),0,1)-normcdf(.5*ZT(2:NnT-1)+.5*ZT(1:NnT-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

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

ExpnOrder=10;
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)

    %The function below calculates drift and its derivatives at w(Z=0);    
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
    [wMu0dt,dwMu0dtdw] =BesselDriftAndDerivatives02(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder);
    
    a(7:10)=0.0;
    %The function below converts derivatives of drift into coefficients of
    %Z-series.
    [b0,b] = CalculateDriftbCoeffs04(wMu0dt,dwMu0dtdw,a);
%    b0
    b(7: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] = EvolveSDESeriesSolution04(a0,a,b0,b,sigma0,dt);
    c(7:10)=0.0;
    if(gamma<=.65)
    for mm=1:5
        if((sign(c(mm))~=sign(c(mm+1)))&&(abs(c(mm))<abs(c(mm+1)*5)))
        %if((abs(c(mm))<abs(c(mm+1)*8)))    
            %c(mm+1)=c(mm+1)*abs(c(mm)/c(mm+1)/11);
            %c(mm+1)=sign(c(mm+1))*abs(c(mm)/4.5);
            c(mm+1)=sign(c(mm+1))*abs(c(mm+1)/5.5);
        end
    end
    end
    
% for mm=1:4
%         if((abs(c(mm))<abs(c(mm+2)*16)))
%             %c(mm+1)=c(mm+1)*abs(c(mm)/c(mm+1)/11);
%             %c(mm+1)=sign(c(mm+1))*abs(c(mm)/4.5);
%             c(mm+2)=sign(c(mm+2))*abs(c(mm)/16);
%         end
%     end

    
    w0=c0;
    
%     tt
%     a0
%     a
%     b0
%     b
%     c0
%     c
%     tt
    %str=input('Look at solution');
    
    end

    %The newly calculated solution is assigned to previous solution
    %coefficients so that we could calculate further future solution with 
    %these coefficients.
    a0=c0;
    a(1:6)=c(1:6);
    
    

end    
wnStart=1;
    w(wnStart:Nn)=c0+c(1)*Z(wnStart:Nn)+c(2)*Z(wnStart:Nn).^2+c(3)*Z(wnStart:Nn).^3+ ...
        c(4)*Z(wnStart:Nn).^4+c(5)*Z(wnStart:Nn).^5+c(6)*Z(wnStart:Nn).^6+ ...
        c(7)*Z(wnStart:Nn).^7+c(8)*Z(wnStart:Nn).^8+c(9)*Z(wnStart:Nn).^9+ ...
        c(10)*Z(wnStart:Nn).^10;
 c   
for nn=1:Nn-1
    if(w(nn)>w(nn+1))
        w(nn)=0;
        wnStart=nn+1;
    end
end

%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(wnStart:Nn) =(c(1)+2*c(2)*Z(wnStart:Nn)+3*c(3)*Z(wnStart:Nn).^2+ ...
%         4*c(4)*Z(wnStart:Nn).^3+5*c(5)*Z(wnStart:Nn).^4+6*c(6)*Z(wnStart:Nn).^5+ ...
%         7*c(7)*Z(wnStart:Nn).^6+8*c(8)*Z(wnStart:Nn).^7+9*c(9)*Z(wnStart:Nn).^8+ ...
%         10*c(10)*Z(wnStart:Nn).^9)*dNn;

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
yy0


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(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
.
.
The function to calculate drift and its derivatives has changed and I am posting it again here. Rest of the functions are old.
.
function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=9;%excluding dt^3 terms

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;


YqCoeff(1)=mu1*dt;
YqCoeff(2)=mu2*dt;
YqCoeff(3)=-.5*sigma0^2*gamma*dt;

YqCoeff(4)=mu1^2*(beta1-gamma)*dt^2;
YqCoeff(5)=mu2^2*(beta2-gamma)*dt^2;
YqCoeff(6)=-.25*sigma0^4*gamma^2*(1-gamma)*dt^2;
YqCoeff(7)=(mu1*mu2*(beta1-gamma)+mu1*mu2*(beta2-gamma))*dt^2;
YqCoeff(8)=.5*mu1*sigma0^2*gamma*(1-gamma)*dt^2;
YqCoeff(9)=.5*mu2*sigma0^2*gamma*(1-gamma)*dt^2;


% YqCoeff(10)=YqCoeff0(1,1,4,1)*dt^3;
% YqCoeff(11)=YqCoeff0(1,2,3,1)*dt^3;
% YqCoeff(12)=YqCoeff0(2,1,3,1)*dt^3;
% YqCoeff(13)=YqCoeff0(1,3,2,1)*dt^3;
% YqCoeff(14)=YqCoeff0(2,2,2,1)*dt^3;
% YqCoeff(15)=YqCoeff0(3,1,2,1)*dt^3;
% YqCoeff(16)=YqCoeff0(1,4,1,1)*dt^3;
% YqCoeff(17)=YqCoeff0(2,3,1,1)*dt^3;
% YqCoeff(18)=YqCoeff0(3,2,1,1)*dt^3;
% YqCoeff(19)=YqCoeff0(4,1,1,1)*dt^3;

Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;

Fp(4)=2*(beta1-gamma)/(1-gamma)-1;
Fp(5)=2*(beta2-gamma)/(1-gamma)-2;
Fp(6)=-3;
Fp(7)=(beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-1;
Fp(8)=(beta1-gamma)/(1-gamma)-2;
Fp(9)=(beta2-gamma)/(1-gamma)-2;

% Fp(10)=Fp1(1,1,4,1);
% Fp(11)=Fp1(1,2,3,1);
% Fp(12)=Fp1(2,1,3,1);
% Fp(13)=Fp1(1,3,2,1);
% Fp(14)=Fp1(2,2,2,1);
% Fp(15)=Fp1(3,1,2,1);
% Fp(16)=Fp1(1,4,1,1);
% Fp(17)=Fp1(2,3,1,1);
% Fp(18)=Fp1(3,2,1,1);
% Fp(19)=Fp1(4,1,1,1);

Fp2(1:9)=Fp(1:9);%/(1-gamma);


wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)*(1-gamma)/((1-gamma)*w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp2(mm)-(nn-1))*(1-gamma)/((1-gamma)*w0);
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
        
    
end



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

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

November 23rd, 2021, 3:24 am

In my morning walk today, I was thinking of new possibilities that open up with a series(in Z) representation of SDEs. Though I am sure experts in stochastics can tell about many more possibilities, here is what I could think apart from obvious possibilities.

1. Some work would have to be done but we could tell about transition probabilities between evolution of an SDE between two points in time when we have their series representation and values of Z at both points is given. I will come back to it in a few days with examples.

2. A series representation has ramifications for calculation of pricing with a numeraire and measure change.
Suppose an asset 
[$]X(Z)=a_1 \, Z+ a_2 \, Z^2 + a_3 \, Z^3+ ...[$]
while numeraire has a representation (with common Z in case of one factor model)
[$]\eta(Z)=b_1 \, Z+ b_2 \, Z^2 + b_3 \, Z^3+ ...[$]

Then numeraire discounted distribution could be given as
[$]\frac{X(Z)}{\eta(Z)} \,P(Z) [$]  where [$]P(Z)[$] is density of standard normal.
which could be integrated to get the point estimate of expectation and variances of the numeraire discounted payoffs etc.

measure change could be done easily and its distribution could be found as 
[$]\frac{X(Z)}{\eta_1(Z)} \, \frac{\eta_1(Z)}{\eta_2(Z)} \,P(Z) [$] 


3. There would be new possibilities in calculation of path integrals and other stochastic integrals of the diffusions. Whenever I tried this earlier, I was hampered since SDE measure for path integrals would change from one time to other time. It would be more interesting to remove the layer of complexity associated with probability density and directly represent the evolution of path integrals and stochastic integrals in terms of a series in Z. I hope that this would help in easy calculation of path integrals and other related stochastic integrals in time.

4. Lot of intelligent friends are working with monte carlo simulations with Artificial intelligence. I think it will be far easier to be able to represent the time evolution of SDE variables in terms of a series in Z and try to find the time dependent coefficients of this series with AI. If the coefficient representation is good and AI works well, it will give a nice smooth distribution of SDE variable. I think this new approach would be far easier as compared to approximating the full blown distribution of SDE with AI.

But a lot of work remains to be done how to calculate the coefficients correctly and it would be more difficult in multiple dimensions. And research has to be done on relative structure and time growth of coefficients and their proper bounds for valid densities etc.

I think I will write and distribute a new program that takes an SDE distribution and finds coefficients of its Z-series using optimization so that a Z-series could be found when needed. Such a function would be very useful. 
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 24th, 2021, 9:11 am

Before writing this post, I want to apologize to good jews but want to tell them there are enough evil jews who are doing every thing to tarnish the great goodwill jews have so rightly earned through generations of hard work and being nice and kind to other people in the society. I want to request good jews to Please force these evil(jewish) people to end the bad practices that are starting to take root in American and other societies on their behest and before it is too late. I also really hope that my mind control would finally end(as I have been begging for past twenty years) and I would never have to write a post like this again.

Friends I have been writing on this forum since its very starting years. And a lot of people people know me from very early years also because I had another web site LIBOR Market Model web that hosted almost all the papers (more than 40-50) about the then cutting edge LIBOR market model and people freely downloaded papers and code from my web site. I continued to mention about my mind control persecution from time to time but only started writing regularly about it when I arrived in UK on a Highly Skilled Migrant work visa. Mind control agencies did everything in their power to fail me in UK. People were asked to not give me a job. People I contacted were routinely told brazen lies about me totally shamelessly. When I realized that I had no future in UK, I returned to Pakistan after a six month stint in UK. Since my hard earned money had thoroughly insulted mind control infrastructure in UK as they had to drug entire cities with mind control chemicals in UK (read briefly about it in my application to HRW below) for a nobody like me and American government had to move all their grand and influential puppets in UK (again for a nobody like me) to force people to drug the manufactured food. British Banks went to the extent to steal my money to bleed me. Read here: https://forum.wilmott.com/viewtopic.php?f=15&t=102298#p857494. British police forcefully stopped me from seeing law firms(Brief details are given in my application to HRW attached as a link below).  Since my money had given mind control agencies a very tough time, mind control agencies decided to never let me make significant money after that. My old Japanese employer was told to not hire me after I returned to Pakistan. Whoever I would contact someone in Europe for job or consulting, mind control agencies would move American embassy infrastructure in that country and somebody highly influential (or army or secret services of that particular country)  in that particular country would approach my target employer and again brazen lies and fabricated stories would be told about me shamelessly. I am a fighter in life and I continued my research despite all economic, coercive and other hardships. All odds were against me but I tried my best to continue.  I have been presenting my research to friends on this forum for past five or six years. At this point, I want to request good American and friends from other countries who have come to know me through my research to please use their influence to end my mind control forever.

Since things might heat up after my writing this and my persecution might increase, I want to ask European embassies in Pakistan to keep a vigil. I want to tell Europeans that my writing about mind control is in their great interest since the very same people who instigate retarding of intelligent muslims are behind retarding of dozens of intelligent europeans and once these people are exposed I am sure they will be far more careful in openly retarding European talent(that they do not like). 

I had told friends earlier that a jewish billionaire is Godfather of mind control. When my persecution started (in late 1997) at the university in New York which is a big time player in mind control of Muslims, blacks and foreigners, many university officials were very pleased that they were getting written in good books of the billionaire jewish Godfather. Though there are far richer billionaires in US who are staunchly opposed to mind control, they are not ready to shell out hundreds of millions of dollars to bribe anyone to end mind control while mind control Godfather  continues to give large bribes to people to continue mind control. Many of the people related to defense in mind control are given extremely high paying jobs in Godfather and his other jewish friends' companies after they retire from mind control and Godfather is extremely generous towards them and it is well established and well known to people in mind control that after retarding muslims, blacks and others in their career, they will move to paradise after they retire. No wonder there is no let up in mind control whatever good people try.

We know in almost all classic folk lore stories evil is too strong and too sure of itself while good is mostly weak and just trying to protect itself from evil who is bent on damaging and eradicating the weak and the poor good but it is the weak and poor good that finally prevails. I want to tell good American friends that evil at this stage has become too bold and thousands of talented people are retarded every year. And it is just not muslims and blacks that are on the line, several hundreds of brilliant whites are retarded across United States by mind control agency every year that has become too emboldened due to lack of any accountability. Thousands of foreigners are also kept retarded in dozens of countries including many European countries. 

My only intention behind this post is that I want mind control to end from our societies. I do not have any joy in mentioning that there are some bad jewish actors behind it. In American society there are a large number of jews who have excelled in the professions of science, education, medicine and have made great contributions to our societies. There are a very large number of jewish professors each of whom would have taught tens of thousands of students on top of making extraordinary contributions to body of human knowledge through their great research. Then there are accomplished jewish scientists whose work, research and inventions make our everyday life better. There are great jewish doctors and surgeons who would treat countless sick people and make them recover from disease. The list of contributions of jewish people to humanity is very long and it continues. Then there are a large number of ordinary jewish folks who are known for their goodness and civility to other people in the society. When we know all of this, we do not want to say anything anti-jewish to hurt the feelings of so many good jewish people. But I still have to say there are really some bad people who are unconnected with most of the nice and accomplished jews and they instigate and continue cruel mind control torture on talented people of other ethnicities, religions and countries based on their personal dislike and hatred. My purpose is only to end the mind control victimization and torture on otherwise nice and talented people and I do not have any hidden intentions beyond that. And I am sure American society knows very well about the contributions of good Jews to both their country and broader humanity and my exposing some bad actors would never have any bad effect of any kind on these good and nice jews who we all respect and who are our role models.I want to tell friends something I have said before here that I am sure American society would eventually end the menace of mind control from their and other countries. I want to request all good Americans and Europeans to come forward and invite any investigative journalist you know to read my threads and ask them to help end the menace of mind control in our societies through their investigative journalism. Without any accountability evil will continue to grow stronger and then mind control would start to threaten even those people who consider themselves safe due to affiliations, status or wealth.

I want to request CNN, New York Times, Washington Post and large reputed European media outlets to run a comprehensive story on mind control torture, retarding and victimization of intelligent people of US and other nations by crooks in US army on behest of some powerful people in United States and due to rightwing extremist biases of these crooks. If you would like to do investigative research about animal practices of US army, one great resource would be mainland European embassies in Muslim countries who keep a detailed account of mind control persecution of intelligent muslims in these countries by crooks in US army. Many of the staff in mainland European embassies are very good human beings who abhor such practices and would love to cooperate with good journalists in exposing the evil animal practices of US army crooks. Only the accounts of people at mainland European embassies in muslim countries thorough animal practices of American army's mind control wing would be enough to drop a great bombshell in the media and general public all across the world. I want to warn all the good people who try to have a civilized dialogue with crooks in American army to end their evil practices that their being nice with crooks is a very misguided approach. Since good people do not have the power to forcefully end the evil practices of US army crooks, only way to end the evil practices of US army crooks is by exposing them openly in US public and all across the world. 
It seems that crooks in American army and mind control  are unable to pay a heed to civilized calls by good American people and others to end their animal cruel practices and remain bent on mischief. If crooks in American army and mind control remain adamant on continuing evil practices, the only solution seems that some brave investigative journalists expose them by running a comprehensive story on mind control practices by American army. It seems that doing a civilized dialogue with hardened crooks in mind control gives them a strong sense of weakness of good people and makes the crooks even more adamant to continue their evil practices.  Reminds me of  Abu-Gharib. Very similarly, Crooks in defense had no conception that they were doing any wrong thing when there was a culture of openly urinating on human captives. It was an open thing in army and most crooks thought it was indeed a very right thing to do( and I am sure some of those crooks still believe that it was a right thing they did even after being reprimanded by the broader American society). It was not until brave people at CNN did a daring story against animal practices by crooks that evil practices stopped. Though American army crooks retard intelligent people of all color and creed, these ultra right wing army crooks love to retard blacks and muslims with great relish. Blacks are only 8-10% of US population but make a very large proportion of victims in united states since many ultra right wing crooks in US army would rather die than let those intelligent blacks succeed in American society in a big way.
For the context of my persecution, please see
In this thread I have given some account of my persecution spread over past 23 years.

My Letter to HRW About My Human Rights Abuse : https://forum.wilmott.com/viewtopic.php?f=15&t=102298

In the thread below I have given some discrete hints about why American defense is so rabid to persecute me and other muslims like me(and blacks) (and why no president finds it politically expedient to end mind control on innocent human beings.)

Muslims, mind control, Muslim-Jewish relationship and American defensehttps://forum.wilmott.com/viewtopic.php?f=15&t=101453

In the thread below I tried in vain to request prime minister of my country to end mind control here but he is helpless before crook generals in Pakistan army who take hundreds of millions of dollars in bribes from CIA. 
My Open Letter to Prime Minister of Pakistan, Imran Khan, to Intervene In Order to Protect the Interests of Pakistan  https://forum.wilmott.com/viewtopic.php?f=15&t=94796&start=1125#p858119

I want to end this post by re-iterating my apology to good jews but want to say again there are enough evil jews who are doing every thing to tarnish the great goodwill jews have so rightly earned through generations of hard work and being nice and kind to other people in the society. I want to request good jews to Please force these evil(jewish) people to end the bad practices that are starting to take root in American and other societies on their behest and before it is too late.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 24th, 2021, 9:28 am

This is another post I am copying from before for reference of friends.

Thank you Paul for your post. I will look at the website you mentioned. But I want to really thank you for letting me write bold things on your forum that most forum owners would not allow. Though I might have different ideas than you on a lot of issues, I am a staunch supporter of freedom of speech just like you are and I am sure both of us agree on one slogan,"Je suis Charlie". I consider myself lucky to have been writing on your forum since 2003.

I want to write this post to ask American tax payers to please cut funding of crooks in US army. I was reading that there has not been a successful audit of spending by US army in past decade. You would probably ask me why you should do anything to bridle spending by crooks in army. Here is why.
Crooks in US defense have spent several billion dollars in attempts to retard me over past twenty years. I went to Hong Kong in 2000 to find a job in finance and derivatives and I stayed there for three months. Crooks in US army continues to drug food and water at large number of places in Hong Kong and Kowloon and continued to play all sort of tricks with me. I still recall when I went to office of South China Morning Post to tell them how badly they had drugged the city of Hong Kong, they had already set up someone who received me and told me after hearing my story that he did not believe me. 
They gave me electric shocks twenty days after I returned from Hong Kong.
Later I was asked by my Japanese employer to come to Tokyo (2005) and work in their office there(I was earlier working from Lahore). They drugged large number of places in Tokyo city and mad e my life absolute hell and I resigned only after one month of stay in Tokyo. I did not want to return to Pakistan due to huge persecution by my family and the army so I decided to go somewhere where US military could not go and got Iranian visit visa from Tokyo and went to Tehran
Crooks in US defense gave huge amount of money to Iranian crooks(crooks everywhere know each other well) and they used huge microwave and gases torture on me. I still recall on my first night in an Iranian hotel, I left the hotel after midnight and randomly ran in streets of Tehran to somehow avoid pain due to directed microwave torture. Just after two or three days in Tehran, I was trying to find refuge in embassies like Vatican in Tehran asking them to help me on human grounds. I left Tehran after one week and flew back to Pakistan. I still recall that I arrived at Tehran airport a lot earlier than flight departure. I had previously protested somewhere against unnecessary removal of shoes of Muslim passengers and it was known to mind control agencies. They asked Iranians security on airport who approached me inside the departure lounge and took me to a room and asked me to remove my shoes and had my body search.
Since during all this, mind control agencies could not control me well and several of my neurotransmitters were coming back, they had extremely huge amount of gases in washrooms in Dubai airport where I would stop for transit while going to Karachi. No wonder ultra conservative crooks in US army align so well with Saudi(who dissolve their criticizers in acid) and Dubai.
Later I got British HSMP work visa and went to London in 2010. Crooks in US defense drugged several large cities of UK on a large scale(it was impossible to get good beverages or bottled water. They even drugged beer, wine and pork at a lot of places.) You can ask people in UK secret services(many of whom would love to tell you everything in detail) about supernatural stories how they drugged food in cities of UK in 2010. I still recall that I wrote several page account of my persecution and I would print hundred of copies of it and distribute them in London. When I was distributing those papers outside of some London Newspapers, administration of those newspapers went to extreme length to point out that certain part of sidewalk was their property and they would call police if I stepped there. I was forcefully removed from the compound in London(after I had passed all security checks) when I went there to ask law firms for help and the law firms on main street would refuse to meet me and say that I had to either email or send the letter by post and they would not accept any papers from me(somebody would approach me at the door and tell me about it and ask me to leave.). 
They continued to drug food and beverages in entire city of Lahore numerous times in past ten years.
Net worth of several mind control agents who continued to retard me in Pakistan is north of several hundred million dollars each depending upon their seniority(with most senior agent worth 800 million USD).
And this drugging of entire cities by crooks of US defense never ends. Now they are drugging the twin cities of Rawalpindi/Islamabad.
I am just telling about myself. There are tens of thousands of mind control victims all over the world.
Please cut their funding before it is too late. They are not retarding just Muslims(I have no religion but of course I am Muslim by birth) only, they retard a large number of totally innocent and talented Americans. Crooks love to bite the same hand that feeds them(It is typical of crooks and large unbridled armies in human history)
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 26th, 2021, 11:14 am

I want to tell friends that I am on antipsychotics since spring 1999 when I came back to Pakistan( with some breaks for example when I came to London or went to Tokyo and Hong Kong). My family was lured by grand promises made to them and continued to force antipsychotics on me and it was also very easy for them since most Pakistan army generals are bribed tens of millions of dollars each by mind control agencies and they force psychiatrists to declare me mentally ill when I talk about mind control. I do everything normally and perfectly well, I have perfect and very organized language skills. I love my life and have great regard for other people's life and their dignity. I am very cool tempered and almost never get upset or angry and even if I am upset, I have complete restraint over myself and I do not show my anger. And eleven years ago besides all of the above I was making USD 150,000 per year in Pakistan but still psychiatrists under the coercion of Pakistan military declared me severely mentally ill.
Ok, coming to gist of it, I just want to enjoy my life whatever of it is left. But when I am given antipsychotic injection, I do not even like music for many days. I love mathematics and I really enjoy working with SDEs and stochastics and I would have been a miserable person if it were not for love of doing mathematics despite that my facility of thought greatly decreases when I get antipsychotic injection. I have lived with bare minimum money for past eleven years and if it were not for antipsychotics I would still be capable of thoroughly enjoying my life with little money since I enjoyed doing mathematics. Though I think making money should not be a big problem for me in near future since I have successfully worked with fantastic market trading algorithms (which would never have been possible if I had been working on these algorithms all my life. These algorithms became possible since I exercised my brain on mathematics and some ideas spilled over to market trading though the algorithms have little to do with SDEs), even if I have very good money, I want to truly continue doing research on mathematics since I truly enjoy it. I am not in the business of maligning anyone and would love to spend most of the rest of my life doing mathematics. And I would love to continue sharing my research in mathematics with friends. And again I want to enjoy the rest of my life with my full brain without antipsychotics and mind control. If I were given the choice of huge money with antipsychotics and mind control as opposed to full brain without antipsychotics but little money, I can never imagine to choose money. I hope and plan when I would have more money due to my market trading algorithms, I would still want to spend most of my time doing research in mathematics. 
I have absolutely no hard feelings towards anyone despite all the trauma, torture and coercion I have seen in my life. I can imagine things might have been slightly different and I might naturally have had hard feelings because of gross injustice meted out to me if I had not found expression in my research for mathematics. As an aside, I want to suggest if somebody is in the low of his/her life due to failure or ill luck, please try to have a passion for some positive activity that you like and you would be able to simply shrug your ill luck or failures and enjoy the life more. To enjoy life, we have to exercise our brain, and the more connected (within its various parts) it is, the more we would be able to enjoy life. 
I used to have a religion but scientific facts and evidence tells me otherwise so I do not believe in religion any more for quite sometime now. 
Last of all I want to apologize to all good jews who might not have liked my previous posts and want to tell them that I harbor absolutely no resentment against jews and consider them friends. From my irreligious perspective, I want to tell them that Jews and Muslims in Israel and Palestine share the same land and all they have to learn is live side by side peacefully and have regard for each other and for each other's dignity. If you ever have a religion, it should better teach love for other human beings irrespective of their beliefs. 
So I want to say again that I have absolutely no ill will or resentment about any wrong done to me, or injustice or ill fate that I had. I want to move on with my life with continuing research in mathematics and truly enjoy my life and doing mathematics. So I make this request to good Americans and also to good Jews to please do whatever you can to end my mind control and ask the mind control agencies to not force antipsychotics on me anymore so I could live the rest of my life with my human dignity intact and enjoy living it.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 28th, 2021, 10:15 pm

I am here in Kot Addu for past several days. It is our old home where I have come after about two years and every place was not properly covered for mind control. Though I slept late, I had a very great sleep last night. And I was able to work on my research very well for the day after that. It is 3:00 am now in the night and I slept today a bit after ten since I was very tired after working on my research continuously during the day. But mind control agencies continued micro wave torture on me and I had an extremely poor and disturbing sleep due to directed microwave on my head. Therefore I had to write this post in the mid of the night due to extreme frustration.  It seems that there is no plan to decrease torture on me. I want to take this opportunity to warn the evil and malicious jews that if my tone over the past few years remained good and I showed restraint, it was only because of the good people among jews that I greatly regard, I want the evil and malicious jews to rest assured that I would continue to expose their hideous tactics and true  reality for the rest of my life. If the evil billionaire jewish godfather thinks, he would continue to retard muslims and others with impunity, he better think think again. I will come back with detailed posts in a few days again. I never wanted to say this but it was similar practices that evil Jews had started in Germany where rich jews got together and failed and squeezed every non-jew on Berlin stock exchange. As the irony goes all of these evil jews made their way successfully to America (with all their riches)before anything could happen to them and all the good jews who could never believe to hurt anyone or wrong anyone were left behind to burn in concentration camps(to pay for the crimes of evil jews) while all the evil ones(who were cornering others) had already found their way out to America with their riches and remain unscathed. I want to tell good jews that mind control victimization in today's America is a truth and please stop (the evil jewish crooks) it.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 29th, 2021, 2:40 am

Friends, I have written the program series solution to SDEs and their density that is accurate to dt^2. I continued to work on it yesterday and was able to complete it. I want to make a few more checks on the program today and will post it later during the day or tonight. I have changed the structure of the program. Though program is very simple, I still want to give enough explanations so that friends could relate equations in my next post to implementation of those equations in the program. 
I hope to post the new program this evening. 
Exciting thing is that we can write equations for the greeks associated with the SDE and parameters of the SDE and run their series solution in parallel to the original calculation of the solution of  SDE and be able to relate the solution of greeks to the solution of SDE in a pointwise sense through a common Z.
 
User avatar
Amin
Topic Author
Posts: 2904
Joined: July 14th, 2002, 3:00 am

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

November 29th, 2021, 11:45 am

Friends I am posting some graphs of solution of densities of the SDEs using new Z-series method. In all of these graphs step size is 1/32. I could even make most of the graphs below with step size 1/16 but then means would be off by another 3-5 basis points. The SDEs I have chosen are those that do not have extremely high volatility and do not go to zero. In some cases we may still need a step size close to 1/128. Volatility exponents close to .5 near zero are still a problem but I am sure they can be handled after some further research. For example, I was able to devise several ad-hoc methods that would make the Z-series solution of densities close to zero stable but those methods still were not perfect. For example we could restrict absolute value of (n+1)th coefficient of series to a limit of 1/5 of absolute value of nth coefficient of the Z-series solution on each time step and do not let higher powers coefficients become bigger than a fraction of consecutive lower power coefficients of the Z-series. Such methods would work many times but still would fail sometimes. Since my approach was just ad-hoc, some research has to be done for  proper growth restrictions on coefficients required for stability. Another thing I noticed was if we decreased series order from six to four and decrease time step a bit, many caases(though not all) of volatility exponents close to .5 near zero would become stable.
I will be posting the program with comments in a few hours.
On my very slow computer a typical series solution out to five years was taking a time between .1 to .2 seconds in matlab (with loops everywhere) when step size was 1/32 as in graphs below.


Image

Image

Image

Image

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

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

November 29th, 2021, 1:28 pm

I am posting the program used to create above graphs. I will follow with explanation in the next post. Once you run this program, it will produce the last graph in previous post.
.
function [] = FPESeriesCoeffsSolutionWilmottOdt2()

%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/2;
TtM=T/dtM;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46;  % 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=.1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.5;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=4.0;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.55;%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)-3.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,ExpnOrder+2);
     
     
     %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] = CalculateDriftbCoeffs04(wMu0dta,dwMu0dtdwa,a);
     [b10,b1] = CalculateDriftbCoeffs04(wMu1dt,dwMu1dtdw,a);

     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);
    c(ExpnOrder+1:10)=0.0;
%     if(gamma<=.75)
%     for mm=1:7
%         %if((sign(c(mm))~=sign(c(mm+1)))&&(abs(c(mm))<abs(c(mm+1)*5)))
%         if((abs(c(mm))<abs(c(mm+1)*5)))    
%             %c(mm+1)=c(mm+1)*abs(c(mm)/c(mm+1)/11);
%             %c(mm+1)=sign(c(mm+1))*abs(c(mm)/4.5);
%             c(mm+1)=sign(c(mm+1))*abs(c(mm+1)/5.5);
%         end
%     end
%     end
    

    
    w0=c0;
    
    end

    a0=c0;
    a(1:ExpnOrder)=c(1:ExpnOrder);
    a(ExpnOrder+1:10)=0;
    

end    
wnStart=1;

w(wnStart:Nn)=c0;
for nn=1:ExpnOrder
    w(wnStart:Nn)=w(wnStart:Nn)+c(nn)*Z(wnStart:Nn).^nn;
end
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(wnStart:Nn) =(c(1)+2*c(2)*Z(wnStart:Nn)+3*c(3)*Z(wnStart:Nn).^2+ ...
%         4*c(4)*Z(wnStart:Nn).^3+5*c(5)*Z(wnStart:Nn).^4+6*c(6)*Z(wnStart:Nn).^5+ ...
%         7*c(7)*Z(wnStart:Nn).^6+8*c(8)*Z(wnStart:Nn).^7+9*c(9)*Z(wnStart:Nn).^8+ ...
%         10*c(10)*Z(wnStart:Nn).^9)*dNn;

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
yy0


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(2*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
.
.
.
function [wMu0dt,dwMu0dtdw,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=3;

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;


YqCoeff(1)=mu1;
YqCoeff(2)=mu2;
YqCoeff(3)=-.5*sigma0^2*gamma;
Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;



wMu0dt0=0;
dwMu0dt(1:ExpnOrder+1)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder+1)=0.0;


for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp(mm);
    for nn=1:ExpnOrder+1
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp(mm)*(1-gamma)/((1-gamma)*w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp(mm)-(nn-1))*(1-gamma)/((1-gamma)*w0);
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder+1
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
wMu1dt=dwMu0dtdw(1);
dwMu1dtdw(1:ExpnOrder)=dwMu0dtdw(2:ExpnOrder+1);
    
end

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



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


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


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




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

b(9:10)=0.0; 
 
 
 % 
% 
% 
% 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 [c0,c] = EvolveSDESeriesSolution08New(a0,a,b0,b,b10,b1,sigma0,dt,ExpnOrder)


SeriesOrder=ExpnOrder;
d1(1:ExpnOrder)=0;
d10=a(1);
for nn=1:ExpnOrder-1
    d1(nn)=(nn+1)*a(nn+1);
end
d2(1:ExpnOrder)=0;
d20=2*a(2);
for nn=1:ExpnOrder-2
    d2(nn)=(nn+1)*(nn+2)*a(nn+2);
end
d3(1:ExpnOrder)=0;
d30=6*a(3);
for nn=1:ExpnOrder-3
    d3(nn)=(nn+1)*(nn+2)*(nn+3)*a(nn+3);
end

[d10Sq,d1Sq] =SeriesProduct(d10,d1,d10,d1,SeriesOrder);
[d10Cube,d1Cube] =SeriesProduct(d10Sq,d1Sq,d10,d1,SeriesOrder);
[d10Quad,d1Quad] =SeriesProduct(d10Cube,d1Cube,d10,d1,SeriesOrder);
[d10Inv,d1Inv] = SeriesReciprocal(d10,d1,SeriesOrder);
[d10SqInv,d1SqInv] = SeriesReciprocal(d10Sq,d1Sq,SeriesOrder);
[d10CubeInv,d1CubeInv] = SeriesReciprocal(d10Cube,d1Cube,SeriesOrder);
[d10QuadInv,d1QuadInv] = SeriesReciprocal(d10Quad,d1Quad,SeriesOrder);
[d20Sq,d2Sq] =SeriesProduct(d20,d2,d20,d2,SeriesOrder);

Z_d1Inv0=0;
Z_d1Inv(1)=d10Inv;
Z_d1Inv(2:SeriesOrder)=d1Inv(1:SeriesOrder-1);

[d2_d1SqInv0,d2_d1SqInv] =SeriesProduct(d20,d2,d10SqInv,d1SqInv,SeriesOrder);
T_Odt0=b0+.5*sigma0^2*d2_d1SqInv0;
T_Odt(1:SeriesOrder)=b(1:SeriesOrder)+.5*sigma0^2*(Z_d1Inv(1:SeriesOrder)+d2_d1SqInv(1:SeriesOrder));


[d2_d1CubeInv0,d2_d1CubeInv] =SeriesProduct(d20,d2,d10CubeInv,d1CubeInv,SeriesOrder);

Z_d2_d1CubeInv0=0;
Z_d2_d1CubeInv(1)=d2_d1CubeInv0;
Z_d2_d1CubeInv(2:SeriesOrder)=d2_d1CubeInv(1:SeriesOrder-1);

[d3_d1CubeInv0,d3_d1CubeInv]=SeriesProduct(d30,d3,d10CubeInv,d1CubeInv,SeriesOrder);
[d2Sq_d1QuadInv0,d2Sq_d1QuadInv]=SeriesProduct(d20Sq,d2Sq,d10QuadInv,d1QuadInv,SeriesOrder);

T_Odt2(1:SeriesOrder)=b1(1:SeriesOrder)+.5*sigma0^2*(d1SqInv(1:SeriesOrder)-Z_d2_d1CubeInv(1:SeriesOrder)+ ...
    d3_d1CubeInv(1:SeriesOrder)-2*d2Sq_d1QuadInv(1:SeriesOrder));

T_Odt20=b10+.5*sigma0^2*(d10SqInv+d3_d1CubeInv0-2*d2Sq_d1QuadInv0);



[Term20,Term2]=SeriesProduct(T_Odt0,T_Odt,T_Odt20,T_Odt2,SeriesOrder);

c0=a0+T_Odt0*dt+Term20*dt^2/2;
c(1:SeriesOrder)=a(1:SeriesOrder)+T_Odt(1:SeriesOrder)*dt+Term2(1:SeriesOrder)*dt^2/2;


end
.
.
.
function [c0,c] = SeriesReciprocal(a0,a,SeriesOrder)

c0=1/a0;
c(1:SeriesOrder)=0;
for nn=1:SeriesOrder
    c(nn)=c(nn)-c0/a0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)-c(kk)/a0*a(nn-kk);
    end
end
        
end
.
.
.
function [c0,c] =SeriesProduct(a0,a,b0,b,SeriesOrder)


c0=a0*b0;
c(1:SeriesOrder)=0.0;
for nn=1:SeriesOrder
    c(nn)=c(nn)+a0*b(nn);
    c(nn)=c(nn)+b0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk);
    end
end

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

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

November 30th, 2021, 12:49 pm

Friends, I was able to somewhat improve the performance of the program for SDEs close to zero. I placed a simple constraint on the growth of derivatives of the drift that ratio between absolute values of consecutive derivatives cannot be larger than ratio of absolute values of drift itself and its first derivative. This constraint helped make the program stable for slightly larger range of SDEs close to zero. Another good thing about this constraint was that it did not seem to (negatively) affect the performance of program for densities away from zero that were already good without this constraint. I made some graphs for SDEs close to zero and I am attaching them. I was not able to find the solution for densities of SDEs with the parameters in these graphs with yesterday's program. I will post the new program in next post. For all of these graphs maximum order of series was only four and time step size was still 1/32 and terminal time is five years.

Image



Image

Image

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

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

November 30th, 2021, 1:00 pm

Here is the new program for calculation of derivatives of drift in which I have implemented simple bounds.
function [wMu0dt,dwMu0dtdw,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08Bounded(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=3;

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;


YqCoeff(1)=mu1;
YqCoeff(2)=mu2;
YqCoeff(3)=-.5*sigma0^2*gamma;
Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;



wMu0dt0=0;
dwMu0dt(1:ExpnOrder+1)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder+1)=0.0;


for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp(mm);
    for nn=1:ExpnOrder+1
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp(mm)*(1-gamma)/((1-gamma)*w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp(mm)-(nn-1))*(1-gamma)/((1-gamma)*w0);
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder+1
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end

    
end


Bound=1/abs(wMu0dt/dwMu0dtdw(1));

for nn=1:ExpnOrder
   
    if(abs(dwMu0dtdw(nn+1))>abs(dwMu0dtdw(nn)*Bound))
        dwMu0dtdw(nn+1)=sign(dwMu0dtdw(nn+1))*abs(dwMu0dtdw(nn))*Bound;    
    end
end



wMu1dt=dwMu0dtdw(1);
dwMu1dtdw(1:ExpnOrder)=dwMu0dtdw(2:ExpnOrder+1);

%wMu0dt
%dwMu0dtdw

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

end
,
Here is the main program that now works well with four series order.
.
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;
TtM=T/dtM;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46;  % 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=.10;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=.0150;%mean reversion target
sigma0=.5;%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)-3.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=4;


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] = CalculateDriftbCoeffs04(wMu0dta,dwMu0dtdwa,a);
     [b10,b1] = CalculateDriftbCoeffs04(wMu1dt,dwMu1dtdw,a);

     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);
    c(ExpnOrder+1:10)=0.0;
%      for mm=1:ExpnOrder-1
%          %if((sign(c(mm))~=sign(c(mm+1)))&&(abs(c(mm))<abs(c(mm+1)*5)))
%          if((abs(c(mm))<abs(c(mm+1)/c0/(1-gamma)*(mm-1))))    
%              %c(mm+1)=c(mm+1)*abs(c(mm)/c(mm+1)/11);
%              c(mm+1)=sign(c(mm+1))*abs(c(mm)*c0*(1-gamma)/(mm-1));
%              %c(mm+1)=sign(c(mm+1))*abs(c(mm+1)/7);
%          end
%      end
    T1=sqrt(((1-gamma)^2+gamma^2)/(gamma*(1-gamma)))*(2);
    
     for mm=1:ExpnOrder-1
    %T1=(1/(1-gamma))^(mm-1)/sqrt(2);
    T1=1/(1-gamma)*(2);
         %if((abs(c(mm))<abs(c(mm+1)/c0*(mm))))
         if((abs(c(mm))<abs(c(mm+1)/c0*(mm)*T1)))    
             
%             c(mm+1)=sign(c(mm+1))*abs(c(mm)*c0/(mm)/T1);
            
         end
     end
    
%dwMu0dtdwa
    
    
    
    
%     if(gamma<=.75)
%     for mm=1:7
%         %if((sign(c(mm))~=sign(c(mm+1)))&&(abs(c(mm))<abs(c(mm+1)*5)))
%         if((abs(c(mm))<abs(c(mm+1)*5)))    
%             %c(mm+1)=c(mm+1)*abs(c(mm)/c(mm+1)/11);
%             %c(mm+1)=sign(c(mm+1))*abs(c(mm)/4.5);
%             c(mm+1)=sign(c(mm+1))*abs(c(mm+1)/5.5);
%         end
%     end
%     end
    

    
    w0=c0;
    
    end

    a0=c0;
    a(1:ExpnOrder)=c(1:ExpnOrder);
    a(ExpnOrder+1:10)=0;
    

end    
wnStart=1;

w(wnStart:Nn)=c0;
for nn=1:ExpnOrder
    w(wnStart:Nn)=w(wnStart:Nn)+c(nn)*Z(wnStart:Nn).^nn;
end
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(wnStart:Nn) =(c(1)+2*c(2)*Z(wnStart:Nn)+3*c(3)*Z(wnStart:Nn).^2+ ...
%         4*c(4)*Z(wnStart:Nn).^3+5*c(5)*Z(wnStart:Nn).^4+6*c(6)*Z(wnStart:Nn).^5+ ...
%         7*c(7)*Z(wnStart:Nn).^6+8*c(8)*Z(wnStart:Nn).^7+9*c(9)*Z(wnStart:Nn).^8+ ...
%         10*c(10)*Z(wnStart:Nn).^9)*dNn;

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
yy0


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(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