difference with fokker-planck solution is that Fokker-planck solution automatically calculates any change in variance due to drift while other monte carlo-like(with random numbers replaced with stochastic integrals) would have to explicitly account for any contraction due to drift and would have to make extra calculations for that.

Friends would recall when I took the Fokker-planck solution and "backed out" standard normal variance increments between any two nodes at time tn and time tn+1. We were excited when we could accurately apply Ito change of variable formula using those Zt increments. But other SDEs were not working any more with those transition Zt increments and other dt-integrals were losing variance. Then we constructed an integrated Z(tn) from transition standard variance increments using transition probabilities. I was comparing variance of these reconstructed Z(tn) with variance of integrated values of standard normals used in monte carlo simulations but I was seriously puzzled that variance of reconstructed Z(tn) on fokker-planck grid using transition probabilities was far lesser than integrated values of normal constructed from integration of driving standard normals in monte carlo. Finally I blamed these discrepancies on inaccuracies of the Fokker-planck based transition probabilities grid and started an independent standard normal increments based grid.
Though there were indeed small inaccuracies on fokker-planck based transition probabilities grid but now it makes sense that underlying driving normal distribution had lesse variance because it contracted due to mean reverting drift and our backed out transition normals were quite close to accurate values and therefore Ito change of variable for that particular SDE was working very well. In fact due to variable local contraction due to drift, our retrieved "normal"(with integration of transition normal on fokker-planck based transition probabilities grid) was not even a perfectly straight line. However transition normal increments everywhere had unit variance(approximately different due to slight lack of accuracy). 
No wonder that we could not integrate any other SDE on the grid of our retrieved "normal"(with integration of transition normal on fokker-planck based transition probabilities grid).  And we also noticed that dt-integrals were fast losing variance on this Z-grid which again makes sense since the constructed normal grid had contraction due to drift. 
And now it all seems to make sense since Fokker-planck algorithm automatically calculates local adjustment in variance of underlying integrated normal grid due to mean-reverting and variance changing drift. 
I am sitting back for next 2-3 days brainstorming and thinking of the best way to go ahead on the transition probabilities project since a lot of things have dramatically changed due to better understanding of underlying phenomenon. 
Again our roadmap of research is the same. I want to complete the transition probabilities for Bessel SDEs part in a very good way and then move to caclulation of dt-integrals and dz-integrals(possibly on different sister grids) with transition probabilities. I think we would be able to solve general(not necessarily bessel form) SDEs which are counterparts of their bessel form on the same grids since underlying Z-grid remains the same.
And then we would move to monte carlo and transition probabilities grid generation of more general SDEs where each term would possibly include time functions(exponentials and powers of time and others) multiplying algebraic functions of SDE variable. So I hope that we would continue to do some exciting work in coming months and then move on to inference, parameter estimation, filtering, AI and other interesting things.

Though in our transition probabilities simulation of noises, we have used a normal distribution as underlying base distribution that increased with a unit variance but when mean-reverting or other general drift is added, the local variance changes differently due to non-linear drift and the underlying base distribution is no longer normal. This problem is no simpler than that of fokker-planck equation in bessel form. In fact fokker-planck solution of bessel SDEs is exactly increasing at uniform constant variance(though not unit variance) and contracting/expanding at the rate of drift' derivative with respect to underlying normal. In fact we have three distributions now. The first is underlying normal that we can take as standard. Second is a version of this normal that has taken account of all variance expansions and contractions and the third version is the solution of SDE in bessel form(or possibly also in original coordinates form.)
But again the problem of taking into account local expansion and contraction of underlying normal is almost exactly similar to the fokker solution of the SDE in bessel coordinates and no simpler.
When variance of underlying normal changes uniformly, it is very simple to add variance to it and every local subdivision expands linearly as a result and probability mass transfer across boundaries of various local subdivisions are perfectly symmetric and cancelling to preserve the mass in each local subdivision. But once there are non-linear local changes to initial existing variance, further variance cannot be added linearly by simple local squared addition since there are non-symmetric probability mass transfer across boundaries of various subdivisions due to initial existing non-uniform variance. I believe that we can calculate these non-symmetric mass transfers(with simple gaussian transition probability as our added variance is uniformly gaussian though existing variance is not and has been deformed due to drift) and then alter/calculate the boundaries/width of subdivisions that will keep the probability mass exactly the same in each subdivision as the initial probability mass with non-uniform variance. In a small simulation discretization step of the kind we have been taking in transition probabilities and FPE solution, the local mass transfer calculations will involve just a few neighboring subdivisions and should be fast enough. I want to try this as an alternative to route through solution of fokker-planck equation. I hope this will work very well but I could not be perfectly sure and will try it next few days.
I want also want to request friends to ask the mind control agencies to become better. They do not let me sleep well at night. They had also decreased forcing mind control chemicals in my mouth during my sleep for two or three days after I complained but then started doing the same thing just like old again. They also use gases at the table where I work and also in wash room adjacent to my room. And they still contaminate my food wherever they can as they did today when I had dinner. I was very excited this evening about sitting whole night to try to complete my experiments but when I had contaminated food at a street restaurant in Islamabad, a lot of my enthusiasm fizzled and I was not able to work well. Please force these mind control agencies to be better and stop trying to control me.

