
amin


Total Posts: 279 
Joined: Aug 2005 


Dear friends, I was able to find very interesting applications of my research to nonlinear filtering that would make nonlinear filtering computationally very simple and efficient. The new method would work equally well for discrete time state space models. I will try to explain the idea with Bayesian filtering applied to a very simple set of univariate discrete time state space equations giving observation process equation and the State process equation. Before we work on a complete algorithm, there are a few principles that we need for that purpose.
1. Mapping an arbitrary density on a normal (or other) density. We can map any density on a normal density. Densities have probability mass normalized to one and this mapping requires that corresponding subdivisions on two densities have the same probability mass. Sometimes, we would like to map a density on a normal density so that subdivisions on the normal density are equally spaced.
2. Nonlinear transformation of a density We can easily take nonlinear transformation of a density and this requires scaling the density variable by the appropriate transformation but the density changes everywhere and new density for the transformed variable can be calculated from the pretransformation density by finding the change of variable derivative for the densities.
3. Adding locally nonlinear innovations to a density. Just like we did for SDEs, we can easily add local possibly nonlinear (in the SDE variable) normalbased innovations to a density and calculate the resulting density. For this we first map the density to a normal density and then add innovations linearly along the appropriate value of normal variable associated with the particular subdivision of the density.
4. Calculating Transition probabilities from one density to the other density when both densities have been mapped on a normal density. This is far simpler than people would otherwise think. And it would only require Green's function that takes into account appropriate variance of the innovations and most local nonlinearities can be easily neglected. For example, in the case of SDEs, we could map the two densities of SDE on Brownian motion densities(by equating the probability mass) and then transition probabilities between any two subdivisions would equal the transition probabilities between corresponding two subdivisions of the Brownian motion densities and this requires only appropriate calculation of time elapsed Δt. And any other nonlinearities in the SDE evolution can be safely neglected. Other nonlinearities show up in local scaling of the density but they can very easily be accounted for by a change of variable for densities. I will be explaining it more clearly with equations tomorrow.
5. Calculation of Filtered density once an Observation has been made. Suppose we have an updated state density and we possibly take a transformation of it (as dictated by observation/measurement equation) and then we add observation noise (or observation innovation) to it to calculate what I call the observation/measurement density. Calculation of filtered density from this observation/measurement density simply requires calculation of all the transition probabilities from the updated state density that would result in the observed point estimate(on the measurement/observation density). Again this filtered density is the transition density from all the points on the updated state density to the observed point on the observation/measurement density. This would not in general be a normalized density and we have to normalize it and this normalized density is the filtered density. And then we update this filtered density according to state update equation to calculate the updated state density. This would possibly require a transformation of the filtered density(2), mapping the transformed density on a normal density(1) and then adding local update innovations(3) required to calculate the state update density. And then we will again map the state update density on a normal density(1) and add observation noise/innovations(3) to calculate the measurement density. And then again find filtered density(5) again. Sorry that when I wrote about mapping a density of SDE on a normal density, I forgot to mention another important relevant fact and that is when densities of SDEs are generated by using ItoTaylor algorithm which gives us the SDE variable X,as a function of standard normal, Z, and in that case F(X(Z))=F(Z) where F means CDF of the density. So when SDE variable X is found as a function of standard normal variable Z, the CDF of X is exactly the same as the CDF of corresponding Z. And SDE variable X is a local scaling of the standard normal variable.
we can calculate the transition probabilities of SDEs by mapping the density of SDEs on relevant generating Brownian motion densities and then transition probability between various subdivisions on SDE density would be exactly the same as the transition probability between the corresponding subdivisions of two Brownian motion densities. Here is the link to the paper: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3119980 But I believe that normal density and densities of SDEs share a far stronger result which says that when we divide the evolving normal density or density of SDEs along CDF subdivisions or probability mass subdivisions, then the cumulative transition probabilities between different CDF subdivisions of both normal densities at different time remain the same. We are used to looking at transition probabilities in terms of absolute grids. Greater variance simply expands the CDF subdivisions and on fixed grids, it seems that density is expanding which is, of course, true with respect to fixed grids . But when considered in terms of CDF subdivisions (fixed probability mass subdivisions i.e subdivisions that expand so as to keep the probability mass the same in them), the local subdivisions also expand in a synchronized fashion so that transition probabilities between various expanding subdivisions really remain the same though on fixed grids one would only notice that normal density/SDE density is expanding while it remains constant when considered in terms of locally expanding fixed probability mass density subdivisions while transition probability dynamics remain constant in terms of expanding subdivisions.





