SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

October 12th, 2018, 1:19 pm

Here are some interesting things that have come up in my research. suppose there is an SDE of the form [$]dX(t)=\mu(t) X^{\beta}(t) dt + \sigma(t) X^{\gamma}(t) dz(t)[$] that we want to simulate with this new Ito-Taylor density simulation method. We noticed that a lot of SDEs of this form could easily be simulated with this method after transformation into Bessel coordinates. While some non-linear SDEs like stochastic volatility SDEs could not be simulated in a satisfactory manner. Here I try to explain this anomaly and I will like to state a few observations.
When we transformed SDEs of the type given above in SDE equation into Bessel coordinates, the numerically calculated first derivative of transformed density with respect to normal variable[$]\frac{\partial Y}{\partial Z}[$] when calculated from density mass equality, remains almost constant and does not take a significant second derivative. Here Y represents Bessel transformed variable from SDE variable X. For all those SDEs that become linear after converting to Bessel coordinates, we can use several methods including the method in my original program. Here the criteria for linearity is that [$]\frac{\partial Y}{\partial Z}[$] is almost constant and [$]\frac{{\partial}^{2} Y}{\partial {Z}^2}[$] is almost negligible. In the equations that are non-linear like stochastic volatility equation that takes a non-negligible second derivative with normal variable, we need to use this numerically calculated second or higher derivatives with possibly hermite polynomials to appropriately calculate existing variance and innovations(I earlier called it simulation variance) variance and their ratios. We cannot easily use the original coordinates since in the original coordinates, the SDE variable continues to take many higher derivatives with normal. I am sure with a slightly better understanding of the method, we could also use original coordinates for density simulation as well. again the reason for my being able to successfully simulate the densities of many simpler linear processes was that when I converted them into transformed coordinates, all higher derivatives with respect to normal variable disappeared and only the first derivative remained. In general, after the appropriate transformations, the number of significant higher derivatives with respect to normal variable sharply decreases. I believe that these numerically calculated higher derivatives of transformed densities with normal variable have to be used with variance of hermite polynomials to satisfactorily calculate the existing and innovations variance.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

October 14th, 2018, 5:30 pm

Suppose we have an SDE for X and we have calculated its existing density at a certain time. We suppose that we have divided the density of X into subdivisions that have the same mass as linear subdivisions of the normal density. We want to calculate local variance for each subdivision. In order to do that we numerically calculate the derivatives  [$]\frac{\partial X}{\partial Z}[$] and [$]\frac{\partial^2 X}{\partial Z^2}[$] where Z is a standard normal variable and these derivatives are calculated with probability mass equality for each subdivision of X. For now we limit ourselves to first two derivatives. We can expand X as a function of Z around the center of nth subdivision, [$]x_n[$] as 
[$]X(Z)=x_n + \frac{\partial X}{\partial Z} Z + .5 \frac{\partial^2 X}{\partial Z^2} Z^2[$] 
Using the variance for standard normal and its powers, we can easily calculate the local variance of X at the nth subdivision.
The only problem is that Z is perfectly correlated across time in our Ito-Taylor density simulation method while it is totally uncorrelated across time in reality so we have to make an adjustment in the value of Z that we use for simulation of Ito-Taylor density simulation method.
if [$]X_0[$] denotes the existing density variable and [$]X_1[$] denotes the innovations density variable, the true variance that is the right variance at the certain subdivision after taking the innovation simulation step should be
[$](\frac{\partial X_0}{\partial Z})^2 + (\frac{\partial X_1}{\partial Z})^2+.5 (\frac{\partial^2 X_0}{\partial Z^2} )^2  +.5 ( \frac{\partial^2 X_1}{\partial Z^2})^2 [$] 
while Ito-Taylor simulation assumes complete correlation for the Z's across time and their variance without any adjustment turns out to be far greater. This variance can also be very easily calculated.
Again we have to make sure that Ito-Taylor simulation method variance matches with the true variance of the SDE at that particular subdivision. This can also be easily done. This is not much more difficult than what we did for one derivative case and requires slightly altering the analysis on the same lines. And then we will advance the SDE density with different local scaling of hermite polynomials on each subdivision so that local variance simulated with Ito-Taylor density simulation method matches with the true simulation variance.
I hope to post a full-fledged program that calculates the true density with Ito-Taylor density simulation method in a few days.
I am travelling tomorrow so will post more details on Tuesday.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

October 17th, 2018, 12:33 pm