Since it would still be a few days before I complete the program perfectly, I want to describe for friends the new proposed method to solve for the simulation of SDEs in Bessel coordinates without resorting to FPE method.
Firstly there are three distributions. The SDE solution on nth node of the grid at time t given as [$]w_n(t)[$], second distribution is what I call Bessel distribution variable which keeps track of all local expansions and contractions of variance and is denoted as [$]B_n(t)[$] and the third is underlying standard normal distribution which  on nth node is denoted as [$]Z_n[$] and has no dependence on t.
Bessel distribution variable [$]B_n(t)[$] is only a displaced version of  [$]w_n(t)[$]. if we denote by [$]w0(t)[$] as the value of SDE solution on grid point where Z=0 (median of density) or [$]w0(t)= {w_n(t)|Z_n=0}[$], then Bessel distribution variable is given by [$]B_n(t)=w_n(t)-w0(t)[$]
This will also mean that at median of density [$]B_n(t)=Z_n=0[$] for the particular value of node n that lies at the median.
When the SDE starts from a delta origin, most SDEs(though of course, this may not hold for all SDEs) remain linear for a very small time like t=.025 so we take first few steps of the density with any linear algorithm(So many of them are possible).
After a few linear steps after the start from delta origin, we start our full-fledged "Bessel distribution variable" update algorithm.
We have B_n(t) given and we want to construct B_n(t+1).
We start to construct [$]B_n(t+1)[$] from the node at the median. At median B_n(t+1) always remains zero. To construct Bessel variable to right of median, We have starting(left) boundary of this first subdivision to the right of median exactly at median where B(t+1)=0. We approximate the right boundary of this first subdivision to right of median at [$]\sqrt{({\frac{\partial B}{\partial Z}}^2 + \sigma^2 dt)}[$]. This right boundary is a trial boundary and would work perfectly well if Bessel variable distribution were a normal distribution. Since both boundaries have been set, we can now use Gaussian formula to calculate actual probability mass transfers from all subdivisions at time t, to this subdivision at time t+1. We want the grid to expand so that integrated mass in every subdivision always remains constant at [$]P_n[$]. However the actual integrated probability mass would be very slightly different from the target constant mass. We calculate this difference(excess probability mass) and call it [$]\Delta P_n[$] and use the Taylor formula as [$]\Delta P_n= \frac{\partial P}{\partial B} \Delta B +  \frac{\partial^2 P}{\partial B^2} {(\Delta B)}^2[$] to calculate how much to move the right boundary by [$]\Delta B[$] so that new [$]\Delta P_n=0[$]. We will have to solve it with special form of quadratic formula that is robust for small coefficients on squared term of quadratic equation.
I propose to approximate [$]\frac{\partial P}{\partial B}=\frac{-Z}{\sqrt{({\frac{\partial B}{\partial Z}}^2 + \sigma^2 dt)}} P(Z)[$]
and similarly for second derivative. Please note that all derivatives in above equations are calculated differently for each node and I have not shown dependence of these derivatives on node variable n.
Once we have fixed the right boundary of first subdivision to the right of median, we will take it as left boundary for second subdivision and set a trial right boundary for second subdivision to right of median and repeat the calculation of excess probability mass and then calculate how much to move our right trial boundary so that excess probability mass goes to zero.
I hope to present the worked out program in a few days. There are a lot of interesting things yet to be seen if we could take a large step with this algorithm and we will find out once the algorithm is ready.

Friends, I am still working on the problem outlined in previous post and a lot of things have changed from the plan of attack I described.
I am very sure many smarter friends would have already realized that once we have found the "Bessel distribution variable" by matching its origin with origin of normal distribution, we can very easily add diffusion to it with a PDE and then find a restored version of SDE variable in bessel coordinates and then add drift. PDEs are extremely fast and very accurate for pure diffusions and drift can be later added independently. Decoupling the addition of variance to SDE simulation and and addition of drift to SDE(in bessel coordinates) simulation using a PDE could be a very powerful tool(as handling strong convections is difficult for numerical solution of PDEs while adding gaussian diffusion is probably almost trivial.)
I hope to complete the program with transition of probability mass within another two days. Meanwhile please read here
https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=864253#p864253
 my post with request to new president Joe Biden about ending mind control. In my post, I have also request CNN, New York Times, and Washington Post to highlight this issue of grave human rights abuses. I have also shared my (bad) opinion about previous president obama.
viewtopic.php?f=15&t=94796&start=1200#p864275
viewtopic.php?f=15&t=94796&start=1200#p864276

I want to tell friends that my research in mathematics and my freedom from mind control are two most important things for me at this point in life. I was expecting that my freedom from mind control would be automatic after arrival of a new democratic president in office. Though it may not be very far in future but I am afraid that it may be a protracted battle for quite sometime. I am writing this post to tell friends that keeping focus on my research is extremely important for me. Though I would be spending slightly more time in my attempts at my and other victims' freedom from mind control and request friends to help as well, I will make sure that my major focus on research in mathematics only increases and not decreases. And I want to thank the readers for taking the time to read my posts and I hope to continue to post quality research here.

Here I have again requested President Biden to end mind control. And commented how to keep America in a leader on today's world. 
https://forum.wilmott.com/viewtopic.php?f=15&p=864314#p864314