amin


Total Posts: 279 
Joined: Aug 2005 


Here I made a post on linkedin for friends: https://www.linkedin.com/pulse/newalgorithmsnonlinearfilteringapplicationsmyresearchamin/
New Algorithms For Nonlinear Filtering: Applications of my Research On SDEs (Part I)
I was able to find very interesting applications of my research to nonlinear filtering that would make nonlinear filtering computationally very simple and efficient. The new method would work equally well for discrete time state space models. Before we work on a complete algorithm, there are a few principles that we need for that purpose. 1. Mapping an arbitrary density on a normal (or other) density. We can map any density on a normal density. Densities have probability mass normalized to one and this mapping requires that corresponding subdivisions on two densities have the same probability mass. Sometimes, we would like to map a density on a normal density so that subdivisions on the normal density are equally spaced. when densities of SDEs are generated by using ItoTaylor algorithm which gives us the SDE variable X,as a function of standard normal, Z, and in that case F(X(Z))=F(Z) where F means CDF of the density. So when SDE variable X is found as a function of standard normal variable Z, the CDF of X is exactly the same as the CDF of corresponding Z. And SDE variable X is a local scaling of the standard normal variable. 2. Nonlinear transformation of a density We can easily take nonlinear transformation of a density and this requires scaling the density variable by the appropriate transformation but the density changes everywhere and new density for the transformed variable can be calculated from the pretransformation density by finding the change of variable derivative for the densities. 3. Adding locally nonlinear innovations to a density. Just like we did for SDEs, we can easily add local possibly nonlinear (in the SDE variable) normalbased innovations to a density and calculate the resulting density. For this we first map the density to a normal density and then add innovations linearly along the appropriate value of normal variable associated with the particular subdivision of the density. 4. Calculating Transition probabilities from one density to the other density when both densities have been mapped on a normal density. This is far simpler than people would otherwise think. And it would only require Green's function that takes into account appropriate variance of the innovations and most local nonlinearities can be easily neglected. For example, in the case of SDEs, we could map the two densities of SDE on Brownian motion densities(by equating the probability mass) and then transition probabilities between any two subdivisions would equal the transition probabilities between corresponding two subdivisions of the Brownian motion densities and this requires only appropriate calculation of time elapsed Δt. And any other nonlinearities in the SDE evolution can be safely neglected. Other nonlinearities show up in local scaling of the density but they can very easily be accounted for by a change of variable for densities. As I mentioned earlier that we can calculate the transition probabilities of SDEs by mapping the density of SDEs on relevant generating Brownian motion densities and then transition probability between various subdivisions on SDE density would be exactly the same as the transition probability between the corresponding subdivisions of two Brownian motion densities. Here is the link to the paper: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3119980 But I believe that normal density and densities of SDEs based on locally nonlinear scaling of normal noise share a far stronger result which says that when we divide the evolving normal density or density of SDEs along CDF subdivisions or probability mass subdivisions, then the cumulative transition probabilities between different CDF subdivisions of both normal densities at different time remain the same. We are used to looking at transition probabilities in terms of absolute grids. Greater variance simply expands the CDF subdivisions and on fixed grids, it seems that density is expanding which is, of course, true with respect to fixed grids . But when considered in terms of CDF subdivisions (fixed probability mass subdivisions i.e subdivisions that expand so as to keep the probability mass the same in them), the local subdivisions also expand in a synchronized fashion so that transition probabilities between various expanding subdivisions really remain the same though on fixed grids one would only notice that normal density/SDE density is expanding while it remains constant when considered in terms of locally expanding fixed probability mass density subdivisions while transition probability dynamics remain constant in terms of expanding subdivisions. 5. Calculation of Filtered density once an Observation has been made. Suppose we have an updated state density and we possibly take a transformation of it (as dictated by observation/measurement equation) and then we add observation noise (or observation innovation) to it to calculate what I call the observation/measurement density. Calculation of filtered density from this observation/measurement density simply requires calculation of all the transition probabilities from the updated state density that would result in the observed point estimate(on the measurement/observation density). Again this filtered density is the transition density from all the points on the updated state density to the observed point on the observation/measurement density. This would not in general be a normalized density and we have to normalize it and this normalized density is the filtered density. And then we update this filtered density according to state update equation to calculate the updated state density. This would possibly require a transformation of the filtered density(2), mapping the transformed density on a normal density(1) and then adding local update innovations(3) required to calculate the state update density. And then we will again map the state update density on a normal density(1) and add observation noise/innovations(3) to calculate the measurement density. And then again find filtered density(5) again. I have written about this post in equal detail at post # and post of W****tt forum here: https://forum.W****tt.com/viewtopic.php?f=4&t=99702&start=720 And at nuclearphynance forum here:http://nuclearphynance.com/Show%20Post.aspx?PostIDKey=180518 



