This program is different from the previous program in the sense that cumulants being added due to diffusion term are explicitly calculated out of diffusion terms consisting of first and second hermite polynomials and their coefficients. Please let me explain. 
Let us say that higher order diffusion term in Bessel diffusion is of the form
[$]\sigma \sqrt{dt} \, Z+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, + \, c_3 \, w(t)^{\alpha_3} \, \, \Big] \, Z[$]
[$]+\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, + \, d_3 \, w(t)^{\beta_3} \, \, \Big] \, (Z^2 \, - 1)\, [$]
where c's and d's are coefficients that are found out of Ito-Hermite expansion as in high order monte carlo, while alphas and betas are powers on w that are found using Ito-Hermite expansion as in high order monte carlo.
We expand the w-terms in the form of a series in [$]Z_0[$] as
[$]\sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, + \, c_3 \, w(t)^{\alpha_3} \, \, \Big][$]
[$]=g_0 \, + \, g_1 Z_0 + \, g_2 {Z_0}^2 \, + \, g_3 {Z_0}^3 \, + \, g_4 {Z_0}^4[$]
here [$]Z_0[$] is standard normal associated with the existing diffusion prior to addition of volatility term and [$]Z_0[$] and [$]Z[$] are independent of each other.
while 
[$]\sigma \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, + \, d_3 \, w(t)^{\beta_3} \, \, \Big][$]
[$]=h_0 \, + \, h_1 Z_0 + \, h_2 {Z_0}^2 \, + \, h_3 {Z_0}^3 \, + \, h_4 {Z_0}^4[$]

So our initial volatility term is expanded as
[$]\sigma \sqrt{dt} \, Z+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, + \, c_3 \, w(t)^{\alpha_3} \, \, \Big] \, Z[$]
[$]+\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, + \, d_3 \, w(t)^{\beta_3} \, \, \Big] \, (Z^2 \, - 1)\, [$]
[$]=\sigma \sqrt{dt} \, Z \, + \Big[g_0 \, + \, g_1 Z_0 + \, g_2 {Z_0}^2 \, + \, g_3 {Z_0}^3 \, + \, g_4 {Z_0}^4 \, \Big] \, Z[$]
[$]+\Big[ \, h_0 \, + \, h_1 Z_0 + \, h_2 {Z_0}^2 \, + \, h_3 {Z_0}^3 \, + \, h_4 {Z_0}^4 \, \Big] \, (Z^2-1)[$]

Now we find value of cumulants associated with the diffusion term being added by taking appropriate powers of the above series expression on RHS and replacing the powers of two different and independent standard normals, [$]Z[$] and [$]Z_0[$] by their expected values. This is how I find the value of cumulants being added due to diffusion term. Mean of this volatility term is obviously zero.

Friends I recently posted a cumulant matching program that advances the density of SDE diffusions (represented in the form of a Z-series) by adding the cumulants of existing density with cumulants of innovations(diffusion term) and then finding another SDE density in Z-series so that its cumulants match the added cumulants of innovations and existing density. This technique and its derivatives can be applied in interesting ways at a lot of places in applied probability theory and I am sure many friends would like to use it extensively. Apart from a lot of simple applied uses, it can be used to solve for example the densities where the driver noise in not normal but for example is a gamma density variable.   
However backing out coefficients of Z-series by Newton method so that its cumulants match with desired cumulants might not always work easily for many diffusions. I was thinking of an alternative easier method that can find the Z-series representation with given target cumulants when a "nearby density" whose Z-series representation is known is input to the numerical method. By "nearby density" I mean a density that is not very different in cumulants/moments from target density and shares the support of the target density. Sharing the support is important since we want to find a new density by multiplying the original density with a polynomial and if original density has smaller tails/support, there is no way to find a good density beyond those tails by any possible polynomial multiplication where the support of original input density is zero. But we describe how to increase the support of the original input density by linearly scaling it when target density has far larger variance so the method could still be used. 
We want to multiply the original "nearby density" with a polynomial whose coefficients are as yet unknown and then numerically calculate first few moments by integrating over the density. Variables get integrated and by equating each moment, we are left with a linear equation for each moment in terms of as yet unknown coefficients. Zeroth moment is the probability distribution itself whose integral is equated to unity to find an extra linear equation in terms of unknown coefficients. We can solve for example six linear equations by gaussian elimination to retrieve the as yet unknown coefficients of the polynomial that multiplies the original "nearby density" and results in a density with deaired cumulants. 
Though I explain things again, please read this post 1023 on page 69 of this thread where I have used this method earlier: Breakthrough in the theory of stochastic differential equations and their simulation - Page 69 - wilmott.com.