Friends, my work on the program(for simulation of SDEs in transition probabilities related framework) slowed down but I have almost completed coding and my logic is sound. I will be debugging and running the complete program in another day. I hope the program has desired accuracy that we want to see. 
Once this program is complete, I will come back with a second program with numerical PDE type solution as I mentioned two posts ago. I also have a very good idea now why our simulation of (dz and dt) integrals of SDE was not working well on the same transition probabilities grid. We will have to make different independent sister-grid for each dz-integral or dt-integral. However each grid will be related to main SDE grid based on corresponding Z-lines(or Z-grid). If main SDE variable is w, We will copy the value of  f(w(Z_n))  from nth node of the parent SDE grid and use it for simulation on nth node of sister grid so both parent and sister grids would be related with respect to Z-lines(or density CDF). This way we would be able to simulate the dt-integrals and dz-integrals of functions of SDE on a single dimensional transition probabilities grid(or possibly have a numerical PDE type method). 
Friends would recall that we could simulate the direct functions of SDE when I tried to impose transition probabilities on FPE solution however we could not simulate any integrals of the SDE on the same grid. Both grids for SDE and its (dz and dt) integrals have very different variances with their own contraction and expansion effect of their drifts. The (dz and dt) integrals grid will only borrow the value of SDE variable along Z-lines from the original SDE grid and use it on corresponding Z-lines on sister grid and then have independent evolution of drift and variance on sister grid according to existing variance of the (dz- and dt integral) sister grid.
I am sorry that my work slowed down a bit but I hope to speed up in next few days if things do not get worse.
I also want to tell friends that crooks in US defense have hacked my website since it was a wonderful website and gave a great image and they did not want many other visitors to see my website as it would give them a good image of me. God, When are these immature crooks in US Army ever be able to grow up? 

Dear friends, I have almost completed my program for simulation of SDEs that I mentioned in previous few posts and its main algorithm works very well. I will be making a few small changes and post the program in another day or two. 
I want to take this opportunity to thank those American friends who insisted on mind control agencies to not force injections on me and therefore I was able to complete this program relatively quickly and reasonably. Though mind control is continuing and mind control agencies are still trying to give me drugged food, not forcing injections on me is a huge step for me.
Anyway, I hope to post the program in a day or two and I hope friends would like it. I also want to thank friends for helping me in making it possible. 
I will post a generic full-fledged working algorithm first. Since it is a new algorithm, I think a lot of optimizations and changes are possible and therefore I hope to come up with a more optimized version in a few more days. Though I am sure that original prototype will also work very well.

Friends, I have worked out the algorithm to update the density for the simulation of SDEs in Bessel coordinates without resorting to FPE method. I am going to give the details of the new modified algorithm below. I have copied some of the material from a previous post to give the right context and avoid having to write the same introduction again.

The SDE solution on nth node of the grid at time t given as [$]w_n(t)[$], second distribution is what I call Bessel distribution variable which keeps track of all local expansions and contractions of variance and is denoted as [$]B_n(t)[$] and the third is underlying standard normal distribution which  on nth node is denoted as [$]Z_n[$] and has no dependence on t.
Bessel distribution variable [$]B_n(t)[$] is only a displaced version of  [$]w_n(t)[$]. if we denote by [$]w0(t)[$] as the value of SDE solution on grid point where Z=0 (median of density) or [$]w0(t)= {w_n(t)|Z_n=0}[$], then Bessel distribution variable is given by [$]B_n(t)=w_n(t)-w0(t)[$]
This will also mean that at median of density [$]B_n(t)=Z_n=0[$] for the particular value of node n that lies at the median.
When the SDE starts from a delta origin, most SDEs(though of course, this may not hold for all SDEs) remain linear for a very small time like t=.025 so we take first few steps of the density with any linear algorithm(So many of them are possible).
After a few linear steps after the start from delta origin, we start our full-fledged "Bessel distribution variable" update algorithm.
We have B_n(t) given and we want to construct B_n(t+1).
We start to construct [$]B_n(t+1)[$] from the node at the median. At median B_n(t+1) always remains zero. To construct Bessel variable to right of median, We have starting(left) boundary of this first subdivision to the right of median exactly at median where B(t+1)=0. We approximate the right boundary of this first subdivision to right of median at [$]\sqrt{({\frac{\partial B}{\partial Z}}^2 + \sigma^2 dt)}[$]. This right boundary is a trial boundary and would work perfectly well if Bessel variable distribution were a normal distribution. Since both boundaries have been set, we can now use Gaussian formula to calculate actual probability mass transfers from all subdivisions at time t, to this subdivision at time t+1. We want the grid to expand/contract so that integrated mass in every subdivision always remains constant at [$]P_n[$]. However the actual integrated probability mass would be very slightly different from the target constant mass. We calculate this difference(excess probability mass) and call it [$]\Delta P_n[$].
What follows is different from what I described in previous post.
We take the right trial boundary and calculate the transition probability from each subdivision at previous time at this trial boundary. From a cell m at time t, to the right boundary of cell n, this transition probability is denoted as p(t,m,t+1,n). If the center of mth cell at t is given by [$]B_m[$] and the right trial boundary of nth  cell at time t+1 is given as[$]\bar{B_n}[$], the transition increment between the two would be [$]\bar{B_n}-B_m[$]. In case of bessel SDE setting, this transition increment would be locally normal with volatility [$]\sigma_m[$]  if the time step is not too large. So  [$] p(t,m,t+1,n)=p_m=N(\bar{B_n}-B_m,0,\sigma_m)[$]
Since probability mass in original subdivision is not one and is equal to [$]P_m[$], we will have to multiply the transition probability by [$]P_m[$]. In what follows, I will write [$] p(t,m,t+1,n)=p_m=P_m N(\bar{B_n}-B_m,0,\sigma_m)[$]
First derivative of this transition probability at the right trial boundary is given as[$]\frac{dp_m}{d\bar{B}_n}=P_m \frac{-1}{\sigma_m} \frac{(\bar{B_n}-B_m)}{\sigma_m} N(\bar{B_n}-B_m,0,\sigma_m)[$]
Second derivative of this transition probability at the right trial boundary is given as[$]\frac{d^2p_m}{d{\bar{B}_n}^2}=P_m \frac{(-1)^2}{{(\sigma_m)}^2} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^2-1) N(\bar{B_n}-B_m,0,\sigma_m)[$]
Third derivative of this transition probability at the right trial boundary is given as[$]\frac{d^3p_m}{d{\bar{B}_n}^3}=P_m \frac{(-1)^3}{{(\sigma_m)}^3} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^3-3(\frac{(\bar{B_n}-B_m)}{\sigma_m})) N(\bar{B_n}-B_m,0,\sigma_m)[$]
The above are transition probability and its derivatives and if we want to convert them to CDF/probability mass, we would have to integrate them. 
if we move the right trial boundary by a very small amount [$]\Delta B[$], the change in probability mass coming into nth cell at time t+1 from mth cell at time t, would be approximated after integration as