amin


Total Posts: 279 
Joined: Aug 2005 


The program has the right algorithm for linearlike SDEs but my algorithm using sqrt(1/T) *N for each step of normal density evolution probably does not seem to hold for very nonlinear 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 ItoTaylor 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 nonlinear 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)^(alpha1) * mu1 * X(t)^beta1 *dt +alpha*X(t)^(alpha1) * mu2 * X(t)^beta2*dt + .5* alpha *(alpha1)* X(t)^(alpha2) * sigma^2 * X(t)^(2*gamma) *dt + alpha * X(t)^(alpha1) * 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 *(alpha1)* 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 precompute everything but we have to calculate the part of above expression where power of x is evaluated. If we have four terms in Itohermite 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 Itohermite expansion to number of terms in polynomial expansion. For example Itohermite 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 precalculate 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 itohermite 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=(m11)*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=(m11)*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)^(alpha1) * mu1 * X(t)^beta1 *dt % +alpha* X(t)^(alpha1) * mu2 * X(t)^beta2 *dt % + .5* alpha *(alpha1)* X(t)^(alpha2) * sigma^2 * X(t)^(2*gamma) *dt % + alpha * X(t)^(alpha1) * 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 *(alpha1)* 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 precompute 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 Itohermite 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 Itohermite expansion to %number of terms in polynomial expansion. For example Itohermite 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 precalculate 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 itohermite %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)=((nn1)*dNn4.0)*1/sqrt(Tt); end
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1/sqrt(Tt)); ZProb(Nn)=1normcdf(.5*Z(Nn)+.5*Z(Nn1),0,1/sqrt(Tt)); for nn=2:Nn1 ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1/sqrt(Tt))normcdf(.5*Z(nn)+.5*Z(nn1),0,1/sqrt(Tt)); end
%ZProb contains probaility mass in each cell.
% Precompute 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)^21; HermiteP(4,nn)=Z(nn)^33*Z(nn); HermiteP(5,nn)=Z(nn)^46*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 %precalcualted 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^(k1); end if mu1 ~= 0 mu11(k)=mu1^(k1); end if mu2 ~= 0 mu22(k)=mu2^(k1); end if sigma0~=0 sigma22(k)=sigma0^(2*(k1)); 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;
%Precompute the time and power exponent values in small multidimensional 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^((l11) + (l21) + (l31) + .5* (l41)); Fp(l1,l2,l3,l4) = (alpha + (l11) * beta1 + (l21) * beta2 + (l31) * 2* gamma + (l41) * gamma ...  (l11)  (l21)  2* (l31)  (l41)); end end end end
%Below call the program that calculates the ItoTaylor 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)=((nn1)*dNn4.0); end
for nn=2:Nn1 Dfx(nn) = (x(nn + 1)  x(nn  1))/(B(nn + 1)  B(nn  1)); end fx(1:Nn)=0; for nn = 2:Nn1 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*beta11) * dt0^2/2.0 + ... mu2^2 * beta2 * X(1:paths).^(2*beta21) * dt0^2/2 + ... mu1*mu2* beta1* X(1:paths).^(beta1+beta21) *dt0^2/2.0 + ... mu1*mu2* beta2* X(1:paths).^(beta1+beta21) *dt0^2/2.0 + ... sigma0^2 * gamma * X(1:paths).^(2*gamma1) .* dt0/2.0 .* (Random1(1:paths).^21) + ... mu1 * sigma0 *gamma* X(1:paths).^(beta1+gamma1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ... mu2 * sigma0 *gamma* X(1:paths).^(beta2+gamma1) .* sqrt(1/3) .* dt0^1.5 .* Random1(1:paths) + ... sigma0 * mu1 *beta1* X(1:paths).^(beta1+gamma1) .* (1sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths) + ... sigma0 * mu2.*beta2* X(1:paths).^(beta2+gamma1) .* (1sqrt(1/3)) .* dt0^1.5 .* Random1(1:paths); end
sum(X(:))/paths %monte carlo average sum(x(:).*ZProb(:)) % ItoTaylor density average(good enough but slightly approximate)
BinSize=.00125; MaxCutOff=20; [XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti(X,paths,BinSize,MaxCutOff); plot(x(2:Nn1),fx(2:Nn1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
end 
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=(m11)*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=(m11)*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 dzintegral 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 Itohermite 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 Itohermite 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 polynomiallike 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=(m11)*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 Itohermite 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=CoeffDX11; ArrIndex0=(m11)*mDim+m2; ArrIndex=((m11)*mDim+(m21))*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 Itohermite expansion level starts and it is a triplynested 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 %nonzero 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=CoeffDX11; ArrIndex0=((m11)*mDim+(m21))*mDim+m3; %ArrIndex=(((m11)*mDim+(m21))*mDim+(m31))*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

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)=11/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)=(11/sqrt(3))*1/2*(11/sqrt(5)); IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(11/sqrt(5)); IntegralCoeff0(1,3,3,2)=1/2*(1sqrt(2)/2); IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5); IntegralCoeff0(1,3,2,3)=(11/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*(11/sqrt(7)); IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(11/sqrt(5))*1/3*(11/sqrt(7)); IntegralCoeff0(3,2,2,2)=(11/sqrt(3))*1/2*(11/sqrt(5))*1/3*(11/sqrt(7)); IntegralCoeff0(3,3,2,2)=1/2*(1sqrt(2)/2)*1/2*(1sqrt(2)/sqrt(6)); IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1sqrt(2)/sqrt(6)); IntegralCoeff0(3,2,3,2)=(11/sqrt(3))*1/sqrt(2)*1/2*1/2*(1sqrt(2)/sqrt(6)); IntegralCoeff0(3,3,3,2)=1/6*(1sqrt(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*(11/sqrt(5))*1/sqrt(2)*1/sqrt(6); IntegralCoeff0(3,2,2,3)=(11/sqrt(3))*1/2*(11/sqrt(5))*1/sqrt(2)*1/sqrt(6); IntegralCoeff0(3,3,2,3)=1/2*(1sqrt(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)=(11/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 dtintegral. %This dtintegral means that you are calculating the Itohermite expansion %of Integral X(t)^alpha dt with X(t) dynamics given by the SDE %here IntegralCoeffdz indicates the coefficients of a dzintegral. %This dzintegral means that you are calculating the Itohermite 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)*(1sqrt(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)*(1sqrt(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)*(1sqrt(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)*(1sqrt(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)*(1sqrt(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)*(1sqrt(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)*(1sqrt(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(XmaxXmax=real(X(p)); end end
IndexMax=floor((XmaxXmin)/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:(IndexMax1))*BinSize;
end