Suppose we have an original nearby density f(x) on a grid. We can easily calculate all moments/cumulants by numerically integrating over the density. We want to find a target density g(x) that has certain cumulants usually different from the cumulants of nearby density. We multiply the original density with a polynomial with as yet unknown coefficients to find the target density.

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

We do not know the coefficient c's and want to find them so that cumulants or equivalently the moments we have backed out from cumulants are matched with cumulants of g(x), the target density. 
We write as

Since resulting density g(x) has to be a valid probability distribution, our first condition is that integral of density should  to one.

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)

if  [$]\mu_1[$] is the value of first raw moment backed out from target cumulants that we want to match, we must have 

[$]  \int  \,x \, g(x)  \, dx =\int  \, x\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_1  \, [$]   Equation(2)

similarly if if  [$]\mu_2[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^2 \, g(x)  \, dx =\int  \, x^2\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_2  \, [$]   Equation(3)

if  [$]\mu_3[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^3 \, g(x)  \, dx =\int  \, x^3\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_3  \, [$]   Equation(3)

similarly we can write equations of desired raw moments backed out from central moments to fourth or sixth or eighth order.
From equation (1), we have 

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)
[$] = c_0  \, \int  \, f(x)  \,  \, dx+ c_1  \int  \, f(x) \, x  \, dx + c_2  \int  \, f(x) \, x^2 \, dx + c_3 \int  \, f(x)  \, x^3 \, dx+ c_4  \int  \, f(x) \, x^4 \, dx  = \, 1  \, [$]

if we use the notation [$]E[x^n]=\int  \, f(x) \, x^n \, dx =  \zeta_n [$], we can write the moment equations as

[$]  c_0  \, + c_1 \, \zeta_1 + c_2 \, \zeta_2 + c_3 \, \zeta_3+ c_4 \, \zeta_4  = \, 1  \, [$]                            Equation(1)
[$]  c_0  \, \zeta_1 + c_1 \, \zeta_2 + c_2 \, \zeta_3 + c_3 \, \zeta_4+ c_4 \, \zeta_5  = \,\mu_1  \, \, [$]      Equation(2)
[$]  c_0  \, \zeta_2 + c_1 \, \zeta_3 + c_2 \, \zeta_4 + c_3 \, \zeta_5+ c_4 \, \zeta_6  = \,\mu_2  \, \, [$]      Equation(3)
[$]  c_0  \, \zeta_3 + c_1 \, \zeta_4 + c_2 \, \zeta_5 + c_3 \, \zeta_6+ c_4 \, \zeta_7  = \,\mu_3  \, \, [$]      Equation(4)

Since [$]E[x^n]=  \zeta_n [$]  are all known from numerical integration on original nearby density, we are left with a number of linear equations that we need to solve by Gaussian elimination to find out the coefficients c's of multiplying polynomial such that resulting density has cumulants equated with target cumulants. Once these equations are solved we find the resulting density as

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

There is one check that you will need. The density has to be positive everywhere. For well-stated moments with proper density support, this should not be a problem at all. But when bad cumulants are specified or support of the density does not match the support of target density, there could be a possibility that you would get a target density that is negative at some places as a result of solution of  above linear equations.

I will follow with one or two posts to complete description of the method. In new posts, I want to further explain how to use the above method when we have our original density f(x) as a Z_series and we want to find the Z-Series of target density g(x) such that cumulants of target density match the desired cumulants. I also want to mention simple scaling that can be used to increase the support of original density when its variance is smaller so its support matches the support of target density with larger variance. 

.
.
This was a good post. I would request friends to please disregard the post 1372 I made after this. But I thought more about the ideas and have something to share with friends.
In the copied post, we start with an initial density and we solve a set of linear equations to find a target density that matches target moments. I want to try to represent the target density variable as a transformation (given by power series) of the initial input density variable. Once we know this transformation, hopefully, we would even be able to evaluate expectations of functions of target density variable on initial density and  vice versa.

Suppose the initial density of variable X is given 
[$]p_X(x) =  f(x)  \, [$]

and we have the new target density 

[$]p_Y(x) =  f(x)   \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]=  f(x)   \, Q(x) [$] 
where
[$]Q(x)=[c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4][$]

We suppose that we can expand both f(x) and Q(x) as power series.
We want to find a transformation from X to Y,  [$]Y(X)[$] so that given any value(realization) of random variable X, the CDF of [$]Y(X)[$] on target density exactly matches the CDF of X on initial density. 
For that we write the CDF equations as

[$]  \int_{-\infty}^X  \, f(x) \, dx =  \, \int_{-\infty}^{Y(X)}  \, f(x)  \,Q(x) \, dx [$]  Equation(1)

X is natural variable for initial density so [$]\frac{dX}{dx}=1[$]

Differentiating our integral equations by x on both sides, we get