[$] p_m \Delta B+ \frac{dp_m}{d\bar{B}_n} { \Delta B}^2/2 + \frac{d^2p_m}{d{\bar{B}_n}^2} { \Delta B}^3/6+ \frac{d^3p_m}{d{\bar{B}_n}^3} { \Delta B}^4/24=\Delta P_n[$]
The above was equation for change in probability mass coming from mth cell at time t as we moved the trial boundary by a very slight bit. All the values of transition probability and their derivatives are calculated exactly at the trial boundary.
Now we want one single equation for total change in probability mass coming from all of the M cells at time t into nth cell at time t+1 as we move the right trial boundary by a little bit. This is very simple and a mere addition over all M cells and can be written as
[$]\Delta B \sum_{1}^{M} p_m +{ \Delta B}^2/2  \sum_{1}^{M} \frac{dp_m}{d\bar{B}_n}  +{ \Delta B}^3/6 \sum_{1}^{M} \frac{d^2p_m}{d{\bar{B}_n}^2} +{ \Delta B}^4/24 \sum_{1}^{M} \frac{d^3p_m}{d{\bar{B}_n}^3} =\Delta P_n[$]
I have omitted the summation sign on the excess probability mass in RHS. Obviously This has to be taken as summation of excess prob mass after taking into account diffusion from all M cells on grid at time t.
So once we know the total excess probability mass with the trial boundary given as [$]\Delta P_n[$], we can easily calculate how much do we have to move bcak/forward [$]\Delta B[$] the trial boundary so that the excess mass is exactly reversed. We would have to solve the above quartic equation as
[$]a \Delta B  +b { \Delta B}^2   +c { \Delta B}^3 +d { \Delta B}^4 =-\Delta P_n[$]
[$]a= \sum_{1}^{M} p_m[$]
[$]b=\sum_{1}^{M} \frac{dp_m}{d\bar{B}_n}/2[$]
[$]c=\sum_{1}^{M} \frac{d^2p_m}{d{\bar{B}_n}^2}/6[$]
[$]d=\sum_{1}^{M} \frac{d^3p_m}{d{\bar{B}_n}^3}/24[$]
Sometimes special care has to be taken in choice of signs as we move right or left away from the median and if excess probability mass is positive or negative. I calculated excess probability mass by calculating transition CDF from each of the M cells at time t to both left boundary(that was fixed) and to the right boundary of nth cell at time t+1 and subtracting the CDF on both boundaries.

I was earlier doing the above analytics with just two terms with solving a quadratic equation. Though results were very good and I was able to closely calculate the required change in the boundary of each cell within the Zn= -4 SD to Zn=+4 SD but outside to +5SD, there were some problems since two term quadratic approximation was not perfectly good and small inaccuracies would add up.
In my algorithm, I have also taken special measures for time t grid cells overlapping with t+1 grid cells. If a grid cell at time t is such that a part of it is to the right of the trial boundary(of t+1 cell) and part of it is to the left of the trial boundary, I divide the cell into two parts so that one part is completely to the left of the boundary and the other is completely to the right of the boundary so as to be able to calculate valid transition probabilities. My grid is also large (200 cells) so it seems to work well. 
I am now solving a four term quartic equation and I solved it in mathematica with general coefficients and later extremely simplified it with repeated algebraic substitutions and it is quite cheap to calculate all four roots of quartic equation. But I still have to do some work with robust selection of right root from four roots.  
Quartic equation is particularly interesting since it locally approximates the true change in probability mass extremely well and Fokker-planck solution also included a third derivative. I hope we would be able to take a reasonably large step with this algorithm.
Once we have fixed the right boundary of first subdivision to the right of median, we will take it as left boundary for second subdivision and set a trial right boundary for second subdivision to right of median and repeat the calculation of excess probability mass and then calculate how much to move our right trial boundary so that excess probability mass goes to zero.
I hope to post the new program in another few days.
I also want to take this opportunity to tell friends that mind control agents are doing their full effort to give me drugged food. I have been able to get good food and water for past few days but it might become difficult to outsmart them every day and that might effect my work. Please force them to take it easy and let me continue to devote time actively to my research.
I have also told friends before that my research  and programs are given to some people in academia who are friendly to mind control agencies and if you know somebody who claims to do exactly the same thing, please weigh the chances of their doing it on their own. I really think that one reason behind frenetic attempts to drug my food in past days was to not let me work after having drugged food and ask other people to complete the work in the meantime and present it as their own.

