SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

July 2nd, 2018, 10:56 am

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);
%YdtCoeff= ItoTaylorCoeffsNew(alphadt,beta1,beta2,gamma);
%YdzCoeff= ItoTaylorCoeffsNew(alphadz,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
 
User avatar
Amin
Topic Author
Posts: 1969
Joined: July 14th, 2002, 3:00 am

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

July 11th, 2018, 9:49 am

Amin wrote:
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.  
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...