[$]    \, f(X) \, =  \,   \, f(Y(X))  \,Q(Y(X)) \, \frac{dY}{dX} [$]     Equation (A)

We had supposed that we could find power series of  f(X) and Q(X). Now we suppose that transformation relationship between X and Y is also given by a power series as
[$]Y(X)\,=\,a_0 \,+\, a_1  \,X \,+\, a_2\, X^2 \,+\, a_3\, X^3\, +\, a_4 \,X^4[$] Equation (B)

We notice that all terms in Equation(A) can be expanded as a power series. Equality between both sides of the equation would hold when power series expansions on both sides have equal coefficients on powers of X in which we are expanding the terms. We want to turn this equality of coefficients into equations(for each power of X) in terms of unknown coefficients [$]a_0\,, \,a_1 , \,a_2 , \,a_3 , \,a_4[$] . These equations would be algebraic expressions that can be solved possibly by root search. We can write rhs of equation (A) as an expansion as

[$]  \,   \, f(Y(X))  \,Q(Y(X)) \, \frac{dY}{dX} = \alpha_0(a_0,a_1,a_2,a_3,a_4) + \alpha_1(a_0,a_1,a_2,a_3,a_4) X + \alpha_2(a_0,a_1,a_2,a_3,a_4) X^2[$]
[$]+ \alpha_3(a_0,a_1,a_2,a_3,a_4) X^3+ \alpha_4(a_0,a_1,a_2,a_3,a_4) X^4 [$]  

while lhs of equation(A) can be expanded as

[$]  \,   \, f(X) = b_0 + \,b_1 X + \,b_2 X^2+ \,b_3 X^3+ \,b_4 X^4 [$]
which gives us five(or more) equations by matching coefficients of each power of X as  

[$]\alpha_0(a_0,a_1,a_2,a_3,a_4)=b_0 [$]
[$]\alpha_1(a_0,a_1,a_2,a_3,a_4)=b_1 [$]
[$]\alpha_2(a_0,a_1,a_2,a_3,a_4)=b_2 [$]
[$]\alpha_3(a_0,a_1,a_2,a_3,a_4)=b_3 [$]
[$]\alpha_4(a_0,a_1,a_2,a_3,a_4)=b_4 [$]
Each alpha above is an algebraic equation.
Once we solve the above equations by root search, we can find the transformation in the form of a power series that gives relationship between variable X of initial distribution with the variable Y of the target distribution as 

[$]Y(X)\,=\,a_0 \,+\, a_1  \,X \,+\, a_2\, X^2 \,+\, a_3\, X^3\, +\, a_4 \,X^4[$] Equation (B)

Once this relationship is known it could be possible to evaluate expectation (and its moments)  of variable Y(X) of the target density by integrating its series transformation relationship on the initial density as

[$]E[Y]=  \int_{-\infty}^{\infty} \,Y(X) \, f(X) \, dX =  \, \int_{-\infty}^{\infty} \, X\, f(X)  \,Q(X) \, dX[$]
or for its moments
[$]E[Y^n]=  \int_{-\infty}^{\infty} \,[Y(X)]^n \, f(X) \, dX =  \, \int_{-\infty}^{\infty} \, X^n\, f(X)  \,Q(X) \, dX[$]

So in the copied post we understood how to find a new target density whose cumulants match the target cumulants. And in this post, we have seen how to find a transformation between the two variables related to initial density and a different target density respectively.
I hope this post will be useful but please pardon any errors. 
In a few days I want to come back with a program that uses this method to simulate the densities of SDEs in time.

Friends, I continued to work yesterday after writing my previous post. And all through rest of the day yesterday, mind control agents continued to focus waves on my eyes(my eyes literally continued to tear/water for several hours). They also continued to focus waves on my feet, my penis and also my back. It was very difficult to work but I continued to bear the torture and did not mention it. But today before I started to work mind control agents started to force needles in my penis even though I was urinating in a bottle inside my room. They did it multiple times and for the second time when they forced needles into my penis, it was so painful that I was literally shouting with pain and then I decided to write this post. 
Mind control agents here in Pakistan are visibly disturbed that I want to continue my research. I want to tell friends that there is a very close gang of Pakistani American mind control agents and almost all of them are from brooklyn. These agents each make several million dollars every year and they are very very afraid that if mind control ends in Pakistan, huge money they swindle and divide among each other every year other would disappear. These agents remain in touch with Jews who are behind mind control in Pakistan and play upon fear and hatred of those jews by telling them about scenarios how fast this nation could develop if mind control ends in Pakistan. The malicious Jews in turn make American defense support them with all kind of money and embassy assistance. And this symbiotic relationship continues where mind control of hundreds of innocent Pakistanis like me continues and each mind control agent continues to mint millions of dollars every year. And of course mind control could not continue in Pakistan if it were not for crook generals in Pakistani military each of whom takes several tens of millions of dollars to continue mind control in Pakistan. 