In previous post, I forgot to scale the transition probability from mth subdivision at time t to nth subdivision at time t+1 with probability mass in original mth subdivision given as [$]P_m[$]. This would give slightly  wrong idea about the algorithm. I altered one paragraph in previous post and I am copying the altered paragraph here.

We take the right trial boundary and calculate the transition probability from each subdivision at previous time at this trial boundary. From a cell m at time t, to the right boundary of cell n, this transition probability is denoted as p(t,m,t+1,n). If the center of mth cell at t is given by [$]B_m[$] and the right trial boundary of nth  cell at time t+1 is given as[$]\bar{B_n}[$], the transition increment between the two would be [$]\bar{B_n}-B_m[$]. In case of bessel SDE setting, this transition increment would be locally normal with volatility [$]\sigma_m[$]  if the time step is not too large. So  [$] p(t,m,t+1,n)=p_m=N(\bar{B_n}-B_m,0,\sigma_m)[$]
Since probability mass in original subdivision is not one and is equal to [$]P_m[$], we will have to scale the transition probability by [$]P_m[$]. In what follows, I will write [$] p(t,m,t+1,n)=p_m=P_m N(\bar{B_n}-B_m,0,\sigma_m)[$]
First derivative of this transition probability at the right trial boundary is given as[$]\frac{dp_m}{d\bar{B}_n}=P_m \frac{-1}{\sigma_m} \frac{(\bar{B_n}-B_m)}{\sigma_m} N(\bar{B_n}-B_m,0,\sigma_m)[$]
Second derivative of this transition probability at the right trial boundary is given as[$]\frac{d^2p_m}{d{\bar{B}_n}^2}=P_m \frac{(-1)^2}{{(\sigma_m)}^2} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^2-1) N(\bar{B_n}-B_m,0,\sigma_m)[$]
Third derivative of this transition probability at the right trial boundary is given as[$]\frac{d^3p_m}{d{\bar{B}_n}^3}=P_m \frac{(-1)^3}{{(\sigma_m)}^3} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^3-3(\frac{(\bar{B_n}-B_m)}{\sigma_m})) N(\bar{B_n}-B_m,0,\sigma_m)[$]

 

In previous post, I forgot to scale the transition probability from mth subdivision at time t to nth subdivision at time t+1 with probability mass in original mth subdivision given as [$]P_m[$]. This would give slightly  wrong idea about the algorithm. I altered one paragraph in previous post and I am copying the altered paragraph here.

We take the right trial boundary and calculate the transition probability from each subdivision at previous time at this trial boundary. From a cell m at time t, to the right boundary of cell n, this transition probability is denoted as p(t,m,t+1,n). If the center of mth cell at t is given by [$]B_m[$] and the right trial boundary of nth  cell at time t+1 is given as[$]\bar{B_n}[$], the transition increment between the two would be [$]\bar{B_n}-B_m[$]. In case of bessel SDE setting, this transition increment would be locally normal with volatility [$]\sigma_m[$]  if the time step is not too large. So  [$] p(t,m,t+1,n)=p_m=N(\bar{B_n}-B_m,0,\sigma_m)[$]
Since probability mass in original subdivision is not one and is equal to [$]P_m[$], we will have to scale the transition probability by [$]P_m[$]. In what follows, I will write [$] p(t,m,t+1,n)=p_m=P_m N(\bar{B_n}-B_m,0,\sigma_m)[$]
First derivative of this transition probability at the right trial boundary is given as[$]\frac{dp_m}{d\bar{B}_n}=P_m \frac{-1}{\sigma_m} \frac{(\bar{B_n}-B_m)}{\sigma_m} N(\bar{B_n}-B_m,0,\sigma_m)[$]
Second derivative of this transition probability at the right trial boundary is given as[$]\frac{d^2p_m}{d{\bar{B}_n}^2}=P_m \frac{(-1)^2}{{(\sigma_m)}^2} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^2-1) N(\bar{B_n}-B_m,0,\sigma_m)[$]
Third derivative of this transition probability at the right trial boundary is given as[$]\frac{d^3p_m}{d{\bar{B}_n}^3}=P_m \frac{(-1)^3}{{(\sigma_m)}^3} ({(\frac{(\bar{B_n}-B_m)}{\sigma_m})}^3-3(\frac{(\bar{B_n}-B_m)}{\sigma_m})) N(\bar{B_n}-B_m,0,\sigma_m)[$]

 

Many friends would have noticed that [$]P_m[$] in above equations is a source of error. I am assuming as if all the mass in mth subdivision at time t is concentrated as delta mass at its centre. Obviously this is wrong. Distribution of mass within the subdivision is not linear or symmetric and transition probability is also not linear at all. Because of both these reasons there will be some error in previous analysis. But I have done analytics to calculate the correction(due to both of the mentioned factors) that when added will take our previous calculations to exact values. I will be sharing it in a few hours or by tomorrow at the latest.