Here is the correlated sum of existing variance and innovation/simulation variance of Ito-Taylor density simulation method in terms of original Z's. We will weight [$]Z_1[$] with [$]\rho[$] in LHS which is the correlated evolution and take expectations and equate with true uncorrelated (and unweighted) sum of existing and innovation/simulation variance on RHS and solve this quartic equation for [$]\rho[$] and this calculated [$]\rho[$] will be used for scaling the Z's in hermite polynomials in Ito-Taylor density simulation method. I am writing the correlated sum of existing and innovation variance with [$]Z_1[$] weighted with [$]\rho[$] in RHS and true variance after one step evolution  in RHS. After taking expectations for Z, we will solve these resulting quartic equation in [$]\rho[$] and it is [$]\rho[$] that will scale the Z's in hermite polynomials in the density simulation method. In the first part of the equation or LHS, both [$]Z_0[$] and [$]Z_1[$] are perfectly correlated when they appear as a product but otherwise they are standard normal with unit variance and [$]\rho[$] is simply a scalar multiplier.

[$](\frac{\partial X_0}{\partial Z_0})^2 E(Z_0 )^2+ (\frac{\partial X_1}{\partial Z_1})^2 E(\rho Z_1)^2 + 2 \frac{\partial X_0}{\partial Z_0} \frac{\partial X_1}{\partial Z_1} E[Z_0 (\rho Z_1)][$]
[$]+ .25 (\frac{\partial^2 X_0}{ \partial {Z_0}^2})^2  E (Z_0)^4-.25 (\frac{\partial^2 X_0}{\partial {Z_0}^2})^2  ({E(Z_0)}^2)^2 [$]
[$] +.25 (\frac{\partial^2 X_1}{ \partial {Z_1}^2})^2  E (\rho Z_1)^4-.25 (\frac{\partial^2 X_1}{\partial {Z_1}^2})^2  ({E(\rho Z_1)}^2)^2 [$]
[$]+.5 \frac{\partial^2 X_0}{\partial {Z_0}^2} \frac{\partial^2 X_1}{\partial {Z_1}^2} E[{Z_0}^2 (\rho Z_1)^2] -.5 \frac{\partial^2 X_0}{\partial {Z_0}^2} \frac{\partial^2 X_1}{\partial {Z_1}^2} (E[ Z_0 (\rho Z_1)])^2[$]
=[$](\frac{\partial X_0}{\partial Z_0})^2 E(Z_0 )^2+ (\frac{\partial X_1}{\partial Z_1})^2 E( Z_1)^2 [$]
[$]+ .25 (\frac{\partial^2 X_0}{ \partial {Z_0}^2})^2  E (Z_0)^4-.25 (\frac{\partial^2 X_0}{\partial {Z_0}^2})^2  ({E(Z_0)}^2)^2 [$]
[$] +.25 (\frac{\partial^2 X_1}{ \partial {Z_1}^2})^2  E ( Z_1)^4-.25 (\frac{\partial^2 X_1}{\partial {Z_1}^2})^2  ({E( Z_1)}^2)^2 [$]
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

October 18th, 2018, 7:15 pm

OK friends, I tried the above with explicit solution of quartic equation and applied it to stochastic volatility type SDEs. My fit was generally remarkably better than before but I could see some minor discrepancy in the right tail of the stoch  vol SDEs. Since in asset pricing, we want to be precise, so more care has to be taken. I am very sure any discrepancy in the right tail in the skew is due to small amount of third derivative that my analysis till 2nd order did not cover. Only problem is that sextic equation may not have a closed form solution. I will redo the above analysis till fourth order and then write a program that does Newton iteration to solve the eighth order polynomial equation for the scaling factor. I Am sure we can get away with around hundred points in the SDE density and then take a quarter year time step and do all newton iterations for the whole 100 point grid to get the solution of eighth order polynomial in parallel. I hope I will complete all of that in 7-10 days and post a new program here after that.
I am hundred percent certain after looking at the results of stoch vol SDEs densities that above analysis works perfectly fine since there was a remarkably marked improvement in results from first order analysis. Anyway, for researchers interested in looking at it on their own, I hope to post a matlab program with just 2nd order solution from analysis in previous post in next 2-3 days. I hope this program will be helpful for many interested friends to do experiments on their own.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

October 25th, 2018, 6:10 pm

