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

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

The program has the right algorithm for linear-like SDEs but my algorithm using  sqrt(1/T) *N  for each step of normal density evolution probably does not seem to hold for very non-linear SDEslike stoch vol SDEs.
But I am giving the simulation algorithm in the hope that it would be helpful to provide a simple framework for people to experiment on their own and try their ideas where Ito-Taylor coefficients and simulation program has been calculated and normal variance structure or hermite polynomial type might have to be altered. My first thought would be that trying Brownian motion based hermite polynomials (orthogonal hermite polynomials that take Brownian motion variance) might really fix the problem with very non-linear SDEs and I hope that this simulation framework would be helpful for people to try brownian motion based hermite polynomials and other ideas they have. And it would be easier for many people to slightly change the existing program and try the ideas they have.
I have not been able to work for past few weeks for many reasons and would love to try many of these ideas on my own. But, of course, my personal inability to work must not stop other friends from working on their research ideas. And I would request friends to give me credit for the original basic work I have already done and would absolutely welcome all good research from other friends.

Here is the basic program explanation.
In the expansion of X(t)^alpha
With SDE
dX(t)=mu1 * X(t)^beta1 dt + mu2* X(t)^beta2 dt + sigma * X(t)^gamma dz(t)
The first order expansion of X^alpha with the above SDE is
X(t+1)^alpha= X(t)^alpha +alpha* X(t)^(alpha-1) * mu1 * X(t)^beta1 *dt
+alpha*X(t)^(alpha-1) * mu2 * X(t)^beta2*dt
+ .5* alpha *(alpha-1)* X(t)^(alpha-2) * sigma^2 * X(t)^(2*gamma) *dt
+ alpha * X(t)^(alpha-1) * sigma * X(t)^gamma *dz(t)
In my program, for the above four term expansion, I use a four dimensional array notation that shows the order of differentiation on each of these four types of terms. So if l1 represents drift one, l2 represents drift two and l3 represents quadratic variation term and l4 represents volatility term, we can use the notation Y(l1,l2,l3,l4) for the coefficients of expansion of the above term. So in the above first order expansion
Y(1,0,0,0)=alpha * mu1
Y(0,1,0,0)=alpha *mu2
Y(0,0,1,0)=.5* alpha *(alpha-1)* sigma^2
Y(0,0,0,1)= alpha * sigma

For higher orders of expansion, we will have to carefully compute the coefficients separately and then lump them together in Y according to array indices. Once we have done that, we can write one term of the general expansion as
Y(l1,l2,l3,l4) * x^(alpha+l1 * beta1 +l2* beta2 +l3* 2* gamma +l4 * gamma
- l1 -l2 - 2*l3 -l4) * HermiteP(l4) * dt^(l1 + l2 +l3 + .5*l4);

In my program, since matlab does not take indices starting from zero, I have slightly adjusted the above notation.  In the density simulation program, we pre-compute everything but we have to calculate the part of above expression where power of x is evaluated.
If we have four terms in Ito-hermite SDE expansion at first order, In every next level, each of the four terms would take four descendent terms and there will be sixteen terms on second level, and on third level there will be 64 descendent terms and 256 descendent terms on fourth level. However we can reduce the number of terms in Ito-hermite expansion to number of terms in polynomial expansion. For example Ito-hermite expansion of SDE with four terms on first order can be reduced to the number of terms in expansion of a polynomial with four terms on each level. To do this, we will have to pre-calculate the Coefficients carefully and the lump them together as in a multinomial like expansion.
Again the above array Y(l1,l2,l3,l4) is relatively sparse and does not take all possible indices. Just like in polynomial expansion, on a certain expansion order, it takes only indices so that sum of indices goes to expansion order. For example on second order, only those indices could be found such that l1+l2+l3+l4=2. And similarly for third ito-hermite expansion level only those indices would be populated such that l1+l2+l3+l4=3.
In order to do that I have done simple loops in a fashion shown below

for k = 0 : (Order4)
for m = 0:k
l4 = k - m;
for n = 0 : m
l3 = m - n;
for j = 0:n
l2 = n - j;
l1 = j;
First loop variable is for expansion order. And then rest of the loop variables take values such that their sum equals the expansion order. When we would add a fifth term in the SDE expansion for example, we will have to add another loop variable with the same logic that sum of loop variables equal to expansion order. In my matlab program, I have given offset to l4, l2,l3, l1 since array indices cannot start at zero.

In the coefficient calculation program which calculates Y(l1,l2,l3,l4), I have used four levels of looping each for relevant expansion order. The first loop takes four values and second loop takes 16 values and third loop takes 64 values and so on. And then each coefficient term can be individually calculated while carefully accounting for path dependence.

So for example in a nested loop structure

m1= 1:mDim
m2=1:mDim
m3=1:mDim

l(m1)=l(m1)+1;
l(m2)=l(m2)+1;
l(m3)=l(m3)+1;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

in the above looping loop takes values from one to four with one indicating the first drift term, two indicating the second drift term and three indicating quadratic variation term and four indicating the volatility term. And with this looping structure we can
So in the above looping m1=1 would mean that all terms are descendent of first drift term and m2=4 would mean that all terms are descendent of first drift term on first expansion order and descendent of volatility term on the second order and so we can keep track of path dependence perfectly.

And on each level, we individually calculate the descendent terms. While keeping track of path dependence and calculating the coefficients with careful path dependence consideration, we update the appropriate element in our polynomial like expansion coefficientarray .Here l(1) denotes l1 but written as l(1) so it can be conveniently updated with the loop variable when the loop variable takes value one indicating first drift term . And l(2) could be conveniently updated when the loop variable takes value equal to two indicating second drift term and so on.

Here is the part of code snippet for that

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

The first four lines update the array indices according to the parent term. And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.
ArrIndex0=m1;calculates the array index of the parent term
And
ArrIndex=(m1-1)*mDim;calculates the array index of the descendent terms

And coefficient of the drift and volatility descendent terms is calculated by multiplying the coefficient of the parent term by
Coeff1st=Y1(ArrIndex0)*CoeffDX1;

And coefficient of the quadratic variation descendent terms is calculated by multiplying the coefficient of the parent term by

Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;

And then each of the four descendent terms are updated with Coeff1st if they are drift or volatility descendent terms or Coeff2nd if they are quadratic variation descendent terms.
Here Y1 indicates the temporary coefficient array with parent terms on first level. And Y2 denotes temporary coefficient array with parent terms on second level and Y3 indicates temporary coefficient array with parent terms on third level.

Here is the main simulation Code after the coefficients have been calculated
function [] = ItoTaylorDensitySimulation()

% In the expansion of X^alpha
% With SDE
% dX(t)=mu1 * X(t)^beta1 dt + mu2* X(t)^beta2 dt + sigma * X(t)^gamma dz(t)
% The first order expansion of X^alpha with the above SDE is
% X(t+1)^alpha= X(t)^alpha +alpha* X(t)^(alpha-1) * mu1 * X(t)^beta1 *dt
% +alpha* X(t)^(alpha-1) * mu2 * X(t)^beta2 *dt
% + .5* alpha *(alpha-1)* X(t)^(alpha-2) * sigma^2 * X(t)^(2*gamma) *dt
% + alpha * X(t)^(alpha-1) * sigma * X(t)^gamma *dz(t)

% In my program, for the above four term expansion, I use a four dimensional
%array notation that shows the order of differentiation on each of these
%four types of terms. So if l1 represents drift one, l2 represents drift
%two and l3 represents quadratic variation term and l4 represents
%volatility term, we can use the notation Y(l1,l2,l3,l4) for the
%coefficients of expansion of the above term. So in the above first order expansion