In this post, I want to write how to calculate exact Gaussian transition CDF, Gaussian transition probability and derivatives of Gaussian transition probability from mth Bessel grid cell  [$]B_m[$] at time t, to any point on the new grid at time t. Please note that we are not calculating Gaussian CDF and transition probability merely between two points. Rather we have an initial Bessel grid where we take an arbitrary mth cell and we want to calculate exact mass transfer CDF, transition probability and its derivatives from the entire mth cell at time t to nth cell at time t+1. As I indicated in my previous post there is an obvious difference between mass transfer(its CDF and pdf) between two points and mass transfer from a cell of certain width to another point and both of these are different and replacing mass transfer between a cell and a point by mass transfer between two points can lead to increasing accuracies in any transition probabilities algorithm.

Let us fix the notation first which is slightly different from what I have been using earlier. We have a Bessel grid and we indicate an arbitrary cell by [$]B_m[$] and the center of this cell is given by [$]B_{m0}[$]. This center [$]B_{m0}[$] will correspond to point [$]Z_m[$] on Z-grid and in general both boundaries of the cell(that will correspond to [$]Z_m+\Delta Z/2[$] and [$]Z_m-\Delta Z/2[$]) will not be equidistant from this center of the cell [$]B_{m0}[$]. We will represent the left boundary of cell [$]B_m[$] corresponding to [$]Z_m-\Delta Z/2[$] as  [$]B_{m1}[$] and right boundary corresponding to  [$]Z_m+\Delta Z/2[$]  as [$]B_{m2}[$]
The transition probability between center of mth cell [$]B_{m0}[$] to arbitary point [$]B_n[$] on next grid is denoted as [$] p_{m0,n}=N(B_n-B_{m0},0,\sigma_m)[$].
Its CDF is defined by uppercase of PDF as [$] P_{m0,n}=\int_{-\infty}^{(B_n-B_{m0})} N(B_n-B_{m0},0,\sigma_m) dN [$]
I want to emphasize that above PDF, CDF and their derivatives below are between two points(The centre of [$]B_m[$] given by [$]B_{m0}[$] and point [$]B_n[$]
First derivative of this transition probability [$]p_{m0,n}[$] at point [$]B_{m0}[$] is given as[$]\frac{dp_{m0,n}}{dB}=(-1) \frac{-1}{\sigma_m} \frac{(B_n-B_{m0})}{\sigma_m} N(B_n-B_{m0},0,\sigma_m)[$]
Second derivative of this transition probability [$]p_{m0,n}[$] at point [$]B_{m0}[$] is given as[$]\frac{d^2p_{m0,n}}{d{B}^2}=\frac{(-1)^2}{{(\sigma_m)}^2} ({(\frac{(B_n-B_{m0})}{\sigma_m})}^2-1) N(B_n-B_{m0},0,\sigma_m)[$]
Third derivative of this transition probability [$]p_{m0,n}[$] at point [$]B_{m0}[$] is given as[$]\frac{d^3p_{m0,n}}{d{B}^3}=(-1)*\frac{(-1)^3}{{(\sigma_m)}^3} ({(\frac{(B_n-B_{m0})}{\sigma_m})}^3-3(\frac{(B_n-B_{m0})}{\sigma_m})) N(B_n-B_{m0},0,\sigma_m)[$]

We define the probability density function at center of mth cell in Bessel variable [$]B_m[$] as lowercase [$]p_{m0}=p(Z_m) \frac{dZ}{dB_m}[$]
while uppercase will denote its CDF given as [$]P_{m0}=P(Z_m)[$]
The probability mass in mth cell which is the difference of CDF at boundaries is denoted as [$]\Delta P_{m}[$] 
The derivatives of  probability density function at center of mth cell in Bessel variable [$]B_m[$] are given as[$]\frac{dp_{m0}}{dB}=\frac{dp(Z_m)}{dZ} {(\frac{dZ}{dB_m})}^2+p(Z_m) \frac{d^2Z}{d{B_m}^2} [$]
We can similarly calculate higher derivatives [$]\frac{d^2p_{m0}}{d{B}^2}(B_{m0})[$] and [$]\frac{d^3p_{m0}}{d{B}^3}[$] by further differentiating the above function
We want to calculate the transition CDF, transition probability, and its derivatives from a cell m at time t, to an arbitrary point n on the grid at time t+1.  
We first calculate the CDF from entire mth cell to nth point on grid at time t+1 and write it as an integral as
[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$] Equation(A)
Here [$]p_m(B)[$] indicates the probability density function at arbitrary point B within the mth subdivision(not necessarily at center).
Here [$]P_{m,n}(B)[$] indicates the CDF of transition probability from an arbitrary point B within the mth subdivision(not necessarily from its center) to arbitrary point [$]B_n[$] on next grid.
Our strategy to solve the above integral is to take derivatives of integrand at center of the grid cell [$]B_{m0}[$] and expand it as Taylor series and then solve it. Our expansion is basically constant values calculated at center of the mth subdivision(which are handily available) multiplied by an integral of the powers of difference of distance of two boundaries from the center. I have already indicated previously how to calculate these derivatives at center to any higher order.
So we have 
[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$]
[$]=\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{dp_{m0}}{dB} P_{m0,n}+ p_{m0} p_{m0,n})  (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^2p_{m0}}{d{B}^2} P_{m0,n}+2 \frac{dp_{m0}}{dB} p_{m0,n}+  p_{m0} \frac{dp_{m0,n}}{dB})  \frac{{(B - B_{m0})}^2}{2}  dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^3p_{m0}}{d{B}^3} P_{m0,n}+3 \frac{d^2p_{m0}}{d{B}^2} p_{m0,n}+3 \frac{dp_{m0}}{dB} \frac{dp_{m0,n}}{dB}+  p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \frac{{(B - B_{m0})}^3}{6}   dB[$]  Equation(B)