This is something interesting so I thought I will share it. I was thinking of the relationship between derivatives of the SDE variable X with respect to Z and the coefficients of hermite polynomials and I set up a small experiment in which I constructed a random variable X from hermite polynomials along a density evolution using Ito-Taylor density evolution method.
if [$]X(t+1)=c_1(X(t)) H_1(Z) + c_2(X(t)) H_2(Z) + c_3(X(t)) H_3(Z) [$]
If we differentiate [$]X(t+1)[$]  with respect to Z using density mass equality, we get that first derivative of X with respect to Z equals the coefficient [$]c_1(X(t))[$] of first hermite polynomial i.e. [$] \frac{dX(t+1)}{dZ} = c_1(X(t))[$] and second derivative follows the relationship  [$].5 \frac{d^2 X(t+1)}{dZ^2}=c_2(X(t))[$] 
and the third derivative followed [$]-\frac{1}{6}  \frac{d^3 X(t+1)}{dZ^3} =c_3(X(t))[$] 
so these weighted derivatives with respect to standard normal density (that are coefficients of Taylor expansion of X with respect to Z) have a direct interpretation as coefficients of hermite polynomials. This should hold when X is a continuous random variable like we have in the density of SDEs. 
To give a concrete idea, 1st derivative of X with respect to Z at nth point on the grid is simply calculated by finite differences as [$]\frac{X(n+1)-X(n-1)}{Z(n+1)-Z(n-1)}[$] on the Ito-Taylor density grid (to 2nd order accuracy).
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

October 27th, 2018, 6:13 pm

sorry for being so wrong as I have been in the last post. But here I re-did my analysis for the relationship between derivatives of X with respect to Z and the coefficients of hermite polynomials.
1st four hermite polynomials are
[$]H_1(Z)=Z[$]
[$]H_2(Z)=Z^2-1[$]
[$]H_3(Z)=Z^3-3Z[$]
[$]H_4(Z)=Z^4-6Z^2+3[$]
We suppose we have a function/variable that is composed of first four hermite polynomials and given as
[$]X(Z)=c_1 H_1(Z) + c_2 H_2(Z)+c_3 H_3(Z) + c_4 H_4(Z)[$]
differentiating the above appropriately, we get
[$]\frac{\partial^4 X}{\partial Z^4}=24 c_4[$]
[$]\frac{\partial^3 X}{\partial Z^3}=6 c_3+ 24 c_4 Z[$]
[$]\frac{\partial^2 X}{\partial Z^2}=2 c_2+6 c_3 Z + c_4(12 Z^2 -12)[$]
[$]\frac{\partial X}{\partial Z}=c_1 + 2 c_2 Z + c_3 (3 Z^2-3) + c_4 (4 Z^3-12 Z)[$]

Above equations can easily be solved but at Z=0, these equations are particularly easy to solve and we get
[$]\frac{\partial^4 X}{\partial Z^4}=24 c_4[$]
[$]\frac{\partial^3 X}{\partial Z^3}=6 c_3[$]
[$]\frac{\partial^2 X}{\partial Z^2}=2 c_2-12 c_4[$]
[$]\frac{\partial X}{\partial Z}=c_1 -3 c_3[$]

So we are making a slow but steady progress.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

November 3rd, 2018, 8:22 pm

On this weekend I will post a new program for simulation. Currently, it has small instabilities but I am sure numerical analysis experts and many other intelligent people could easily overcome them. The logic is very simple: I add drift in a linear fashion and volatility term in an orthogonal fashion (square root of squared sum with previous increment from the mean). It turns out to be very simple and there is just very basic mathematics in it. Though there might still be cases with problems but I hope this will be very helpful. And again there are minor instabilities that grow with time and some way still have to be worked out how to tackle them but I still decided to post the program. I will be posting it tomorrow or on Sunday after writing explanations on the program.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

November 4th, 2018, 6:34 pm

Here is the new matlab program.In this program, I am simulating the SDE given as
       dy(t)=kappa*(theta-y(t))*dt + y(t)^gamma * dz(t);
I have not directly simulated the SDE but simulated the transformed Bessel process version of the SDE and then changed coordinates to retrieve the SDE in original coordinates. The present program will directly simulate only the Bessel Process version of the SDE in transformed coordinates.
The program has instability in 3-4 subdivisions in the mid of the density and I will be releasing a program in another few days taking care of that and fixing the problem. You can do it on your own as well.
In this program, We are not taking any cross terms for the drift and volatility and have just added one constant volatility term that occurs in the Bessel transformed SDE. The only volatility term is sigma0 *sqrt(t)*H1(N).
Here is the program. Sorry that I have not optimized the code in any way and this is just to demonstrate the idea of the program. Will post a complete worked out program in a few days with complete optimizations.

function [] = ItoTaylorDensityWilmottSimplev1()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=kappa(theta-y(t))*dt + y(t)^gamma dz(t);
%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.
%The program has instability in 3-4 subdivisions in the mid of the density
%and I will be relasing a program in another few days taking care of that
%and fixing the problem. You can do it on your own as well.
dt=.125/2;   % Simulation time interval.
Tt=8;     % Number of simulation levels. Terminal time= Tt*dt;  
OrderA=4;  %expansion order for the analytic SDE.
OrderM=4;  %expansion order for the monte carlo simulation.
dNn=.05/1;   % Normal density width. would change with number of subdivisions
Nn=161;  % No of normal density subdivisions
NnMidn=81;%Median cell/subdivision.
NnMid=4.0;
x0=.30;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.950;   % volatility power. Before you change this to values where 
            %diffusion goes to zero, you will have to change the 
            %calculation of the mean of the density in the body of the program
            %do not decrease it more than .75 before changing that.