I want to explain the calculation of variance of dt-integral with a toy example of four time periods. 
The gist of analytics below is that variances summation can only be applied on simulation intervals that are orthogonal to each other.
Let us first get oriented with variables and their notations.
We have an SDE for [$]X(t)[$] with four simulation periods on time points [$]t_0[$], [$]t_1[$], [$]t_2[$], [$]t_3[$], [$]t_4[$],
We suppose that SDE is such that its evolution within every time period is orthogonal to all other time periods. This holds for CEV noises and simple SDEs. For more difficult SDEs, you will have to take covariances into account carefully into this analysis.

We denote the evolution within each time period as
[$]\int_{t_0}^{t_1} dX(t) = X(t_1) \, - \, X(t_0) \,=\, X_{0,1} [$]
similarly for other simulation periods for example
[$]\int_{t_1}^{t_2} dX(t) = X(t_2) \, - \, X(t_1) \,=\, X_{1,2} [$]

please note that [$]X_{0,1} [$], [$]X_{1,2} [$], [$]X_{2,3} [$], [$]X_{3,4} [$], are all orthogonal to each other. 

Now we write X(t) at end of all simulation periods in terms of above orthogonal increments as

[$] X(t_1) \,= \, X_0  \,+  \, X_{0,1}[$]
[$] X(t_2) \,= \, X_0  \,+  \, X_{0,1}\,+  \, X_{1,2}[$]
[$] X(t_3) \,= \, X_0  \,+  \, X_{0,1}\,+  \, X_{1,2}\,+  \, X_{2,3}[$]
[$] X(t_4) \,= \, X_0  \,+  \, X_{0,1}\,+  \, X_{1,2}\,+  \, X_{2,3}\,+  \, X_{3,4}[$]

Now (this is approximate but I am sure that we can make it exact with proper integration over each time interval)
[$]\, \int_{t_0}^{t_4} \, X(t) \, dt \, \approx \, \sum_{t_n=t_0}^{t_4} \, X(t_n) \, \Delta t \,[$]
[$]=\, \Delta t \, X(t_1) \,+\,\Delta t \, X(t_2) \,+\,\Delta t \, X(t_3) \,+\,\Delta t \, X(t_4) \,[$] 
We want to convert above summation into orthogonal intervals so that we could apply summation of variances. Doing that
[$]=\, \Delta t \, [\, X_0  \,+  \, X_{0,1}] \,+\,\Delta t \, [\, X_0  \,+  \, X_{0,1}\,+  \, X_{1,2} ] \,+\, \Delta t \, [\, X_0  \,+  \, X_{0,1}\,+ \, X_{1,2}\,+  \, X_{2,3} ]\,[$]
[$]+\,\Delta t \, [\, X_0  \,+  \, X_{0,1}\,+  \, X_{1,2}\,+  \, X_{2,3}\,+  \, X_{3,4}] [$] 
[$]=4 \, \Delta t \, \, X_0  \,+4 \, \Delta t \, X_{0,1}\,+\, 3 \Delta t   \, X_{1,2} \,+2\, \Delta t  \, X_{2,3} \,+\,\Delta t  \, X_{3,4} [$] 
So we have
[$]\, \int_{t_0}^{t_4} \, X(t) \, dt \, \approx \, \sum_{t_n=t_0}^{t_4} \, X(t_n) \, \Delta t \,[$]
[$]=4 \, \Delta t \, \, X_0  \,+4 \, \Delta t \, X_{0,1}\,+\, 3 \Delta t   \, X_{1,2} \,+2\, \Delta t  \, X_{2,3} \,+\,\Delta t  \, X_{3,4} [$] 

Now matching variances (considering we are evolving from a delta point at [$]X_0[$] which has no variance
[$]\, Var[\int_{t_0}^{t_4} \, X(t) \, dt] \, \approx \, Var[\sum_{t_n=t_0}^{t_4} \, X(t_n) \, \Delta t \,][$]
[$]=16 \, {\Delta t}^2 \, Var[X_{0,1}]\,+\, 9 \, {\Delta t}^2   \, Var[X_{1,2}] \,+4\, {\Delta t}^2  \, Var[X_{2,3}] \,+\, {\Delta t}^2  \, Var[X_{3,4}] [$] 

I used the above logic to calculate the variances of dt-integrlas in the graphs for simple CEV SDEs that I attached a few posts earlier. I used simple dt summation in monte carlo with the same simulation period that I used for analytic calculations with hermite polynomials.
In next post I explain how I made the actual variance calculations with hermite polynomials.