We might have to expand to fourth order many times but I am not writing the higher order terms to make it easy to type here.
We re-arrange terms in the above equation(B) to write as

[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$]
[$]=\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{dp_{m0}}{dB} P_{m0,n})  (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^2p_{m0}}{d{B}^2} P_{m0,n})  \frac{{(B - B_{m0})}^2}{2}  dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^3p_{m0}}{d{B}^3} P_{m0,n})  \frac{{(B - B_{m0})}^3}{6}   dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  ( p_{m0} p_{m0,n})  (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (2 \frac{dp_{m0}}{dB} p_{m0,n}+  p_{m0} \frac{dp_{m0,n}}{dB})  \frac{{(B - B_{m0})}^2}{2}  dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (3 \frac{d^2p_{m0}}{d{B}^2} p_{m0,n}+3 \frac{dp_{m0}}{dB} \frac{dp_{m0,n}}{dB}+  p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \frac{{(B - B_{m0})}^3}{6}   dB[$] Equation(C)

We notice that first four terms in equation(C) after the equality sign and their series continuation equals the probability mass [$]\Delta P_m[$] in the mth cell(difference of CDF at boundaries) multiplied by transition CDF between the center of mth cell [$]B_{m0}[$] and point [$]B_n[$] on next time grid and we write those four terms again as

[$]\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{dp_{m0}}{dB} P_{m0,n})  (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^2p_{m0}}{d{B}^2} P_{m0,n})  \frac{{(B - B_{m0})}^2}{2}  dB[$]
[$]+\int_{B_{m1}}^{B_{m2}}  (\frac{d^3p_{m0}}{d{B}^3} P_{m0,n})  \frac{{(B - B_{m0})}^3}{6}   dB[$]
[$]=(P_{m2}-P_{m1}) P_{m0,n}=\Delta P_m P_{m0,n}[$] Equation(D)