kappa=1.0;   %mean reversion parameter.
theta=.2;   %mean reversion target
mu1=theta*kappa;   % first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=1.0;%Volatility value
alpha=1.0;% x^alpha is being expanded
alpha1=1-gamma;%For the expansion in transformed coordinates.

for nn=1:Nn
    x(nn)=x0^(1-gamma)/(1-gamma); %Initial value is assigned to each 
                                  %subdivision in transformed coordinates. 
end
%Please read notes at the end of the program for understanding of what is
%Y(l1,l2,l3,l4) and what are l1,l2, l3 and l4. In the ItoTaylor expansion 
%l1 is first drift term expansion order index, l2 is for the second drift
%term, l3 is for quadratic variation drift term index and l4 is the
%volatility expansion order index.
%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:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:Tt+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = 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));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
            end
        end
    end
end

for nn=1:Nn
    Z(nn)=((nn-1)*dNn-NnMid);
end

ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%The routine below calls the Ito-Taylor Coeffs in transformed coordinates
%and uses alpha1 as the power where alpha1=1-gamma.
Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
Y(1,1,1,1)=0;
for nn=1:Nn
    Z(nn)=((nn-1)*dNn-NnMid);%*Rho1(nn);
    HermiteP(1,nn)=0;
    HermiteP(2,nn)=Z(nn);%rho1(nn)
    HermiteP(3,nn)=(Z(nn)^2-1);%rho2(nn)
    HermiteP(4,nn)=(Z(nn)^3-3*Z(nn));
    HermiteP(5,nn)=(Z(nn)^4-6*Z(nn)^2+3);
