I am sharing this program and this could be a possible starting point for more advanced and precise algorithms. The program works for simple SDEs that do not take auto-correlation structure in integrated variance and you would have to make modifications in this program if you want to simulate mean-reverting or other SDE that take auto-correlation in time though this modification can be easily done. I am sure you would have to modify it for many other SDEs but I hope that it would be useful as a good starting point.

I have simulated the SDE given as

[$]dX(t)=\mu(t) X(t) dt + \sigma(t) X(t) dz(t)[$]

The program can take any term structure of volatilities and variable simulation steps. But it simulates one single unique SDE and you would have to make some changes in the program if you want to change the SDE half way along the simulation.

I have used the following simple algorithm to calculate the variance scaling of hermite polynomials.

I have equated the Ito-Taylor simulation method correlated integrated variance with a scaling factor for simulation standard deviation and equated it with uncorrelated true variance of the SDE evolution by altering the scaling factor that makes both sides of the equations equal i.e. this scaling factor equates the integrated variance of the Ito-Taylor density simulation method with the true integrated variance. This is the scaling that is applied to hermite polynomials. In case of auto-correlation in true variance evolution, you would have to account for the autocorrelation effect that I have not done. Here are the equations

if [$]\sigma_0^2[$] is the variance of existing density and [$]\sigma_1^2[$] is the variance of simulation step. Integrated variance of uncorrelated true SDE evolution is on the RHS of the equation below. And Ito-Taylor density simulation method integrated variance is on the LHS of equation below. To the LHS of the equation which is Ito-Taylor density integrated variance, I have applied a scaling denoted by [$]\rho[$] to the simulation step variance (its squre-root) [$]\sigma_1[$] which would equate both sides of the equation.

[$]{\sigma_0}^2 + 2 \sigma_0 (\rho \sigma_1) + {(\rho \sigma_1)}^2 = {\sigma_0}^2 + {\sigma_1}^2 [$]

so we get the equation for rho as

[$]{\rho}^2 + \frac{2 \sigma_0}{\sigma_1} \rho -1=0[$]

Though we could also calculate absolute variances in normal terms but we just need to find the ratio of standard deviations(meaning square roots of variances) which is done by transforming the densities into normal terms and taking the ratios of both transformed existing density and transformed simulation density subdivision interval.

[$]\rho = \frac{-b + \sqrt{b^2 +4}}{2}[$] where [$] b= \frac{2 \sigma_0}{\sigma_1}[$]

rho is the scaling of hermite polynomials that is needed to match the Ito-Taylor simulation method density with true density on each simulation step in time. Again here [$]\sigma_0[$] is square root of variance associated with existing density and [$]\sigma_1[$] is square root of variance associated with simulation step. And this ratio is calculated in "transformed normal density terms" and I would like to compare the normal density transformation as dimensionless-coordinates.

This is b in the above equation where we need to find the ratio of standard deviations.

Here is the matlab program. This is the first quick version and I hope to update it with a better program in a few days.

