Here is a program that will calculate the density of an SDE using hermite polynomials. In the following program, I have used 24 hermite polynomials for expansion of driftless CEV SDE density at two years. I have kept the drift at zero and this is just an initial experimentation. I have not carefully checked everything and posting the program in the hope that it would be helpful for others. I will stay away from this project for a few days and possibly come back next week and make more changes. Here is the code.

```
In[83]:= Clear[x,t,x0,Y,P1,P2,l1,l2,Z,mu,sigma,mu0,sigma0,alpha, beta, gamma,a0,k,k1,m,L1,L2,YAns,H,Hz,X1,fX1,DfX1,Nn,n,HCoeff,ZCoeff,DZCoeff,zz];
In[84]:= P1=24;
In[85]:= P2=24;
In[86]:= L1=0;
In[87]:= L2=0;
In[88]:= Array[Y,{P1+L1+1,P1+L1+1,P2+L2+1},{0,0,0}];
In[89]:= For[k=0,k<=P1+L1,k++,(For[m=0,m<=P1+L1,m++,(For[n=0,n<=P2+L2,n++,(Y[k,m,n]=0);]);]);];
In[90]:= a0=1;
In[91]:= mu0=0;
In[92]:= sigma0=.75;
In[93]:= alpha=1;
In[94]:= gamma=.85;
In[95]:= beta=0;
In[96]:= Y[L1,L1,L2]=a0;
In[97]:= For[k=0,k<P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;
For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;Y[l1+1,l2,l3]=Y[l1+1,l2,l3]+SetPrecision[Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*mu0,50];
Y[l1,l2+1,l3]=Y[l1,l2+1,l3]+SetPrecision[Y[l1,l2,l3]*.5*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2)-1)*sigma0^2,50];
If[l3<L2+P2,(Y[l1,l2,l3+1]=Y[l1,l2,l3+1]+SetPrecision[Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*sigma0,50]),0];);];)];);];
In[98]:= Array[HCoeff,P2+L2+1,0];
In[99]:= For[k=0,k<P2+L2+1,k++,HCoeff[k]=0;]
In[100]:= x=1;
In[101]:= t=2;
In[102]:= For[k=0,k<=P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;
For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;
HCoeff[l3]=HCoeff[l3]+SetPrecision[Y[l1,l2,l3]*x^(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(Sqrt[t^(2*l1+2*l2-2*L1+l3)]*Sqrt[(2*l1+2*l2-2*L1)!])/(Sqrt[(2*l1+2*l2-2*L1+l3)!] *((l1+l2-L1)!))/Sqrt[l3!],50];)];)];);];
In[103]:= Array[Hz,{P2+L2+1,P2+L2+1},{0,0}];
In[104]:= For[k=0,k<P2+L2+1,k++,For[m=0,m<P2+L2+1,m++,Hz[k,m]=0;]];
In[105]:= Hz[0,0]=1.0;
In[106]:= For[k=0,k<P2+L2,k++,(For[m=0,m<=k,m++,(Hz[k+1,m+1]=SetPrecision[Hz[k,m],50];
If[m>=1,Hz[k+1,m-1]=SetPrecision[Hz[k+1,m-1]-m*Hz[k,m],50],0];)])];
In[107]:= Array[ZCoeff,P2+L2+1,0];
In[108]:= For[k=0,k<P2+L2+1,k++,ZCoeff[k]=0];
In[109]:= For[k=0,k<P2+L2+1,k++,(For[m=0,m<=k,m++,ZCoeff[m]=ZCoeff[m]+SetPrecision[Hz[k,m]*HCoeff[k],50]])];
In[110]:= Array[DZCoeff,P2+L2,0];
In[111]:= For[k=0,k<P2+L2,k++,DZCoeff[k]=0];
In[112]:= For[k=0,k<P2+L2,k++,DZCoeff[k]=SetPrecision[ZCoeff[k+1]*(k+1),50]];
In[113]:= Nn=81;
In[114]:= Array[X1,Nn,0];
In[115]:= Array[DfX1,Nn,0];
Array[fX1,Nn,0];
Array[zz,Nn,0];
In[118]:= For[n=0,n<Nn,n++,(X1[n]=0;DfX1[n]=0;fX1[n]=0;zz[n]=0;)];
In[119]:= For[n=0,n<Nn,n++,(zz[n]=n*.1-4;X1[n]=0;DfX1[n]=0;fX1[n]=0;
X1[n]=ZCoeff[0];DfX1[n]=DZCoeff[0];For[k=1,k<P2+L2+1,k++,(X1[n]=X1[n]+SetPrecision[ZCoeff[k]*zz[n]^k,50];If[k<P2+L2,DfX1[n]=DfX1[n]+SetPrecision[DZCoeff[k]*zz[n]^k,50],0];)];
fX1[n]=SetPrecision[PDF[NormalDistribution[0,1],zz[n]]/Abs[DfX1[n]],50];
)];
In[120]:= data1=Table[{X1[n],fX1[n]},{n,0,80}];
In[122]:= ListLinePlot[data1]
```