end
for tt=1:Tt
    for nn=1:Nn
        x1(nn)=x(nn);
        x_mean(nn)=0;
        for k = 0 : OrderA
            for m = 0:k
                l4 = k - m + 1;
                for n = 0 : m
                    l3 = m - n + 1;
                    for j = 0:n
                        l2 = n - j + 1;
                        l1 = j + 1;
                        if(l4<=1)
                        x_mean(nn)=x_mean(nn) + YqCoeff(l1,l2,l3,l4) .* ...
                            mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
                            ((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* Ft(l1,l2,l3,l4);
                        %x_mean calculates the total drift from combination
                        %of three drift terms in their Ito-Taylor
                        %expansion.
                        %((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) is
                        %term for the original untransformed SDE variable.
                        %This is basically ODE for each subdivision where
                        %initial value of ODE is reset at each time level
                        %after addition of diffusion increment fro previous
                        %time levelthat is different for each subdivision.
                        end
                       
                        if((l1<=1)&&(l2<=1)&&(l3<=1)&&(l4==2))
                        
                        
                        x1(nn)=x1(nn) + YqCoeff(l1,l2,l3,l4) .* ...
                            mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
                            ((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* HermiteP(l4,nn) .* Ft(l1,l2,l3,l4);
                        %We are not taking any cross terms for the drift
                        %and volatility and have just added one constant
                        %volatility terms that occurs in the Bessel
                        %transformed SDE. l1=1 means zeroth expansion
                        %order. So we have zero expansion order for all
                        %drifts. I could have written the above volatility 
                        %term only as sigma0 *sqrt(t)*HermiteP(2,nn)but
                        %wrote this long expression due to legacy code.
                        %You can unquote the term below instead of the
                        %above term.
                        %x1(nn)=x1(nn) +sigma0*HermiteP(l4,nn)*sqrt(dt);    
                        end
                    end
                end
            end
        end
    dx(nn)=x1(nn)-x(nn); %diffusion increment.
    end
    if(tt>1)  %For first step when starting from delta origin, we do not 
              %want squared orthogonal addition for volatility. 
        MeanPrev=sum(ZProb(1:Nn).*x(1:Nn));%Mean of the ODE from previous step
        %Please be careful, when working with diffusion exponents that go
        %negative, this simple calculation of the mean will corrupt the
        %whole program when all sort of negative and imaginary values are
        %used to calculate the mean. For that case, you will have to
        %carefully take the mean with only values that are on positive real
        %axis.
        for nn=1:Nn
            if(x(nn)>MeanPrev)    
                x(nn)=MeanPrev+x_mean(nn)+sqrt((x(nn)-MeanPrev).^2+dx(nn)^2);
                %Add the orthogonal increment in a squared fashion if above
                %the mean and add. x_mean is drift increment and is linearly added.
            elseif (x(nn)<MeanPrev)    
                x(nn)=MeanPrev+x_mean(nn)-sqrt((x(nn)-MeanPrev).^2+dx(nn)^2);
                %Add the orthogonal increment to already negative value and
                %subtract.
            end
        end
    else
        x(1:Nn)=x(1:Nn)+x_mean(1:Nn)+dx(1:Nn); %This is linear addition in 
                                %first time step when starting from delta origin. 
    end
end
x(NnMidn)=(x(NnMidn+1)+x(NnMidn-1))/2;  %I just made this very quick and bad fix for the 
                                %unstable value in the mid. Mid three or
                                %four subdivisions have instability.
for nn=1:Nn
    y_x(nn) = ((1-gamma)*x(nn)).^(1/(1-gamma)); %Convert to original coordinates
end
for nn=2:Nn-1
    Dfx(nn) = (x(nn + 1) - x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfy_x(nn) = (y_x(nn + 1) - y_x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end

fx(1:Nn)=0;
fy_x(1:Nn)=0;
for nn = 2:Nn-1
    fx(nn) = normpdf(Z(nn),0, 1)/abs(Dfx(nn)); %Bessel process density
    fy_x(nn) = normpdf(Z(nn),0, 1)/abs(Dfy_x(nn));%Origianl coordinates density
end
str=input('Analytic Values Calculated. Press a key to start Monte Carlo');
%Calculate the density with monte carlo below.
rng(85109634, 'twister')
paths=100000;
X(1:paths)=x0^(1-gamma)/(1-gamma); %Bessel process Monte carlo
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:Tt
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
    %HermiteP1(6,1:paths)=Random1(1:paths).^5-10*Random1(1:paths).^3+15*Random1(1:paths);
    %HermiteP1(7,1:paths)=Random1(1:paths).^6-15*Random1(1:paths).^4+45*Random1(1:paths).^2-15;
    XX0=X;
    YY0=YY;
    for k = 0 : OrderM
        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;
                    if(l4<=5)            
                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
                    mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                    ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
                    HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                    
                %Above is bessel process monte carlo equation.
                
                    YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                    mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                    YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                %Above is original coordinates monte carlo equation.
                    end
                end
            end
        end
    end
end

X(X<0)=0;
sum(X(:))/paths   %monte carlo average for bessel coordinates
sum(x(2:Nn-1).*ZProb(2:Nn-1))  % Ito-Taylor density average for bessel coordinates
sum(y_x(2:Nn-1).*ZProb(2:Nn-1)) %Original process average from coordinates 
                                %derived from bessel process 
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075;
MaxCutOff=30;
[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');
str=input('Look at Bessel graphs comparison between monte carlo density and analytic density');
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_x(2:Nn-1),fy_x(2:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('Look at Original variable graphs comparison between monte carlo and analytic density');
end

% 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. 

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

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

November 4th, 2018, 6:55 pm

This is how you will see the final graph comparison of the stochastic volatility SDE with monte carlo in green and analytic density in red. Parameters are
x0=.3, theta=.2,kappa=1.0, gamma=.95,sigma=1.0; and T=.5 years. You can see the slight instability in the middle that affects middle 3-4 subdivisions.

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

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

November 4th, 2018, 7:48 pm

Sorry, I missed many supporting matlab files for the program in post 773. I am posting them here.
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 second supporting file.
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/(2);
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
And here is the graphing from monte carlo file.
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(X(p)<0)
X(p)=0;
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: 2085
Joined: July 14th, 2002, 3:00 am

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

November 7th, 2018, 11:07 am

I am posting a new version of the program which is quite better than the first version. It still suffers from the instability issue but to a lesser degree. I have slightly changed the framework and my purpose for posting new code is that anybody who wants to totally finesse the program and remove any instability will have a far better starting version. One reminder is that you can easily simulate any SDE of the form
dx(t)=mu1 X(t)^beta1 dt + mu2 X(t)^beta2 dt+ sigma dz(t);
I have simulated the Bessel transformed form of the above SDE and the program automatically converts the above SDE to original coordinates after simulating the Bessel process.
Here is the code.
function [] = ItoTaylorDensityWilmottSimplev3()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=kappa(theta-y(t))*dt + y(t)^gamma dz(t);
%I have not directly simulated the SDE but simulated the transformed 
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.
%The program has instability in 3-4 subdivisions in the mid of the density
%and I will be relasing a program in another few days taking care of that
%and fixing the problem. You can do it on your own as well.
dt=.125/2;   % Simulation time interval.
Tt=8;     % Number of simulation levels. Terminal time= Tt*dt;  
OrderA=4;  %expansion order for the analytic SDE.
OrderM=4;  %expansion order for the monte carlo simulation.
dNn=.1/1;   % Normal density width. would change with number of subdivisions
Nn=80;  % No of normal density subdivisions
%NnMidn=81;%Median cell/subdivision.
NnMidl=40;
NnMidh=41;
NnMid=4.0;
x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.950;   % volatility power. Before you change this to values where 
            %diffusion goes to zero, you will have to change the 
            %calculation of the mean of the density in the body of the program
            %do not decrease it more than .75 before changing that.
kappa=.750;   %mean reversion parameter.
theta=.150;   %mean reversion target
mu1=1*theta*kappa;   % first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=1.0;%Volatility value
alpha=1.0;% x^alpha is being expanded
alpha1=1-gamma;%For the expansion in transformed coordinates.

for nn=1:Nn
    x(nn)=x0^(1-gamma)/(1-gamma); %Initial value is assigned to each 
                                  %subdivision in transformed coordinates. 
end
%Please read notes at the end of the program for understanding of what is
%Y(l1,l2,l3,l4) and what are l1,l2, l3 and l4. In the ItoTaylor expansion 
%l1 is first drift term expansion order index, l2 is for the second drift
%term, l3 is for quadratic variation drift term index and l4 is the
%volatility expansion order index.
%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:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:Tt+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = 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));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
            end
        end
    end
end

for nn=1:Nn
    Z(nn)=((nn-.5)*dNn-NnMid);
end
%Z
%str=input('Show Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);
ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%The routine below calls the Ito-Taylor Coeffs in transformed coordinates
%and uses alpha1 as the power where alpha1=1-gamma.
Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
Y(1,1,1,1)=0;
CC(1:Nn)=0;
for nn=1:Nn
    Z(nn)=((nn-.5)*dNn-NnMid);%*Rho1(nn);
    HermiteP(1,nn)=0;
    HermiteP(2,nn)=Z(nn);%rho1(nn)
    HermiteP(3,nn)=(Z(nn)^2-1);%rho2(nn)
    HermiteP(4,nn)=(Z(nn)^3-3*Z(nn));
    HermiteP(5,nn)=(Z(nn)^4-6*Z(nn)^2+3);
end
for tt=1:Tt
    for nn=1:Nn
        x1(nn)=x(nn);
        x_mean(nn)=0;
        for m = 0 : OrderA    
           for n = 0 : m
               l3 = m - n + 1;
               for j = 0:n
                   l2 = n - j + 1;
                   l1 = j + 1;
                   x_mean(nn)=x_mean(nn) + YqCoeff(l1,l2,l3,l4) .* ...
                        mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4).* ...
                        ((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* Ft(l1,l2,l3,l4);
                        %x_mean calculates the total drift from combination
                        %of three drift terms in their Ito-Taylor
                        %expansion.
                        %((1-gamma)*x(nn)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) is
                        %term for the original untransformed SDE variable.
                        %This is basically ODE for each subdivision where
                        %initial value of ODE is reset at each time level
                        %after addition of diffusion increment fro previous
                        %time levelthat is different for each subdivision.
                  %      end
               end    
           end   
        end
        dx(nn)= +sigma0*HermiteP(2,nn)*sqrt(dt);
                        %We are not taking any cross terms for the drift
                        %and volatility and have just added one constant
                        %volatility terms that occurs in the Bessel
                        %transformed SDE.   
    
    end
    if(tt>1)
        
        
        %For first step when starting from delta origin, we do not 
              %want squared orthogonal addition for volatility. 
        MeanPrev=sum(ZProb(1:Nn).*x(1:Nn));%Mean of the ODE from previous step
        x_meanAvg=sum(ZProb(1:Nn).*x_mean(1:Nn));
        %MeanPrev=.5*x(NnMidl)+.5*x(NnMidh);
        %MeanPrev=.5*MeanPrev+.5*MeanPrev0;
        %Please be careful, when working with diffusion exponents that go
        %negative, this simple calculation of the mean will corrupt the
        %whole program when all sort of negative and imaginary values are
        %used to calculate the mean. For that case, you will have to
        %carefully take the mean with only values that are on positive real
        %axis.
        MeanPrev=sum(ZProb(1:Nn).*x(1:Nn));
        
        %MeanPrev
        for nn=1:Nn
            if(x(nn)>=MeanPrev)
                sign_x(nn)=1;
            else
                sign_x(nn)=-1;
            end
            if(dx(nn)>=0.0)
                sign_dx(nn)=1;
            else
                sign_dx(nn)=-1;
            end
            if(x_mean(nn)>=x_meanAvg)
                sign_x_mean(nn)=1;
            else
                sign_x_mean(nn)=-1;
            end
            if((x(nn)-MeanPrev)+dx(nn)>=0)
                sign_Sum(nn)=+1;
            else
                sign_Sum(nn)=-1;
            end
            sign_Sum(nn)=sign_x(nn);    
        end
        
        x_0(1:Nn)=x(1:Nn);
        for nn=1:Nn
            %This is the squared addition for t>1 time steps.
            x(nn)=MeanPrev+x_mean(nn)+sign_Sum(nn).*sqrt(abs(sign_x(nn)*  ...
                (x(nn)-MeanPrev).^2+sign_dx(nn).*dx(nn).^2));
        end
        
        
    else
        x(1:Nn)=x(1:Nn)+x_mean(1:Nn)+dx(1:Nn); %This is linear addition in 
                                %first time step when starting from delta origin. 
    end
end
%x(NnMidn)=(x(NnMidn+1)+x(NnMidn-1))/2;  %I just made this very quick and bad fix for the 
                                %unstable value in the mid. Mid three or
                                %four subdivisions have instability.
for nn=1:Nn
    y_x(nn) = ((1-gamma)*x(nn)).^(1/(1-gamma)); %Convert to original coordinates
end
for nn=2:Nn-1
    Dfx(nn) = (x(nn + 1) - x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfy_x(nn) = (y_x(nn + 1) - y_x(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end

fx(1:Nn)=0;
fy_x(1:Nn)=0;
for nn = 2:Nn-1
    fx(nn) = normpdf(Z(nn),0, 1)/abs(Dfx(nn)); %Bessel process density
    fy_x(nn) = normpdf(Z(nn),0, 1)/abs(Dfy_x(nn));%Origianl coordinates density
end
str=input('Analytic Values Calculated. Press a key to start Monte Carlo');
%Calculate the density with monte carlo below.
rng(85109634, 'twister')
paths=100000;
X(1:paths)=x0^(1-gamma)/(1-gamma); %Bessel process Monte carlo
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:Tt
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
    %HermiteP1(6,1:paths)=Random1(1:paths).^5-10*Random1(1:paths).^3+15*Random1(1:paths);
    %HermiteP1(7,1:paths)=Random1(1:paths).^6-15*Random1(1:paths).^4+45*Random1(1:paths).^2-15;
    XX0=X;
    YY0=YY;
    for k = 0 : OrderM
        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;
                    if(l4<=5)            
                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
                    mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                    ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
                    HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                    
                %Above is bessel process monte carlo equation.
                
                    YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                    mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                    YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                %Above is original coordinates monte carlo equation.
                    end
                end
            end
        end
    end
end

X(X<0)=0;
sum(X(:))/paths   %monte carlo average for bessel coordinates
sum(x(2:Nn-1).*ZProb(2:Nn-1))  % Ito-Taylor density average for bessel coordinates
sum(y_x(2:Nn-1).*ZProb(2:Nn-1)) %Original process average from coordinates 
                                %derived from bessel process 
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075;
MaxCutOff=30;
[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');
str=input('Look at Bessel graphs comparison between monte carlo density and analytic density');
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_x(2:Nn-1),fy_x(2:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('Look at Original variable graphs comparison between monte carlo and analytic density');
end

% 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. 

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

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

December 11th, 2018, 7:43 pm

Sorry Bad post.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

January 14th, 2019, 10:06 pm

We know that normal distribution and its transformations are known by Normal distribution's standard deviations. We can find CDF of any transformation of normal density by knowing the associated standard deviation.
We have that
                    Standard deviations=[$]Z(x,t)=\frac{x-\int_{t_0}^{t} \mu(x) dt}{\sqrt{\int_{t_0}^t {\sigma(x)}^2 dt}}[$]
When we are moving on a particular standard deviation line we have
                                                    [$]\Delta Z(x,t)=0[$]
For that to happen, we must have 
                                                  [$]\Delta Z(x,t)=\frac{\partial Z}{\partial x} \Delta x + \frac{\partial Z}{\partial t} \Delta t=0[$] 
 or in order to remain on the same standard deviation line, we can move in x direction in time [$]\Delta t[$] by magnitude given as
                                                 [$]\frac{\Delta x}{\Delta t}=-\frac{\frac{\partial Z}{\partial t}}{\frac{\partial Z}{\partial x}}[$]
where
            [$]\frac{\partial Z}{\partial x}=\frac{1-\frac{\partial \int_{t_0}^{t} \mu(x) dt}{\partial x}}{\sqrt{\int_{t_0}^t {\sigma(x)}^2 dt}}-(\frac{x-\int_{t_0}^{t} \mu(x) dt}{({\int_{t_0}^t {\sigma(x)}^2 dt})^{(-3/2)}}) \frac{\partial {\int_{t_0}^{t} {\sigma(x)}^2 dt}}{\partial x}[$]
and 
             [$]\frac{\partial Z}{\partial t}=\frac{-\frac{\partial \int_{t_0}^{t} \mu(x) dt}{\partial t}}{\sqrt{\int_{t_0}^t {\sigma(x)}^2 dt}}-(\frac{x-\int_{t_0}^{t} \mu(x) dt}{({\int_{t_0}^t {\sigma(x)}^2 dt})^{(-3/2)}})\frac{\partial {\int_{t_0}^{t} {\sigma(x)}^2 dt}}{\partial t}[$]
all the integrals of the type [$]\frac{\partial \int_{t_0}^{t} \mu(x) dt}{\partial x}[$], [$]\frac{\partial \int_{t_0}^{t} \mu(x) dt}{\partial t}[$], [$] \frac{\partial {\int_{t_0}^{t} {\sigma(x)}^2 dt}}{\partial x}[$] and [$] \frac{\partial {\int_{t_0}^{t} {\sigma(x)}^2 dt}}{\partial t}[$]
can be calculated from previous stochastic analysis.
Soon I hope to post a program based on the above analysis.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

January 21st, 2019, 7:18 pm

Here is derivation of an interesting density simulation scheme for simpler SDEs. In this scheme we continue to calculate drift and integrated vol on Standard deviation or constant CDF lines and continue to project the density into time. Ok here we go
On a certain constant CDF line, we have
[$]Z=\frac{x_2-\int_{0}^{t_2} \mu(x,t) dt - x_0}{\sqrt{\int_0^{t_2} \sigma(x,t)^2 dt}}[$]                 Eq(1)
suppose [$]x_1[$] is another point that lies on the same constant CDF line at time [$]t_1[$] so we must also have
[$]Z=\frac{x_1-\int_0^{t_1} \mu(x,t) dt -x_0}{\sqrt{\int_0^{t_1} \sigma(x,t)^2 dt}}[$]                     Eq(2)
From Eq(1), we do some algebra and we have that Eq(1) can also be written as
[$]Z=\frac{x_2-\int_0^{t_1} \mu(x,t) dt - \int_{t_1}^{t_2}\mu(x,t) dt + x_1 -x_1 -x_0}{\sqrt{\int_{0}^{t_2} \sigma(x,t)^2 dt}}[$]
We can further write the above expression as
[$]\frac{x_2-\int_{t_1}^{t_2} mu(x,t) dt}{\sqrt{\int_0^{t_2}\sigma(x,t)^2 dt}} - \frac{x_1-\int_0^{t_1} mu(x,t) dt -x_0}{\sqrt{\int_0^{t_2}\sigma(x,t)^2 dt}} =Z[$]
We divide the numerator and denominator of second term both by [$]\sqrt{\int_0^{t_1} \sigma(x,t)^2 dt}[$] and multiply both sides of the total equation by [$]\sqrt{\int_0^{t_2} \sigma(x,t)^2 dt}[$] to get the expression
[$]x_2-\int_{t_1}^{t_2} \mu(x,t) dt -x_1 + \sqrt{\int_0^{t_1} \sigma(x,t)^2 dt} \frac{x_1-\int_0^{t_1} \mu(x,t) dt -x_0}{\sqrt{\int_0^{t_1} \sigma(x,t)^2 dt}} =\sqrt{\int_0^{t_2}\sigma(x,t)^2 dt}  Z[$]
but from Eq(2) we also have
[$]x_2 -\int_{t_1}^{t_2} \mu(x,t) dt -x_1 + \sqrt{\int_0^{t_1} \sigma(x,t)^2 dt } Z = \sqrt{\int_0^{t_2} \sigma(x,t)^2 dt } Z[$]
So our density simulation scheme becomes along a particular CDF line Z(n), we have after a minor change of notation
[$]x(t_2)=x(t_1) + \int_{t_1}^{t_2} \mu(x,t) dt - \sqrt{\int_0^{t_1} \sigma(x,t)^2 dt } Z(n) + \sqrt{\int_0^{t_2} \sigma(x,t)^2 dt } Z(n)[$]

This is very handy since we can continue to travel along the nth CDF line and continue to calculate drift and integrated volatility on the same line and continue to evolve into time. In higher dimensions, this just needs integrated volatility net of correlations so very simple to evolve density in higher dimensions as well.
 
User avatar
Amin
Topic Author
Posts: 2085
Joined: July 14th, 2002, 3:00 am

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

January 22nd, 2019, 3:07 pm

What I have posted in the previous message seems to hold true in Bessel Coordinates for diffusions that do not have extreme non-linearities like mean-reverting SDEs and other non-linear SDEs. It does work for non-linear mean-reverting diffusions for smaller times but then starts to lose accuracy with increasing time. It does not seem to hold in the original coordinates for CEV type diffusions. I will soon be posting a program that will give a better idea for friends interested in this.
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...


GZIP: On