% Y(1,0,0,0)= alpha * mu1
% Y(0,1,0,0)=alpha *mu2
% Y(0,0,1,0)= .5* alpha *(alpha-1)* sigma^2
% Y(0,0,0,1)= alpha * sigma
%
% For higher orders of expansion, we will have to carefully compute the
%coefficients separately based on path dependence and then lump them
%together in Y according to array indices. Once we have done that,
%we can write one term of the general expansion as

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

Order4=4;  %expansion order
Nn=161;    % No of normal density subdivisions
dNn=.05;   % Normal density width. would change with number of subdivisions
dt=.250;   % Simulation time interval.
Tt=2;      % simulation level. Terminal time= Tt*dt;
x0=.250;   % starting value of SDE

alpha=1;   % x^alpha is being expanded
beta1=0;   % the first drift term power.
beta2=1;   % Second drift term power.
gamma=1;   % volatility power.
mu1=.25;   % first drift coefficient.
mu2=-1;    % Second drift coefficient.
sigma0=.750;%Volatility value

for nn=1:Nn
x(nn)=x0;
end

for nn=1:Nn
Z(nn)=((nn-1)*dNn-4.0)*1/sqrt(Tt);

end

ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1/sqrt(Tt));
ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1/sqrt(Tt));
for nn=2:Nn-1
ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1/sqrt(Tt))-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1/sqrt(Tt));
end

%ZProb contains probaility mass in each cell.

% Pre-compute the hermite polynomilas associated with each subdivision of
% normal.
for nn=1:Nn
HermiteP(1,nn)=1;
HermiteP(2,nn)=Z(nn);
HermiteP(3,nn)=Z(nn)^2-1;
HermiteP(4,nn)=Z(nn)^3-3*Z(nn);
HermiteP(5,nn)=Z(nn)^4-6*Z(nn)^2+3;
end

%These sigma11, mu11, mu22, sigma22 can be lumped with Y(l1,l2,l3,l4) but I
%kept them separate for a general program where these values could change
%on every step while exponents on x in each term would remain same in
%pre-calcualted Y(l1,l2,l3,l4)