```
function [] = ItoTaylorNewWilmottQuadratic03ANew()
%This program simultes the density by making a comparison of
%existing variance of the density and variance to be simulated.
%The Ito-Taylro Density Simulaiton method adjusts the variance
%of hermite polynomials so that simulated variance matches
%with the true integrated variance of the SDE.
%The current program works in the case of SDEs with little
%auto-correlation strucutre across time. So In the case
%of mean-reverting stochastic volatility models or other
%similar models where auto-correlation exists in time,
%the method has to be appropriately modified.
%The first time step is directly simulated
%But for all later steps the variance of hermite
%polynomials is altered so that Ito-Taylor simulation
%method's integrated variance matches the true integrated variance
%The variance is matched by solution of the equation
% sigma1^2 + 2 * sigma1 *(r * sigma2)+ (r * sigma2)^2
%=sigma1^2 + 2 * sigma1 *sigma2+ sigma2^2
%The LHS of the equation is correlated variance
%equation with r which is the scaling factor.
%The RHS of the equation is uncorrelated variance
%evolution which will have to be altered in the
%case of SDEs that take auto-correlation.
% r is the scaling factor that is adjusted so
%that the variance of both sides of the equations are equal
%and this is the scalingfactor that adjusts the
%hermite polynomials's variance
Nn=161; % No of normal density subdivisions
dNn=.05; % Normal density width. would change with number of subdivisions
Tt=20; %No of total simulation time steps;
x0=1.0; % starting value of SDE
X(1:Nn)=x0; %Initialize the density subdivisions with initial value.
beta1=.55; % the first drift term power.
gamma=.850; % volatility power.
mu0=.26; % first drift coefficient.
dt(1)=.05; %I have kept the time steps the same but you can freely change them.
dt(2)=.05;
dt(3)=.05;
dt(4:Tt)=.065
sigma0(1)=.85;
sigma0(2)=.65;
sigma0(3)=.75;
sigma0(4)=.85;
sigma0(5)=.85;
sigma0(6)=.25;
sigma0(7)=.55;
sigma0(8)=.65;
sigma0(9)=.75;
sigma0(10:Tt)=.95;
T1=1.0;%Hermite polynomials have unit varaince for first step.
for tt=1:Tt
if(tt>1)
tt
NnMid=81; %This is the subdivision used for calculation
%of ratio of standard deviations.
for nn=NnMid-1:NnMid+1
Zn(nn)=((nn-1)*dNn - 4.0);
%This is a training step to calculate a ratio of
%one simulation step volatility of the density and
%volatility of exising density.
%Volatility in above line to be understood as square root of variance
%One step simulation density is calculated using a
%single starting point of x0
%but you can freely choose any other starting
%point since the variance of transformed normal
%density is independent of the starting coordinates
%though there could be some better starting points
% from the point of view of accuracy.
%This simulation step uses hermite polynomial with standard variance.
%This simulation step is re-started from time zero
%and not on top of existing density. This is just
%to calculate one step variance(or its square root)
%It works since we are calculating the density in
%normal or "dimensionless" coordinates.
X1(nn)=x0+ ...
mu0* x0.^beta1 * dt(tt) + ...
mu0^2 * beta1* x0.^(2*beta1-1) * dt(tt).^2/2.0 + ...
(sigma0(tt) * x0.^gamma *sqrt(dt(tt)) + ...
sigma0(tt) * mu0 * beta1 * x0.^(beta1 + gamma -1) *(1 - 1 / sqrt(3)) * dt(tt).^1.5 + ...
sigma0(tt) * mu0 * gamma * x0.^(beta1 + gamma -1) * 1 / sqrt(3) * dt(tt)^1.5) * Zn(nn) + ...
sigma0(tt).^2 *gamma* x0.^(2*gamma -1 ) .* dt(tt)/2 .* ( Zn(nn).^2 - 1 );
end
%Calculate the start and ends of mid density subdivisions
%y0h and y0l are start and end of mid of existing density
%before the simulation
y0h=X(NnMid)+(X(NnMid+1)-X(NnMid))/2.0;
y0l=X(NnMid)-(X(NnMid)-X(NnMid-1))/2.0;
%y1h and y1l are start and end of mid of density that has been
%simulated to calculate the variance after the simulation
y1h=X1(NnMid)+(X1(NnMid+1)-X1(NnMid))/2.0;
y1l=X1(NnMid)-(X1(NnMid)-X1(NnMid-1))/2.0;
% We calculate the standard deviation ratio after converting both
% densities to normal density after transformation.
% For a single unique SDE, we need to compare just one density subdivison
% and we just compare the ratio of varaince of middle subdivisions
% of existing and simulated training density
SigmaRatio=(2)*(y0h^(1-gamma)-y0l^(1-gamma))/(y1h^(1-gamma)-y1l^(1-gamma))
rho= (-SigmaRatio + sqrt(SigmaRatio^2 + 4))/2.0;
%This rho is scaling of hermite polynomials.
T1=rho^2*(1);
end
for nn=1:Nn
Z(nn)=((nn-1)*dNn-4.0)*sqrt(T1);
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
for nn=1:Nn
X(nn)=X(nn)+ ...
mu0* X(nn).^beta1 * dt(tt) + ...
mu0^2 * beta1* X(nn).^(2*beta1-1) * dt(tt)^2/2.0 + ...
(sigma0(tt) * X(nn)^gamma *sqrt(dt(tt)) + ...
sigma0(tt) * mu0 * beta1 * X(nn).^(beta1 + gamma -1) *(1 - 1 / sqrt(3)) * dt(tt)^1.5 + ...
sigma0(tt) * mu0 * gamma * X(nn).^(beta1 + gamma -1) * 1 / sqrt(3) * dt(tt)^1.5) * HermiteP(2,nn) + ...
sigma0(tt).^2 *gamma* X(nn).^(2*gamma -1 ) * dt(tt)/2 * HermiteP(3,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
%Monte carlo simulation follows for comparison of our simulated desnsity and
% monte carlo density.
paths=100000;
X2(1:paths)=x0;
Random1(1:paths)=0;
sigma00(1:Tt)=sigma0(1:Tt);
for tt=1:Tt
dt0=dt(tt);
Random1=randn(size(Random1));
X2(1:paths)=X2(1:paths)+ mu0* X2(1:paths).^beta1 *dt0 + ...
sigma00(tt) * X2(1:paths).^gamma .* Random1(1:paths) *sqrt(dt0) + ...
mu0^2 * beta1* X2(1:paths).^(2*beta1-1) * dt0^2/2.0 + ...
sigma00(tt)^2 * gamma * X2(1:paths).^(2*gamma-1) .* dt0/2.0 .* (Random1(1:paths).^2-1) + ...
mu0 * sigma00(tt) *gamma* X2(1:paths).^(beta1+gamma-1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ...
sigma00(tt) * mu0 *beta1* X2(1:paths).^(beta1+gamma-1) .* (1-sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths);
end
sum(X2(:))/paths %monte carlo average
BinSize=.0075;
MaxCutOff=20;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X2,paths,BinSize,MaxCutOff);
plot(X(3:Nn-2),fX(3:Nn-2),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
end
```