We write the above results of Equation(D) and substitute them in equation (C) and take constants out of the integrals to write the equation(E) as 
[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$]
[$]=\Delta P_m P_{m0,n}[$]
[$]+ ( p_{m0} p_{m0,n})  \int_{B_{m1}}^{B_{m2}}  (B - B_{m0}) dB[$]
[$]+ (2 \frac{dp_{m0}}{dB} p_{m0,n}+  p_{m0} \frac{dp_{m0,n}}{dB}) \int_{B_{m1}}^{B_{m2}}  \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+(3 \frac{d^2p_{m0}}{d{B}^2} p_{m0,n}+3 \frac{dp_{m0}}{dB} \frac{dp_{m0,n}}{dB}+  p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \int_{B_{m1}}^{B_{m2}}   \frac{{(B - B_{m0})}^3}{6}  dB[$]  Equation(E)

Please note that first term on RHS is what we had been using earlier as if all mass were concentrated at delta center of mth subdivision. Next three terms are corrections to CDF required for an appropriate and exact calculation of CDF at point [$]B_n[$] on next grid when probability mass from entire cell is considered. Please notice that all derivatives evaluated at the center of mth grid subdivision can be trivially calculated as I explained earlier and we are left with extremely simple integrals of powers of distance of boundaries from the center of subdivision which are again extremely simple to calculate. 
Please note that all integrals with odd powers on [$](B - B_{m0}) [$] would have vanished if both grid boundaries were equidistant from the center. 
We can repeat the steps to calculate the corrections for transition PDF and all its derivatives just like we solved for the CDF. We will get the same equation but every one degree higher derivative of CDF will mean one degree higher derivative(of transition probability between center [$]B_{m0}[$] and the point [$]B_n[$] on next boundary) terms in equation(E)
To find corrections to transition CDF, we solved the integral [$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$] Equation(A)
In order to find corrections to transition PDF, we would solve the integral [$]\int_{B_{m1}}^{B_{m2}} p_m(B) p_{m,n}(B) dB[$] Equation(A2)
and so on for higher derivatives.
Sorry friends, it is a bit late(3:10 am) and I have to wake up early tomorrow so I will make any proper corrections of errors in this note tomorrow. I will remove "Post in progress" at start after I have removed any errors tomorrow.

Dear friends, I wanted to mention again that if friends would like to do more original research based on what I have posted here on this forum, please feel perfectly free to do that. Though I do expect that friends would cite whatever I have written on this forum when they do further research. 
I want to take this opportunity to thank those good people who protested against my mind control and who tried to press on mind control agencies to become better. Though mind control has not decreased and a lot of other problems continue, not forcing antipsychotic injections on me due to protests of good people is a very huge step for me and I would never have been able to do good research that I was recently able to do. I want to thank good people for making it possible for me to do good research and I want to express my huge gratitude to good people for favors they did to me. Though I hope that my mind control and mind control of thousands of innocent victims would eventually end, I really want to thank good people for making my life a lot easier and a whole lot better.
And I also want to tell friends that I seriously want to continue doing good research on stochastics and sharing it with friends on this thread. I will really try doing good research that friends like.

I want to request friends to read my new post on off-topic here:https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=864512#p864512

Post in progress.
I tried a bit to see if I could further improve on the formula for probability mass transfer CDF I wrote a few posts ago. I am copying equation(E) from the relevant earlier post. 
We write the above results of Equation(D) and substitute them in equation (C) and take constants out of the integrals to write the equation(E) as 
[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$]
[$]=\Delta P_m P_{m0,n}[$]
[$]+ ( p_{m0} p_{m0,n})  \int_{B_{m1}}^{B_{m2}}  (B - B_{m0}) dB[$]
[$]+ (2 \frac{dp_{m0}}{dB} p_{m0,n}+  p_{m0} \frac{dp_{m0,n}}{dB}) \int_{B_{m1}}^{B_{m2}}  \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+(3 \frac{d^2p_{m0}}{d{B}^2} p_{m0,n}+3 \frac{dp_{m0}}{dB} \frac{dp_{m0,n}}{dB}+  p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \int_{B_{m1}}^{B_{m2}}   \frac{{(B - B_{m0})}^3}{6}  dB[$]  Equation(E)

We have earlier simplified the equation by collecting some terms and recognizing their series and finding the term that describes the series. We want to try doing that with some other terms as well but it is a bit more difficult now.
Let us put separate some more terms and try to simplify them as

[$] ( p_{m0} p_{m0,n})  \int_{B_{m1}}^{B_{m2}}  (B - B_{m0}) dB[$]
[$]+ (p_{m0} \frac{dp_{m0,n}}{dB}) \int_{B_{m1}}^{B_{m2}}  \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+(p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \int_{B_{m1}}^{B_{m2}}   \frac{{(B - B_{m0})}^3}{6}  dB[$]  equation(a)

The above three terms could  be recognized as series for integral of CDF if we could add first order Taylor term that it is missing. So we could write as
[$]\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m,n}(B) dB[$]
[$] -  \int_{B_{m1}}^{B_{m2}} ( p_{m0} P_{m0,n})  dB[$]
[$] =  \int_{B_{m1}}^{B_{m2}}  ( p_{m0} p_{m0,n})(B - B_{m0}) dB[$]
[$]+  \int_{B_{m1}}^{B_{m2}} (p_{m0} \frac{dp_{m0,n}}{dB}) \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+ \int_{B_{m1}}^{B_{m2}} (p_{m0} \frac{d^2p_{m0,n}}{dB^3})   \frac{{(B - B_{m0})}^3}{6}  dB[$] equation(b)

So we can basically replace the series continuation of the following terms  
[$] ( p_{m0} p_{m0,n})  \int_{B_{m1}}^{B_{m2}}  (B - B_{m0}) dB[$]
[$]+ (p_{m0} \frac{dp_{m0,n}}{dB}) \int_{B_{m1}}^{B_{m2}}  \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+(p_{m0} \frac{d^2p_{m0,n}}{dB^3})  \int_{B_{m1}}^{B_{m2}}   \frac{{(B - B_{m0})}^3}{6}  dB[$] Equation(a)

with the following two terms

[$] p_{m0} \int_{B_{m1}}^{B_{m2}} P_{m,n}(B) dB[$]
[$] - ( p_{m0} P_{m0,n}) \int_{B_{m1}}^{B_{m2}}   dB[$] Equation(c)

The first analytic integral might be hard to solve analytically but we can try to integrate it by parts as 
[$]  \int_{B_{m1}}^{B_{m2}} P_{m,n} dB[$]
[$]=B_{m2} P_{m2,n} - B_{m1} P_{m1,n}  - \int_{B_{m1}}^{B_{m2}}  B p_{m,n}(B) dB[$]  Equation(d)

Substituting Equation(d) in Equation(c), we get
[$] p_{m0} \int_{B_{m1}}^{B_{m2}} P_{m,n}(B) dB[$]
[$] - ( p_{m0} P_{m0,n}) \int_{B_{m1}}^{B_{m2}}   dB[$]
[$]=p_{m0} B_{m2} P_{m2,n} -p_{m0} B_{m1} P_{m1,n} [$]
[$]- p_{m0} \int_{B_{m1}}^{B_{m2}}B p_{m,n}(B) dB[$]
[$] - ( p_{m0} P_{m0,n}) \int_{B_{m1}}^{B_{m2}}   dB[$] Equation(e)

Now we can substitute equation(e) in the original equa.tion(E) and write as
[$]\int_{B_{m1}}^{B_{m2}} p_m(B) P_{m,n}(B) dB[$]
[$]=\Delta P_m P_{m0,n}[$]
[$]+p_{m0} B_{m2} P_{m2,n} -p_{m0} B_{m1} P_{m1,n} [$]
[$]- p_{m0} \int_{B_{m1}}^{B_{m2}}B p_{m,n}(B) dB[$]
[$] - ( p_{m0} P_{m0,n}) \int_{B_{m1}}^{B_{m2}}   dB[$]
[$]+ (2 \frac{dp_{m0}}{dB} p_{m0,n}) \int_{B_{m1}}^{B_{m2}}  \frac{ {(B - B_{m0})}^2}{2} dB[$]
[$]+(3 \frac{d^2p_{m0}}{d{B}^2} p_{m0,n}+3 \frac{dp_{m0}}{dB} \frac{dp_{m0,n}}{dB})  \int_{B_{m1}}^{B_{m2}}   \frac{{(B - B_{m0})}^3}{6}  dB[$] 
[$] + higher order terms [$]