sigma11(1:Order4+1)=0;
mu11(1:Order4+1)=0;
mu22(1:Order4+1)=0;
sigma22(1:Order4+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:(Order4+1)
if sigma0~=0
sigma11(k)=sigma0^(k-1);
end
if mu1 ~= 0
mu11(k)=mu1^(k-1);
end
if mu2 ~= 0
mu22(k)=mu2^(k-1);
end
if sigma0~=0
sigma22(k)=sigma0^(2*(k-1));
end
end

Ft(1:Tt+1,1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;
Fp(1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;

%Pre-compute the time and power exponent values in small multi-dimensional arrays

for k = 0 : (Order4)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
Ft(l1,l2,l3,l4) = dt^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
end
end
end
end

%Below call the program that calculates the Ito-Taylor coeffs.

Y = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
Y(1,1,1,1)=0;

for tt=1:Tt
for nn=1:Nn
x1(nn)=x(nn);
for k = 0 : Order4
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;
x1(nn)=x1(nn) + Y(l1,l2,l3,l4) * ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
x(nn)^Fp(l1,l2,l3,l4) * HermiteP(l4,nn) * Ft(l1,l2,l3,l4);
end
end
end
end
x(nn)=x1(nn);
x(nn)
end
end

for nn=1:Nn
B(nn)=((nn-1)*dNn-4.0);
end

for nn=2:Nn-1
Dfx(nn) = (x(nn + 1) - x(nn - 1))/(B(nn + 1) - B(nn - 1));
end
fx(1:Nn)=0;
for nn = 2:Nn-1
fx(nn) = normpdf(B(nn),0, 1)/abs(Dfx(nn));
end

%x

%Ft(1,1,1,2)
%Calculate the density with monte carlo below.

dt0=dt/8;
paths=100000;
X(1:paths)=x0;
Random1(1:paths)=0;
for tt=1:Tt*8
Random1=randn(size(Random1));
X(1:paths)=X(1:paths)+ mu1* X(1:paths).^beta1 *dt0 + ...
mu2 * X(1:paths).^beta2*dt0 + ...
sigma0 * X(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu1^2 * beta1* X(1:paths).^(2*beta1-1) *  dt0^2/2.0 + ...
mu2^2 * beta2 * X(1:paths).^(2*beta2-1) * dt0^2/2 + ...
mu1*mu2* beta1* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
mu1*mu2* beta2* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
sigma0^2 * gamma * X(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu1 * sigma0 *gamma* X(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
mu2 * sigma0 *gamma* X(1:paths).^(beta2+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
sigma0 * mu1 *beta1* X(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ...
sigma0 * mu2.*beta2* X(1:paths).^(beta2+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths);
end

sum(X(:))/paths   %monte carlo average
sum(x(:).*ZProb(:))  % Ito-Taylor density average(good enough but slightly approximate)

BinSize=.00125;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff);
plot(x(2:Nn-1),fx(2:Nn-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');

end


Here is the Coefficient calculation program
function [Y] = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma)

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

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

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

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

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

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

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

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

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

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

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

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

for m3=1:mDim
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
l(m1)=l(m1)+1;
l(m2)=l(m2)+1;
l(m3)=l(m3)+1;
CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
- (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
CoeffDX2=CoeffDX1-1;
ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
%ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
Coeff1st=Y3(ArrIndex0)*CoeffDX1;
Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%Y4(1+ArrIndex)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(2+ArrIndex)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(3+ArrIndex)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
%Y4(4+ArrIndex)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m3),n1(m2),n1(m1));
end
end
end
end
Here is the stochastic integrals calculation program
function [IntegralCoeff0,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs()

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

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

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

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

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

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

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

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

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

l0(1:2)=1;

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

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

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

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

for m1=1:2
l0(1)=1;
l0(2)=1;
l0(m4)=l0(m4)+1;
l0(m3)=l0(m3)+1;
l0(m2)=l0(m2)+1;
if(m1==1)
IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
end
if(m1==2)
IntegralCoeff(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
l0(m1)=l0(m1)+1;
%IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff4(m4,m3,m2,m1);
IntegralCoeffdt(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
IntegralCoeffdz(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
end
end
end
end
end
function [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti(X,Paths,BinSize,MaxCutOff )
%Processes monte carlo paths to return a series Xdensity as a function of IndexOut. IndexMax is the maximum value of index.
%
Xmin=0;
Xmax=0;
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
XDensity(1:IndexMax)=0.0;

for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);
if(real(index)<1)
index=1;
end
if(real(index)>IndexMax)
index=IndexMax;
end

XDensity(index)=XDensity(index)+1.0/Paths/BinSize;

end

IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;

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

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

It will be a few days before I post a new full fledged program since some of the program architecture has to be completely changed.

But some interesting things about hermite polynomials.

1. If we have two different stochastic integrals or dependent processes of a single SDE and we want to get a density representation of the sum/difference and we know their hermite representations, We can simply add/subtract the respective hermite coefficients as hermite is an orthogonal basis. By coefficients of hermite polynomials, I mean different variable weights that we have assigned to each hermite polynomial when we found the hermite representation. And these weights are not to be confused with the fixed coefficients assigned to different powers in each hermite polynomial.

2. The same applies if we have two different stochastic integrals or dependent processes of two different SDEs driven by a single normal random variable and we want to get a density representation of the sum/difference of stochastic processes of different SDEs driven by the same normal noise and we know their hermite representations, again we can simply add/subtract the respective hermite coefficients/weights as hermite is an orthogonal basis.

3. If we have two different stochastic integrals of totally different SDEs driven by independent normal noises, and we want to add them and get a single density representation for the sum, we could easily do that if we have respective hermite representations of each stochastic integral since there is a map from two different hermite polynomials to a single hermite polynomial of the sum of two processes.

$H_{n}(x+y) = 2^{\frac{-n}{2}} \sum_{k=0}^{n} \binom{n}{k} H_{n-k}(x \sqrt{2}) H_{k}(y \sqrt{2})$

I copied the formula from Wikipedia and I think it is based on a combination of Taylor expansion and Fourier transforms. This is a formula (and its fourier counterparts) that relates hermite representation of a sum of two independent processes in terms of their hermite representation and coefficients of the resulting process can be found from the same formula in terms of coefficients/weights of hermite polynomials of two independent processes. I am sure there can be vast applications of this formula and one very simple application can be finding a single hermite representation of sum of two processes driven by two independent lognormal SDEs or sum of various stochastic integrals of two different lognormal SDEs. And once we have a single hermite representation of a sum of two stochastic processes, we can easily calculate density of such a stochastic process. I am sure there will also be many other interesting applications of this formula.

One simple way to deal with non-linear equations would be to use some variant of the formula
$H_{n}(Z_{m+1} * \sqrt{T}) = Const * \sum_{k=0}^n \binom{n}{k} H_{n-k} (Z_{1} *\sqrt{dt} ) H_{k} (Z_{m} \sqrt{m*dt})$

where $\sqrt{T} = \sqrt{dt + m *dt}$
and all of $Z_m$ for m=0 to M are standard normal variables.

And as I mentioned in the copied post, you could calculate the coefficients of $H_{n}(Z_{m+1} *\sqrt{T})$ from the above formula with the same operations on coefficients of hermite polynomials on RHS as done for the calculation of hermite polynomial on LHS of the above equation from two polynomials on RHS and we just need to calculate the new coefficients and not try to calculate the actual hermite polynomilas as calculated in the above formula. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

Within next three to four days, I hope to post a new program that will also be able to price non-linear stochastic volatility models perfectly. It will also calculate all intermediate densities in a single loop while calculating the final density at terminal point. I will use stochastic calculus of standard deviations technology for that purpose. So it will be a combination of Ito Taylor method and stochastic calculus of standard deviations method. Cuchulainn
Posts: 59205
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

'perfectly' is not a word in mathematician' vocabulary. It is context-dependent. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

I am posting a simple program that uses the stochastic calculus of standard deviations method along with Ito-Taylor expansions. In the preliminary program, I have just used Ito-Taylor expansion till order two only. This program works with linear SDEs since I am using only second order expansion. At each stage of the time loop(or time evolution) in density simulation, you can seamlessly calculate the density of the SDEs. In two to three days, I will post another program that will have Ito-Taylor expansions up to order four and it would be able to solve non-linear SDEs easily. I will also explain the program methodology after the next post.
Here is the code. This is a first preliminary version and I will post another version with Ito-Taylor expansions in two to three days.
function [] = SDEXticsInfinitiTechnologiesBestMR( )
%This program solves SDEs dX(t)=mu1 X(t)^beta1 dt + mu2 * X(t)^beta2 dt +
%epsilon* X(t)^gamma dz(t)
%but first try it with linear SDEs with mu2=0.
%This program uses Ito-Taylor expansions to second order only.
%I will soon be posting a very general program that will use Ito-Taylor
%expansions to fourth order and will tackle non-linear SDEs
X0_max=201;%Total number of SD Fractions.
X0_mid=101;%Mid of the SD Grid. Lies at z0 by definition.
X0_mid0=101;
X0_start=1;
dX0=.04;%*8.0;%Width of SD cells in terms of Standard Deviations.
T_index=40*1;
dt=.025/1.0;
epsilon=0.50
gamma=.650;
%kappa*theta=.2
%kappa*1=1
z0=.250;
mu1=0.0;
mu2=.4;
beta1=0;
beta2=1.0;
%beta=.5;
dz_integral2(1:X0_max)=z0;
iivector(1:X0_max)=1:X0_max;
Zvector(1:X0_max)=((1:X0_max)-1)*dX0-4.0;
paths=500000; %No. of Paths for creation of Monte Carlo density by simulating z, F(z);
%zz4(1:paths)=z0;
prob2=exp(-.5*((iivector-X0_mid).*dX0).^2)./sqrt(22/7.0*2.0)*dX0;
prob2=normcdf(((iivector-X0_mid).*dX0)+.5*dX0)-normcdf(((iivector-X0_mid).*dX0)-.5*dX0);
prob0=exp(-.5*((iivector-X0_mid).*dX0).^2)./sqrt(22/7.0*2.0);
for nn=1:T_index %Loop over time.
t=(nn)*dt;
t0=(nn-1)*dt;
prob=exp(-.5*((iivector-X0_mid).*dX0).^2)./sqrt(22/7.0*2.0)/(sqrt(t)*epsilon); %Probability in each SD Grid cell.
dt=t-t0;
dz_integral2=dz_integral2+...
(mu1.*dz_integral2.^beta1+mu2.*dz_integral2.^(beta2)) * dt+ ...
epsilon.*dz_integral2.^gamma.*(sqrt(t)-sqrt(t0)).*Zvector + ...
epsilon^2.*gamma.*dz_integral2.^(2*gamma-1).* (((t-t0)/2.*((Zvector.^2)-1))- ...
Zvector.*sqrt(t0).*(Zvector.*sqrt(t)-Zvector.*sqrt(t0)))+  ...
epsilon*(mu1.*dz_integral2.^beta1+mu2.*dz_integral2.^beta2).* gamma.*dz_integral2.^(gamma-1) .*Zvector* ...
(1/sqrt(3).*(t^1.5-t0^1.5)- t0*sqrt(t) + t0*sqrt(t0))+ ...
epsilon*(mu1.*beta1.*dz_integral2.^(beta1-1)+mu2.*beta2.*dz_integral2.^(beta2-1)) .* ...
dz_integral2.^(gamma) .*Zvector* ...
((1-1/sqrt(3)).*(t^1.5-t0^1.5)- t*sqrt(t0) + t0*sqrt(t0))+ ...
mu1*mu1*beta1*dz_integral2.^(2*beta1-1) *(.5*t^2-.5*t0^2-t0*t+t0^2) +...
mu2*mu2*beta2*dz_integral2.^(2*beta2-1) *(.5*t^2-.5*t0^2-t0*t+t0^2) +...
mu1*mu2*beta1*dz_integral2.^(beta1+beta2-1) *(.5*t^2-.5*t0^2-t0*t+t0^2) +...
mu1*mu2*beta2*dz_integral2.^(beta1+beta2-1) *(.5*t^2-.5*t0^2-t0*t+t0^2) +...
.5*mu1*epsilon^2.*beta1*(beta1-1).*dz_integral2.^(beta1+2*gamma-2)*(.5*t^2-.5*t0^2-t0*t+t0^2) +...
.5*mu2*epsilon^2.*beta2*(beta2-1).*dz_integral2.^(beta2+2*gamma-2)*(.5*t^2-.5*t0^2-t0*t+t0^2) +...
.5*epsilon^3*gamma*(gamma-1).*dz_integral2.^(3*gamma-2).*Zvector * ...
(1/sqrt(3).*(t^1.5-t0^1.5)- t0*sqrt(t) + t0*sqrt(t0));

end
dt0=dt;
paths=200000;
X(1:paths)=z0;
Random1(1:paths)=0;
for tt=1:T_index
Random1=randn(size(Random1));
X(1:paths)=X(1:paths)+ mu1* X(1:paths).^beta1 *dt0 + ...
mu2 * X(1:paths).^beta2*dt0 + ...
epsilon * X(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu1^2 * beta1* X(1:paths).^(2*beta1-1) *  dt0^2/2.0 + ...
mu2^2 * beta2 * X(1:paths).^(2*beta2-1) * dt0^2/2 + ...
mu1*mu2* beta1* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
mu1*mu2* beta2* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
epsilon^2 * gamma * X(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu1 * epsilon *gamma* X(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
mu2 * epsilon *gamma* X(1:paths).^(beta2+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
epsilon * mu1 *beta1* X(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ...
epsilon * mu2.*beta2* X(1:paths).^(beta2+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ...
epsilon^3*.5*gamma*(gamma-1)* X(1:paths).^(3*gamma-2) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) +...
.5* mu1*epsilon^2.*beta1*(beta1-1)* X(1:paths).^(beta1+2*gamma-2) *dt0^2/2.0 + ...
.5* mu2*epsilon^2.*beta2*(beta2-1)* X(1:paths).^(beta2+2*gamma-2) *dt0^2/2.0 ;
end

effective_sigma=epsilon;
Jacobian_dz1(1:X0_max)=0;
Jacobian_dz2(1:X0_max)=0;
dz_integral2(dz_integral2<0)=0.0;

for ii=2:X0_max-1

Jacobian_dz2(ii)=((dz_integral2(ii+1))-(dz_integral2(ii-1)))/(2*dX0*sqrt((T_index)*dt)*effective_sigma);
if(Jacobian_dz2(ii)~=0.0)
Jacobian_dz2(ii)=1.0/abs(Jacobian_dz2(ii));
end
end
dz_integral2_avg=sum(dz_integral2.*prob2)
dz_integral2_avg0=sum(dz_integral2.*prob0.*dX0)
X_avg=sum(X)./paths
BinSize=.0025;%Here you can change the resolution of monte carlo density. You might have to look at the graph and if it is jagged and noisy, increase the bin size
%If it is made of straight line increments, decrease the bin size.

BinSize=.0125;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff);
plot((dz_integral2(X0_start+2:X0_max-2)),prob(X0_start+2:X0_max-2).*Jacobian_dz2(X0_start+2:X0_max-2),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');

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

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

The evolution equation of nth subdivision of SDE Variable $X_n(t)$ associated with nth Gaussian standard deviation is given below

$X_n(t)=X_n\left(t_0\right)+\text{\epsilon X}_n\left(t_0\right)\int _{t_0}^tdz_n(s)$          Integral(1)
$+\epsilon ^2 \gamma X_n\left(t_0\right){}^{2\gamma -1}\int_{t_0}^t\int _{t_0}^sdz_n(v)dz_n(s)$           Integral(2)
$+.5\epsilon ^3 \gamma (\gamma -1) X_n\left(t_0\right){}^{3\gamma -2}\int _{t_0}^t\int _{t_0}^sdvdz_n(s)$              Integral(3)
$+\epsilon ^3 \gamma (2\gamma -1) X_n\left(t_0\right){}^{3\gamma -2}\int _{t_0}^t\int _{t_0}^s\int _{t_0}^vdz_n(u)dz_n(v)dz_n(s)$          Integral(4)
$+.5\epsilon ^4 \gamma (\gamma -1) (3\gamma -2)X_n\left(t_0\right){}^{4\gamma -3}\int _{t_0}^t\int _{t_0}^s\int _{t_0}^v dz_n(u)dvdz_n(s)$        Integral(5)
$+\epsilon ^4 \gamma ^2(\gamma -1) X_n\left(t_0\right){}^{4\gamma -3}\int _{t_0}^t\int _{t_0}^s\int _{t_0}^v du dz_n(v)dz_n(s)$       Integral(6)

To calculate these iterated integrals, we rewrite the integrals as some of integrals starting from lower integral limits as zero so we could easily calculate them using ito isometry. here H1(n), H2(n) and H3(n) are hermite polynomials of 1st, 2nd and 3rd order associated with nth subdivision of standard normal density. In the accompanying program 1st hermite polynomial is const/prob2; 2nd hermite polynomial is (const2/prob2-1); third hermite polynomial is (const3/prob2-3*const/prob2);These are calculated using standard normal density.

Solution to integral (1)
$\int _{t_0}^tdz_n(s)=z_n(t)-z_n\left(t_0\right)$

Solution to integral (2)
$\int _{t_0}^t\int _{t_0}^sdz_n(v)dz_n(s)=\int _{t_0}^t\int _0^sdz_n(v)dz_n(s)-\int _{t_0}^t\int _0^{t_0}dz_n(v)dz_n(s)$
$=\int _0^t\int _0^sdz_n(v)dz_n(s)-\int _0^{t_0}\int _0^sdz_n(v)dz_n(s)-\int _0^t\int _0^{t_0}dz_n(v)dz_n(s)+$
$\int _0^{t_0}\int _0^{t_0}dz_n(v)dz_n(s)$
$=.5t \text{H2}(n)-.5t_0\text{H2}(n)-z_n(t_0) z_n(t)+z_n(t_0)z_n(t_0)$

Solution to integral (3)
$\int _{t_0}^t\int _{t_0}^sdvdz_n(s)=\int _{t_0}^t\int _0^sdvdz_n(s)-\int _{t_0}^t\int _0^{t_0}dv dz_n(s)$
$=\int _0^t\int _0^sdvdz_n(s)-\int _0^{t_0}\int _0^sdvdz_n(s)-\int _0^t\int _0^{t_0}dvdz_n(s)+\int _0^{t_0}\int _0^{t_0}dvdz_n(s)$
$=\frac{t^{1.5}}{\sqrt{3}} \text{H1}(n)-\frac{t_0{}^{1.5}}{\sqrt{3}} \text{H1}(n)-t_0 z_n(t)+t_0z_n\left(t_0\right)$

Solution to integral (4)
$\int _{t_0}^t\int _{t_0}^s\int _{t_0}^vdz_n(u)dz_n(v)dz_n(s)$
$=\int _{t_0}^t\int _{t_0}^s\int _0^vdz_n(u)dz_n(v)dz_n(s)-\int _{t_0}^t\int _{t_0}^s\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)$
$=\int _{t_0}^t\int _0^s\int _0^vdz_n(u)dz_n(v)dz_n(s)-\int _{t_0}^t\int _0^{t_0}\int _0^vdz_n(u)dz_n(v)dz_n(s)$
$-\int _{t_0}^t\int _0^s\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)+\int _{t_0}^t\int _0^{t_0}\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)$
$=\int _0^t\int _0^s\int _0^vdz_n(u)dz_n(v)dz_n(s)-\int _0^{t_0}\int _0^s\int _0^vdz_n(u)dz_n(v)dz_n(s)$
$-\int _0^t\int _0^{t_0}\int _0^vdz_n(u)dz_n(v)dz_n(s)+\int _0^{t_0}\int _0^{t_0}\int _0^vdz_n(u)dz_n(v)dz_n(s)$
$-\int _0^t\int _0^s\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)+\int _0^{t_0}\int _0^s\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)$
$+\int _0^t\int _0^{t_0}\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)-\int _0^{t_0}\int _0^{t_0}\int _0^{t_0}dz_n(u)dz_n(v)dz_n(s)$
$={H3}(n)\frac{t^{1.5}}{6}-\text{H3}(n)\frac{t_0{}^{1.5}}{6}$
$-{H2}(n)\frac{t_0}{2}+z_n\left(t_0\right)\text{H2}(n)\frac{t_0}{2}$
$-\text{H2}(n)\frac{t}{2}z_n\left(t_0\right)+\text{H2}(n)\frac{t_0}{2}z_n\left(t_0\right)$
$+z_n(t)z_n\left(t_0\right)z\left(t_0\right)-z_n\left(t_0\right)z_n\left(t_0\right)z_n\left(t_0\right)$

Solution to integral (5)
$\int _{t_0}^t\int _{t_0}^s\int _{t_0}^vdz_n(u)dvdz_n(s)$
$=\int _0^t\int _0^s\int _0^vdz_n(u)dvdz_n(s)-\int _0^{t_0}\int _0^s\int _0^vdz_n(u)dvdz_n(s)$
$-\int _0^t\int _0^{t_0}\int _0^vdz_n(u)dvdz_n(s)+\int _0^{t_0}\int _0^{t_0}\int _0^vdz_n(u)dvdz_n(s)$
$-\int _0^t\int _0^s\int _0^{t_0}dz_n(u)dvdz_n(s)+\int _0^{t_0}\int _0^s\int _0^{t_0}dz_n(u)dvdz_n(s)$
$+\int _0^t\int _0^{t_0}\int _0^{t_0}dz_n(u)dvdz_n(s)-\int _0^{t_0}\int _0^{t_0}\int _0^{t_0}dz_n(u)dvdz_n(s)$
$=\text{H2}(n)\frac{t^2}{24}-\text{H2}(n)\frac{t_0{}^2}{24}$
$-z(t)\text{H1}(n)\frac{t_0{}^{1.5}}{\sqrt{3}}+z_n\left(t_0\right)\text{H1}(n)\frac{t_0{}^{1.5}}{\sqrt{3}}$
$-z_n\left(t_0\right)\text{H1}(n)\frac{t^{1.5}}{\sqrt{3}}+z_n\left(t_0\right)\text{H1}(n)\frac{t_0{}^{1.5}}{\sqrt{3}}$
$+z_n\left(t_0\right)t_0z_n(t)-z_n\left(t_0\right)t_0z_n\left(t_0\right)$

Solution to integral (6)
$\int _{t_0}^t\int _{t_0}^s\int _{t_0}^vdudz(v)dz(s)$
$=\int _0^t\int _0^s\int _0^vdudz_n(v)dz_n(s)-\int _0^{t_0}\int _0^s\int _0^vdud_nz(v)dz_n(s)$
$-\int _0^t\int _0^{t_0}\int _0^vdudz_n(v)dz_n(s)+\int _0^{t_0}\int _0^{t_0}\int _0^vdudz_n(v)dz_n(s)$
$-\int _0^t\int _0^s\int _0^{t_0}dudz_n(v)dz_n(s)+\int _0^{t_0}\int _0^s\int _0^{t_0}dudz_n(v)dz_n(s)$
$+\int _0^t\int _0^{t_0}\int _0^{t_0}dudz_n(v)dz_n(s)-\int _0^{t_0}\int _0^{t_0}\int _0^{t_0}dudz_n(v)dz_n(s)$
$=\text{H2}(n)\frac{t^2}{24}-\text{H2}(n)\frac{t_0{}^2}{24}$
$-z(t)\text{H1}(n)\frac{t_0{}^{1.5}}{\sqrt{3}}+z_n\left(t_0\right)\text{H1}(n)\frac{t_0{}^{1.5}}{\sqrt{3}}$
$-t_0\text{H2}(n)\frac{t}{2}+t_0\text{H2}(n)\frac{t_0}{2}$
$+t_0z\left(t_0\right)z_n(t)-t_0z_n\left(t_0\right)z_n\left(t_0\right)$

Please ask any question and let me know how I could explain better.
In the above program, I have tried to explain the method. There were some minor errors though since I was thinking erroneously that all the stochastic integrals commute which was not true and you would have to be careful at some places but the logic remains the same. These integrals are calculated slightly differently than the way stochastic integrals are generally calculated.
I would want to stress again that all these integral in the copied post are calculated from t0 to t where t0 is the start of the previous interval and t0 is the start of the next interval and again, unlike monte carlo, these integrals do not restart from zero.
So unlike monte carlo simulations, stochastic calculus of standard deviations calculates the Ito-Taylor expansions differently in that these stochastic integrals do not restart from zero and have to be calculated from the starting time of each time interval as the lower integration limit and end time of the time interval as the upper integration limit. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

In my earlier program for stochastic volatility SDE density simulation, I have chosen a weight of $\sqrt{\frac{1}{T}}$ with standard normal N but we have to choose this weight adaptively for each subdivision in a way so that it does not affect the integrated variance. In the above line T is the number of time steps. This is why for linear SDEs with constant term structure of volatility this weight of $\sqrt{\frac{1}{T}}$ works great. But in non-linear SDEs or SDEs with different variance at each time step, we have to be careful with the choice of weights that multiply the standard normal and this weight has to be adaptively calculated so that total integrated variance with the choice of appropriate weights along the SDE evolution remains equal to true integrated variance. I hope to be writing a program like that in a few days. Whenever our choice of weights that multiply standard normal increment alters the integrated variance along the SDE evolution, our density will be off from the true density.
Sorry, may be I should not call it integrated variance since variance term and N are evolving together. On each time step, I would try normal increment scaling proportional to variance and this is why the program works for equal variance SDEs. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

Tomorrow, I will be posting a complete program with term structure of variance and how to tackle the term structure of variances with hermite polynomials and how to choose the weights and variances of hermite polynomials. I will post it tomorrow evening with detailed notes and explanations. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

I am posting the simple toy program with term structure of variance in a lognormal distribution SDE. I have used five simulation intervals in a year each with delta_t=.2;  I will explain the program tomorrow in detail.
Volatility term structure for this toy program is given as
sigma0(1)=.4;
sigma0(2)=.75;
sigma0(3)=.2;
sigma0(4)=.05;
sigma0(5)=.135;

I basically just changed the earlier program and turned the second drift to zero for this toy example.
function [] = ItoTaylorDensitySimulation010()

% In the expansion of X^alpha
% With SDE
% dX(t)=mu1 * X(t)^beta1 dt + sigma * X(t)^gamma dz(t)
% The first order expansion of X^alpha with the above SDE is
% I have shown the solution with volatility term structure with simple
% lognormal distribution in this code example.
Order4=4;  %expansion order
Nn=161;    % No of normal density subdivisions
dNn=.05;   % Normal density width. would change with number of subdivisions
dt=.20;   % Simulation time interval.
Tt=5;      % No of simulation intervals. Terminal time= Tt*dt;
x0=.250;   % starting value of SDE

alpha=1;   % x^alpha is being expanded
beta1=1;   % the first drift term power.
beta2=0;   % Second drift term power. zero for this examole.
gamma=1.0;   % volatility power.
mu1=.25;   % first drift coefficient.
mu2=0;    % Second drift coefficient. is zero in this simple program.

%Below is the monte carlo volatility term structure for five periods in a
%year. When you change number of density simulation intervals, you will
%have to change this part.
sigma0(1)=.4;
sigma0(2)=.75;
sigma0(3)=.2;
sigma0(4)=.05;
sigma0(5)=.135;
sigmaT=sigma0(1)^2*dt+sigma0(2)^2*dt+sigma0(3).^2*dt+sigma0(4).^2*dt+sigma0(5).^2*dt;
%sigmaT
%1/sigmaT
T1(1)=sigma0(1)^2*dt/sigmaT/.2/.2;
T1(2)=sigma0(2)^2*dt/sigmaT/.2/.2;
T1(3)=sigma0(3)^2*dt/sigmaT/.2/.2;
T1(4)=sigma0(4)^2*dt/sigmaT/.2/.2;
T1(5)=sigma0(5)^2*dt/sigmaT/.2/.2;

T1(1)
T1(2)
T1(3)
T1(4)
T1(5)
T1T=T1(1)+T1(2)+T1(3)+T1(4)+T1(5);
T1T

for nn=1:Nn
x(nn)=x0;
end
for nn=1:Nn
%    Z(nn)=((nn-1)*dNn-4.0)*1/sqrt(Tt);
end
for tt=1:Tt
for nn=1:Nn
Z(nn,tt)=((nn-1)*dNn-4.0)*1/sqrt(T1(tt));

end
end

%ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1/sqrt(Tt));
%ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1/sqrt(Tt));
%for nn=2:Nn-1
%    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1/sqrt(Tt))-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1/sqrt(Tt));
%end
%ZProb contains probaility mass in each cell.
% Pre-compute the hermite polynomilas associated with each subdivision of
% normal.
for tt=1:Tt
for nn=1:Nn
HermiteP(1,nn,tt)=1;
HermiteP(2,nn,tt)=Z(nn,tt);
HermiteP(3,nn,tt)=(Z(nn,tt)^2-1);
HermiteP(4,nn,tt)=Z(nn,tt)^3-3*Z(nn,tt);
HermiteP(5,nn,tt)=Z(nn,tt)^4-6*Z(nn,tt)^2+3;

end
end
%These sigma11, mu11, mu22, sigma22 can be lumped with Y(l1,l2,l3,l4) but I
%kept them separate for a general program where these values could change
%on every step while exponents on x in each term would remain same in
%pre-calcualted Y(l1,l2,l3,l4)
sigma11(1:Order4+1)=0;
mu11(1:Order4+1)=0;
mu22(1:Order4+1)=0;
sigma22(1:Order4+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:(Order4+1)
for tt=1:Tt
if sigma0~=0
sigma11(k,tt)=sigma0(tt)^(k-1);
end
if sigma0~=0
sigma22(k,tt)=sigma0(tt)^(2*(k-1));
end
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

%sigma11(2)=vol1;
%sigma11(3)=vol2;
end
Ft(1:Tt+1,1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;
Fp(1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (Order4)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
Ft(l1,l2,l3,l4) = dt^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
Y = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
Y(1,1,1,1)=0;
for tt=1:Tt
for nn=1:Nn

Y(1,1,1,2)*x(nn)^Fp(1,1,1,2);
x1(nn)=x(nn);
for k = 0 : Order4
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;

x1(nn)=x1(nn) + Y(l1,l2,l3,l4) * ...
mu11(l1)*mu22(l2)*sigma22(l3,tt)*sigma11(l4,tt)* ...
x(nn)^Fp(l1,l2,l3,l4) * HermiteP(l4,nn,tt) * Ft(l1,l2,l3,l4);

end
end
end
end
x(nn)=x1(nn);
x(nn);
end
end

for nn=1:Nn
B(nn)=((nn-1)*dNn-4.0);
end
for nn=2:Nn-1
Dfx(nn) = (x(nn + 1) - x(nn - 1))/(B(nn + 1) - B(nn - 1));
end
fx(1:Nn)=0;
for nn = 2:Nn-1
fx(nn) = normpdf(B(nn),0, 1)/abs(Dfx(nn));
end

%x
%Ft(1,1,1,2)
%Calculate the density with monte carlo below.
%Below the monte carlo sigma00 coefficients related to time dependent sigma0. Only that there are
%eight times smaller and more intervals in monte carlo.
sigma00(1:8)=sigma0(1);
sigma00(9:16)=sigma0(2);
sigma00(17:24)=sigma0(3);
sigma00(25:32)=sigma0(4);
sigma00(33:40)=sigma0(5);
dt0=dt/8;
paths=100000;
X(1:paths)=x0;
Random1(1:paths)=0;
for tt=1:Tt*8
Random1=randn(size(Random1));
X(1:paths)=X(1:paths)+ mu1* X(1:paths).^beta1 *dt0 + ...
mu2 * X(1:paths).^beta2*dt0 + ...
sigma00(tt) * X(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu1^2 * beta1* X(1:paths).^(2*beta1-1) *  dt0^2/2.0 + ...
mu2^2 * beta2 * X(1:paths).^(2*beta2-1) * dt0^2/2 + ...
mu1*mu2* beta1* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
mu1*mu2* beta2* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
sigma00(tt)^2 * gamma * X(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu1 * sigma00(tt) *gamma* X(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
mu2 * sigma00(tt) *gamma* X(1:paths).^(beta2+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu1 *beta1* X(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu2.*beta2* X(1:paths).^(beta2+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths);
end
sum(X(:))/paths   %monte carlo average
%sum(x(:).*ZProb(:))  % Ito-Taylor density average(good enough but slightly approximate)
theta+(x0-theta)*exp(-kappa*dt*Tt)

BinSize=.0025;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff);
plot(x(2:Nn-1),fx(2:Nn-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
end Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

Here is the output graph comparison of the previous toy program that shows monte carlo density and the analytical density with the term structure described in the previous post that is a lognormal with drift SDE with term structure of variance given as
sigma0(1)=.4;
sigma0(2)=.75;
sigma0(3)=.2;
sigma0(4)=.05;
sigma0(5)=.135;
The toy program  has parameters
dX(t)=mu0 *X(t) + sigma0(t)*X(t)*dz(t)
x0=.25;
mu=.25;
T=1; =(dt*Tt)
dt=.2;
Tt=5;.  Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

Sorry this is a slightly modified program to make the meanings of yesterday's program more clear. It could still be somewhat unclear in yesterday's program how to generalize everything. I have slightly modified that program to make it more clear.
function [] = ItoTaylorDensitySimulation010()

% In the expansion of X^alpha
% With SDE
% dX(t)=mu1 * X(t)^beta1 dt + sigma * X(t)^gamma dz(t)
% The first order expansion of X^alpha with the above SDE is
% I have shown the solution with volatility term structure with simple
% lognormal distribution in this code example.
Order4=4;  %expansion order
Nn=161;    % No of normal density subdivisions
dNn=.05;   % Normal density width. would change with number of subdivisions
dt=.20;   % Simulation time interval.
Tt=10;      % No of simulation intervals. Terminal time= Tt*dt;
x0=.250;   % starting value of SDE

alpha=1;   % x^alpha is being expanded
beta1=1;   % the first drift term power.
beta2=0;   % Second drift term power. zero for this examole.
gamma=1.0;   % volatility power.
mu1=.25;   % first drift coefficient.
mu2=0;    % Second drift coefficient. is zero in this simple program.

%Below is the monte carlo volatility term structure for five periods in a
%year. When you change number of density simulation intervals, you will
%have to change this part.
sigma0(1)=.4;
sigma0(2)=.75;
sigma0(3)=.2;
sigma0(4)=.05;
sigma0(5)=.135;
sigma0(6)=.135;
sigma0(7)=.135;
sigma0(8)=.135;
sigma0(9)=.135;
sigma0(10)=.135;
sigmaT=sigma0(1)^2*dt+sigma0(2)^2*dt+sigma0(3).^2*dt+sigma0(4).^2*dt+sigma0(5).^2*dt + ...
sigma0(6)^2*dt+sigma0(7)^2*dt+sigma0(8).^2*dt+sigma0(9).^2*dt+sigma0(10).^2*dt;
%sigmaT
%1/sigmaT
T1(1)=sigma0(1)^2*dt/sigmaT/(1/Tt)^2;
T1(2)=sigma0(2)^2*dt/sigmaT/(1/Tt)^2;
T1(3)=sigma0(3)^2*dt/sigmaT/(1/Tt)^2;
T1(4)=sigma0(4)^2*dt/sigmaT/(1/Tt)^2;
T1(5)=sigma0(5)^2*dt/sigmaT/(1/Tt)^2;
T1(6)=sigma0(6)^2*dt/sigmaT/(1/Tt)^2;
T1(7)=sigma0(7)^2*dt/sigmaT/(1/Tt)^2;
T1(8)=sigma0(8)^2*dt/sigmaT/(1/Tt)^2;
T1(9)=sigma0(9)^2*dt/sigmaT/(1/Tt)^2;
T1(10)=sigma0(10)^2*dt/sigmaT/(1/Tt)^2;

T1(1)
T1(2)
T1(3)
T1(4)
T1(5)
T1(6)
T1(7)
T1(8)
T1(9)
T1(10)
T1T=T1(1)+T1(2)+T1(3)+T1(4)+T1(5)+T1(6)+T1(7)+T1(8)+T1(9)+T1(10);
T1T

for nn=1:Nn
x(nn)=x0;
end
for nn=1:Nn
%    Z(nn)=((nn-1)*dNn-4.0)*1/sqrt(Tt);
end
for tt=1:Tt
for nn=1:Nn
Z(nn,tt)=((nn-1)*dNn-4.0)*1/sqrt(T1(tt));

end
end

%ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1/sqrt(Tt));
%ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1/sqrt(Tt));
%for nn=2:Nn-1
%    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1/sqrt(Tt))-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1/sqrt(Tt));
%end
%ZProb contains probaility mass in each cell.
% Pre-compute the hermite polynomilas associated with each subdivision of
% normal.
for tt=1:Tt
for nn=1:Nn
HermiteP(1,nn,tt)=1;
HermiteP(2,nn,tt)=Z(nn,tt);
HermiteP(3,nn,tt)=(Z(nn,tt)^2-1);
HermiteP(4,nn,tt)=Z(nn,tt)^3-3*Z(nn,tt);
HermiteP(5,nn,tt)=Z(nn,tt)^4-6*Z(nn,tt)^2+3;

end
end
%These sigma11, mu11, mu22, sigma22 can be lumped with Y(l1,l2,l3,l4) but I
%kept them separate for a general program where these values could change
%on every step while exponents on x in each term would remain same in
%pre-calcualted Y(l1,l2,l3,l4)
sigma11(1:Order4+1)=0;
mu11(1:Order4+1)=0;
mu22(1:Order4+1)=0;
sigma22(1:Order4+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:(Order4+1)
for tt=1:Tt
if sigma0~=0
sigma11(k,tt)=sigma0(tt)^(k-1);
end
if sigma0~=0
sigma22(k,tt)=sigma0(tt)^(2*(k-1));
end
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

%sigma11(2)=vol1;
%sigma11(3)=vol2;
end
Ft(1:Tt+1,1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;
Fp(1:(Order4+1),1:(Order4+1),1:(Order4+1),1:(Order4+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (Order4)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
Ft(l1,l2,l3,l4) = dt^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
Y = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
Y(1,1,1,1)=0;
for tt=1:Tt
for nn=1:Nn

Y(1,1,1,2)*x(nn)^Fp(1,1,1,2);
x1(nn)=x(nn);
for k = 0 : Order4
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;

x1(nn)=x1(nn) + Y(l1,l2,l3,l4) * ...
mu11(l1)*mu22(l2)*sigma22(l3,tt)*sigma11(l4,tt)* ...
x(nn)^Fp(l1,l2,l3,l4) * HermiteP(l4,nn,tt) * Ft(l1,l2,l3,l4);

end
end
end
end
x(nn)=x1(nn);
x(nn);
end
end

for nn=1:Nn
B(nn)=1*((nn-1)*dNn-4.0);
end
for nn=2:Nn-1
Dfx(nn) = (x(nn + 1) - x(nn - 1))/(B(nn + 1) - B(nn - 1));
end
fx(1:Nn)=0;
for nn = 2:Nn-1
fx(nn) = normpdf(B(nn),0, 1)/abs(Dfx(nn));
end

%x
%Ft(1,1,1,2)
%Calculate the density with monte carlo below.
%Below the monte carlo sigma00 coefficients related to time dependent sigma0. Only that there are
%eight times smaller and more intervals in monte carlo.
sigma00(1:8)=sigma0(1);
sigma00(9:16)=sigma0(2);
sigma00(17:24)=sigma0(3);
sigma00(25:32)=sigma0(4);
sigma00(33:40)=sigma0(5);
sigma00(41:48)=sigma0(6);
sigma00(49:56)=sigma0(7);
sigma00(57:64)=sigma0(8);
sigma00(65:72)=sigma0(9);
sigma00(73:80)=sigma0(10);
dt0=dt/8;
paths=100000;
X(1:paths)=x0;
Random1(1:paths)=0;
for tt=1:Tt*8
Random1=randn(size(Random1));
X(1:paths)=X(1:paths)+ mu1* X(1:paths).^beta1 *dt0 + ...
mu2 * X(1:paths).^beta2*dt0 + ...
sigma00(tt) * X(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu1^2 * beta1* X(1:paths).^(2*beta1-1) *  dt0^2/2.0 + ...
mu2^2 * beta2 * X(1:paths).^(2*beta2-1) * dt0^2/2 + ...
mu1*mu2* beta1* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
mu1*mu2* beta2* X(1:paths).^(beta1+beta2-1) *dt0^2/2.0 + ...
sigma00(tt)^2 * gamma * X(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu1 * sigma00(tt) *gamma* X(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
mu2 * sigma00(tt) *gamma* X(1:paths).^(beta2+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu1 *beta1* X(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu2.*beta2* X(1:paths).^(beta2+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths);
end
sum(X(:))/paths   %monte carlo average
%sum(x(:).*ZProb(:))  % Ito-Taylor density average(good enough but slightly approximate)
%theta+(x0-theta)*exp(-kappa*dt*Tt)

BinSize=.0025;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff);
plot(x(2:Nn-1),fx(2:Nn-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
end Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

Dear friends, I believe I have nailed the whole thing down now. I will be following with a detailed program in 2-3 days but I will share the details of the algorithm now. We do not need to know the term structure of variance like I did in the previous program. We just need to know the "existing variance" of the SDE in evolution and the "simulation variance" of the next increment to be simulated. And the $\sqrt(\frac{1}{T_0})$ and $\sqrt(\frac{1}{T_1})$ for simulation would be assigned in the ratio that $T_0$ and $T_1$ would be in the ratio of existing variance and simulation variance. The existing variance would be calculated by a process similar to change of derivative for probability distributions. What I mean is that for a fixed probability mass a smaller interval would mean a smaller variance and a larger interval would mean a larger variance and this ratio would be calculated with respect to probability mass in the density subdivision. In the program, I would define this exactly and clearly show how to calculate the existing variance in the metric/currency of simulation variance so that both could be added and then a proportional simulation interval $T_1$ to $T_0$ for existing variance interval could be defined in proportions of existing and simulation variance. And on each step of evolution, our probability density would exactly match the true probability density. The whole process is very simple and I would be presenting a comprehensive program for this in next 2-3 days. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

The existing variance would be calculated by a process similar to change of derivative for probability distributions. What I mean is that for a fixed probability mass a smaller interval would mean a smaller variance and a larger interval would mean a larger variance and this ratio would be calculated with respect to probability mass in the density subdivision. In the program, I would define this exactly and clearly show how to calculate the existing variance in the metric/currency of simulation variance so that both could be added and then a proportional simulation interval $T_1$ to $T_0$ for existing variance interval could be defined in proportions of existing and simulation variance.
As friends can see that comparison of any density can be made with normal on a standard deviation to standard deviation basis. Change of derivative comparison can be done for the peak of normal or tail of normal on a standard deviation to standard deviation basis. change of derivative for the peak can be compared to peak of another normal/lognormal/some other CEV exponent as long as we are comparing the same standard deviation subdivisions on both densities. I will explain it in detail in the next program. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

We can very easily find the "existing variance" associated with the density of most of the SDEs that take $X^{\gamma} dz(t)$ in variance trem with any possibly non-linear drift combination. when these SDE densities are transformed to $X^{(1-\gamma)}$ they become normal densities and then variance can be directly and trivially found by using a ratio with standard normal density. We can very easily find the variance of density of CEV type SDE with non-linear drift in normal reference variance term just by calculating the values of $X_h^{(1-\gamma)}-X_l^{(1-\gamma)}$ where $X_l$ and $X_h$ are left and right side of the of the density grid whose existing density has to be calculated. This is very remarkable since the density can get so many types of non-linear drifts but since transformed density is normal, its variance can still be calculated in normal variance terms and again very interestingly since the transformed density is normal, the variance of every subdivision/fraction on the density simulated with any single unique SDE turns out to be exactly equal as the  relative scaling of any fraction on a normal density is exactly the same. This is very very interesting and I will be coming back with a program using this in next two days. Amin
Topic Author
Posts: 2100
Joined: July 14th, 2002, 3:00 am

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

This program simulates the density by making a comparison of existing variance of the density in place and the variance to be simulated. This program simulates the density of a single unique SDE that takes arbitrary exponent in drift and a CEV type diffusion term. Can easily be used for densities that take more than one term in drift. I have just used first two hermite polynomials in the simulation of this program with an obviously shorter time step size with just two hermite polynomials and you can easily extend to greater number of hermite polynomials.
function [] = ItoTaylorNewWilmottQuadratic03()
%This program simultes the density by making a comparison of
%existing variance of the density and variance to be simulated.
%The first time step is directly simulated
%But all later steps calculate a ratio of existing variance SigmaExisting
%and the next step terminal variance SigmaOneSTepTermianl.
%This terminal varaince is calculated by advancing the density using a
%normal.
%And then termianl variance is calculated.
% And the T which is varaince of normal in next step hermite polynomials is calculated as
%T=(1+ (Ratio_SigmaT_SigmaE)^2)/(Ratio_SigmaT_SigmaE-1)^2*sqrt(2)
%Using this scaling for hermite polynomials our simulated density remains
%equal to the true density on each simulation step.
Nn=161;    % No of normal density subdivisions
dNn=.05;   % Normal density width. would change with number of subdivisions
Tt=10;      % simulation level. Terminal time= Tt*dt;
x0=1.0;   % starting value of SDE
X(1:Nn)=x0; %Initialize the density subdivisions with initial value.
beta1=.75;   % the first drift term power.
gamma=.650;   % volatility power.
mu0=.5;   % first drift coefficient.
dt(1)=.05;
dt(2)=.05;
dt(3)=.05;
dt(4:10)=.05
sigma0(1)=.85;
sigma0(2)=.65;
sigma0(3)=.75;
sigma0(4)=.85;
sigma0(5)=.85;
sigma0(6)=.95;
sigma0(7)=.55;
sigma0(8)=.65;
sigma0(9)=.75;
sigma0(10)=.95;
T1=1.0;%Hermite polynomials have unit varaince for first step.
for tt=1:Tt
if(tt>1)
tt
for nn=2:Nn-1
Zn(nn)=((nn-1)*dNn-4.0);
%            %This is a training step that advances the density one step using a
% hermite polynomial with standard variance. This is not
%right simulation and only for the purpose of
%calculation of variance.
X1(nn)=X(nn)+ ...
mu0* X(nn).^beta1 * dt(tt) + ...
mu0^2 * beta1* X(nn).^(2*beta1-1) * dt(tt).^2/2.0 + ...
(sigma0(tt) * X(nn).^gamma *sqrt(dt(tt)) + ...
sigma0(tt) * mu0 * beta1 * X(nn).^(beta1 + gamma -1) *(1 - 1 / sqrt(3)) * dt(tt).^1.5 + ...
sigma0(tt) * mu0 * gamma * X(nn).^(beta1 + gamma -1) * 1 / sqrt(3) * dt(tt)^1.5) * Zn(nn) + ...
sigma0(tt).^2 *gamma* X(nn).^(2*gamma -1 ) .* dt(tt)/2 .* ( Zn(nn).^2 - 1 );
end

%Calculate the start and ends of mid density subdivisions
%y0h and y0l are start and end of mid of existing density
%before the simulation
y0h=X(81)+(X(82)-X(81))/2.0;
y0l=X(81)-(X(81)-X(80))/2.0;

%y1h and y1l are start and end of mid of density that has been
%simulated to calculate the variance after the simulation
y1h=X1(81)+(X1(82)-X1(81))/2.0;
y1l=X1(81)-(X1(81)-X1(80))/2.0;

% We calculate the standard deviation ratio after converting both
% densities to normal density after transformation.
% For a single unique SDE, we need to compare just one density subdivison
% and we just compare the ratio of varaince of middle subdivisions
% of existing and simulated training density
Ratio_SigmaExisting_sigmaOneStepTerminal=((y1h^(1-gamma)-y1l^(1-gamma))/(y0h^(1-gamma)-y0l^(1-gamma)));

%This is the varaince scaling of normal in the hermite polynomials to be
%used in the next step.
T1=(1+Ratio_SigmaExisting_sigmaOneStepTerminal^2)/(Ratio_SigmaExisting_sigmaOneStepTerminal-1)^2*sqrt(2);

end
for nn=1:Nn
Z(nn,tt)=((nn-1)*dNn-4.0)*1/sqrt(T1);

HermiteP(1,nn,tt)=1;
HermiteP(2,nn,tt)=Z(nn,tt);
HermiteP(3,nn,tt)=(Z(nn,tt)^2-1);
HermiteP(4,nn,tt)=Z(nn,tt)^3-3*Z(nn,tt);
HermiteP(5,nn,tt)=Z(nn,tt)^4-6*Z(nn,tt)^2+3;

end
for nn=1:Nn
X(nn)=X(nn)+ ...
mu0* X(nn).^beta1 * dt(tt) + ...
mu0^2 * beta1* X(nn).^(2*beta1-1) * dt(tt)^2/2.0 + ...
(sigma0(tt) * X(nn)^gamma *sqrt(dt(tt)) + ...
sigma0(tt) * mu0 * beta1 * X(nn).^(beta1 + gamma -1) *(1 - 1 / sqrt(3)) * dt(tt)^1.5 + ...
sigma0(tt) * mu0 * gamma * X(nn).^(beta1 + gamma -1) * 1 / sqrt(3) * dt(tt)^1.5) * HermiteP(2,nn,tt) + ...
sigma0(tt)^2 *gamma* X(nn).^(2*gamma -1 ) * dt(tt)/2 * HermiteP(3,nn,tt);
end

end
for nn=1:Nn
B(nn)=((nn-1)*dNn-4.0);
end
for nn=2:Nn-1
DfX(nn) = (X(nn + 1) - X(nn - 1))/(B(nn + 1) - B(nn - 1));
end
fX(1:Nn)=0;
for nn = 2:Nn-1
fX(nn) = normpdf(B(nn),0, 1)/abs(DfX(nn));
end
%Monte carlo simulation follows for comparison of our simulated desnsity and
% monte carlo density.
paths=100000;
X2(1:paths)=x0;
Random1(1:paths)=0;
sigma00(1:Tt)=sigma0(1:Tt);
for tt=1:Tt
dt0=dt(tt);
Random1=randn(size(Random1));
X2(1:paths)=X2(1:paths)+ mu0* X2(1:paths).^beta1 *dt0 + ...
sigma00(tt) * X2(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu0^2 * beta1* X2(1:paths).^(2*beta1-1) *  dt0^2/2.0 + ...
sigma00(tt)^2 * gamma * X2(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu0 * sigma00(tt) *gamma* X2(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu0 *beta1* X2(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths);
end
sum(X2(:))/paths   %monte carlo average
BinSize=.005;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X2,paths,BinSize,MaxCutOff);
plot(X(3:Nn-2),fX(3:Nn-2),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
end


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

Xmin=0;
Xmax=0;
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
XDensity(1:IndexMax)=0.0;
for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);
if(real(index)<1)
index=1;
end
if(real(index)>IndexMax)
index=IndexMax;
end
XDensity(index)=XDensity(index)+1.0/Paths/BinSize;
end
IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;
end  