Dear friends, In this post, I will try to explain with a bit more clarity how to calculate the transition CDF.PDF, and their derivatives between a whole grid cell at time t to an arbitrary point at time t+1. I will use the previous notation but with some slight changes, write and explain equations and later(in the same post) post the function code to calculate those CDF, PDF and its derivatives. I will also post some(simple) code that I used to compare results of some analytic integrals in above formulas with numerical integration. Today, I will only post the above function and related code and post the whole working program with explanatory comments tomorrow. And yes, I will post the simplified version of solution to quartic equation today(that will be used by the main code that I will post tomorrow).

I am copying my old post for reference but will try to explain things again.

**Amin:**

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 [$]

EDIT: In some of the integrals below I have used dB and mostly I have used [$]dB_m[$] but both are the same. I forgot the subscript m in some integrals. But I hope that friends would still perfectly understand as I have given reasonable orientation explanation for most of the integrals.

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+1.

Let us fix the notation first. We have a Bessel grid and we indicate points in an arbitrary cell by [$]B_m[$](this is notation, we use when we integrate over the originating cell etc) 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_m/2[$] and [$]Z_m-\Delta Z_m/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_m/2[$] as [$]B_{m1}[$] and right boundary corresponding to [$]Z_m+\Delta Z_m/2[$] as [$]B_{m2}[$]

Please note that gaussian [$]Z_m[$] that underlies the originating bessel grid is different from the transition gaussian between any point on the originating grid cell [$]B_m[$] to some target point [$]B_n[$] on grid at time t+1 which is indicated with [$]Z_t[$]. By definition [$]Z_t=B_n(t+1)-B_m(t)[$]. Since in the following analyticss, we have to calculate derivatives of gaussian [$]Z_t[$] with respect to both [$]B_n[$] and [$]B_m[$], we denote the derivatives with respect to [$]B_n[$] with a simplified notation using apostrophe on the original variable and derivatives with respect to [$]B_m[$] with a fraction notation. Please note that
[$]\frac{\partial B_n}{\partial Z_t}=1[$] and
[$]\frac{\partial B_m}{\partial Z_t}=-1[$]
The transition probability between arbitrary point on mth cell [$]B_{m}[$] to arbitary point [$]B_n[$] on next grid is denoted as [$] p_{m,n}=N(B_n-B_{m},0,\sigma_m)[$].Its CDF is defined by uppercase of PDF as [$] P_{m,n}=\int_{-\infty}^{(B_n-B_{m})} N(B_n-B_{m},0,\sigma_m) dN [$]
When we expand around center of mth grid cell [$]B_{m0}[$], the above quantities change notation and 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 [$]
When we calculate kth derivative of CDF or its other derivatives (with respect to target [$]B_n[$]) further with respect to originating cell [$]B_m[$], we use the simple formula [$]\frac{\partial^k P_{m,n}}{\partial {B_m}^k}[$][$]=\frac{\partial^k P_{m,n}}{\partial {Z_t}^k} {(\frac{\partial Z_t}{\partial B_m})}^k[$][$]={P_{m,n}}^{(k)}(Zt) * {(-1)}^{k}[$] since [$]\frac{\partial Z_t}{\partial B_m}=-1[$]
Now we come to derivatives of PDF in the originating cell [$]B_m[$] with underlying unit gaussian [$]Z_m[$]. Please note that this gaussian [$]Z_m[$] underlying originating cell [$]B_m[$] is different from transition gaussian [$]Z_t[$] between originating cell [$]B_m(t)[$] and target cell [$]B_n(t+1)[$]. In the following paragraph, we are only concerned with probability distribution and its derivatives in the originating cell [$]B_m[$]. In general, This probability distribution in originating cell is not a gaussian probability but this PDF and its derivatives in mth cell are related to underlying gaussian probability [$]Z_m[$] by change of variable formula for distributions.
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. (I calculated these derivatives in the main body of the program since they have to be calculated just once for every transition to target points and calculating them in the function that calculates CDF and its derivatives from originating cell to target points was not appropriate. They are simply input to the function that calculates CDF and its derivatives to target points at t+1 from a particular cell at time t.)
Again please note the difference that [$]P_{m0}[$] and [$]p_{m0}[$] are CDF and its PDF in the originating cell(calculated at its center) while [$]P_{m0,n}[$] and [$]p_{m0,n}[$] are transition CDF and PDF respectively between center of originating cell and the target point [$]B_n[$] on next time grid.
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_m) P_{m,n}(B_m) dB_m[$] Equation(A)
Here [$]p_m(B_m)[$] indicates the probability density function at arbitrary point [$]B_m[$] within the mth subdivision(not necessarily at center).
Here [$]P_{m,n}(B_m)[$] indicates the CDF of transition probability from an arbitrary point [$]B_m[$] within the mth subdivision(not necessarily from its center) to arbitrary point [$]B_n[$] on next (time level) grid such that [$]Z_t=B_n - B_m[$].
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 (in the following equation B indicates [$]B_m[$]. I have omitted m in subscript in the expansion after the first line.)
[$]\int_{B_{m1}}^{B_{m2}} p_m(B_m) P_{m,n} dB_m[$]
[$]=\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} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}} (\frac{d^2p_{m0}}{d{B}^2} P_{m0,n}+2 \frac{dp_{m0}}{dB} \frac{dP_{m0,n}}{dB}+ p_{m0} \frac{d^2P_{m0,n}}{dB^2}) \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} \frac{dP_{m0,n}}{dB}+3 \frac{dp_{m0}}{dB} \frac{d^2P_{m0,n}}{dB^2}+ p_{m0} \frac{d^3P_{m0,n}}{dB^3}) \frac{{(B - B_{m0})}^3}{6} dB[$] Equation(B)
[$]+ higher order terms in Taylor series[$]
Please note again that in the above equation odd derivatives of transition CDF with respect to [$]B_m[$] or B (abbreviated notation) will take a negative sign with respect to regular(plain) transition CDF and its derivatives(with respect to [$]Z_t[$] as [$]\frac{dB_m}{dZ_t}=-1[$].
We rearrange the above equation in further three groups of terms. Then We find analytic approximation to two of those three groups and add simply add the third group as such.
[$]\int_{B_{m1}}^{B_{m2}} p_m(B_m) P_{m,n}(B_m) dB_m[$]
[$]=G1 + G2 + G3[$]
where

[$]G1=\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[$]
and
[$]G2=\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{d^2P_{m0,n}}{dB^2}) \frac{{(B - B_{m0})}^2}{2} dB[$]
[$]+\int_{B_{m1}}^{B_{m2}} (p_{m0} \frac{d^3P_{m0,n}}{dB^3}) \frac{{(B - B_{m0})}^3}{6} dB[$]
[$]G3=\int_{B_{m1}}^{B_{m2}} (2 \frac{dp_{m0}}{dB} \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} \frac{dP_{m0,n}}{dB}+3 \frac{dp_{m0}}{dB} \frac{d^2P_{m0,n}}{dB^2}) \frac{{(B - B_{m0})}^3}{6} dB[$] Equation(C)

We notice that after taking [$]P_{m0,n}[$] common out of terms in G!, the terms in G1(and their taylor series continuation) are integral of a simple Taylor expansion of probability density function in the originating mth Grid at time t and hence represent the integral of probability mass in mth subdivision. Hence we can write this group G1 as

[$]G1= P_{m0,n} \int_{B_{m1}}^{B_{m2}} p_{m}(B_m) dB_m= P_{m0,n} \Delta P_m[$]
where [$]\Delta P_m[$] indicates total (integral of) probability mass in the mth subdivision.
Now we take the second group of terms.
[$]G2=\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]
[$]+\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{d^2P_{m0,n}}{dB^2}) \frac{{(B - B_{m0})}^2}{2} dB[$]
[$]+\int_{B_{m1}}^{B_{m2}} (p_{m0} \frac{d^3P_{m0,n}}{dB^3}) \frac{{(B - B_{m0})}^3}{6} dB[$]

We notice that if we take [$] p_{m0}[$] common, this group of terms (and their taylor series continuation) would equal the (non-probability weighted) integral of transition CDF over the mth subdivision but in order to do that we will have to add one first term in the series that would equal [$]T2= \int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB[$]. Only if we could add T2 to this G2 group of terms, it would be equal to (non-probability weighted) integral of transition CDF over mth subdivision. So we can write
[$]G2= p_{m0} \int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m-\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB [$]
The term T2 we added over G2 is very simple but we would still have to see how to analytically solve (non-probability weighted) integral of CDF over mth subdivision given as [$]\int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m[$]. In order to solve this integral by parts as
[$]\int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m=\int_{B_{m1}}^{B_{m2}} d[B_m P_{m,n}(B_m)] - \int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]
[$]=B_{m2} P_{m2,n} - B_{m1} P_{m1,n} - \int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]
In order to solve the second term in last line, we write as
[$]\int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]
We notice that [$]B_m=B_n-Z_t[$] , [$]\frac{d P_{m,n}}{dB_m}= -p_{m,n}(Z_t)[$] and [$]dB_m=- dZ_t[$], [$]Z_{Bm1}=B_n - B_{m1}[$] and [$]Z_{Bm2}=B_n - B_{m2}[$]
so the integral in question becomes
[$]\int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m =\int_{Z_{Bm1}}^{Z_{Bm2}} ( B_n-Z_t) p_{m,n}(Z_t) dZ_t[$]
dropping the subscript t in the notation of transition gaussian, we write the above equation as
[$]=\int_{Z_{Bm1}}^{Z_{Bm2}} ( B_n-Z) p_{m,n}(Z) dZ= B_n \int_{Z_{Bm1}}^{Z_{Bm2}} p(Z) dZ - \int_{Z_{Bm1}}^{Z_{Bm2}} Z p(Z) dZ[$]
Both of the above integrals can be solved analytically. I checked the whole of this by parts analytic solution by comparing it with numerical integration of the(non-probability weighted) CDF over different subdivisions and found it to be the right formula. In the end today, I will also give some (informal) programs I used to compare this by parts integral with numerical integration.
With the above calculations, the terms in group G2 equal to the expression as
[$]G2= p_{m0} \int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m-\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB [$]
[$]= p_{m0} (B_{m2} P_{m2,n} - B_{m1} P_{m1,n} ) - p_{m0} B_n \int_{Z_{Bm1}}^{Z_{Bm2}} p(Z) dZ + p_{m0} \int_{Z_{Bm1}}^{Z_{Bm2}} Z p(Z) dZ[$]
[$]- p_{m0} P_{m0,n} \int_{B_{m1}}^{B_{m2}} dB [$]
Please note that gaussian Z in above is not a unit gaussian, and it takes variance associated with the transition probability.
We add the terms in group G3 individually without any special analytical treatment.
I am pasting the code in next post.

I have re-written the earlier note so that it accounts for variable drift and variable volatility over the originating cell for calculation of CDF from originating grid cell to a target point at time t+1 when the dynamics are governed by an SDE. I will first outline three differences with respect to previous exposition.

1. I have changed the definition of transition gaussian and defined it to be standard as opposed to previous exposition in which it was not defined to be standard.
[$]Z_t=\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{\sigma(B_m(t))}[$]

2. Earlier the derivatives of transition gaussian with respect to originating cell [$]B_m[$] were very simple due to different defintion, constant volatility and due to lack of drift but now it has changed to
[$]\frac{\partial Z_t}{\partial B_m}=\frac{-1-\frac{\partial\mu(B_m}{\partial B_m}}{\sigma(B_m(t))}-\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{{\sigma(B_m(t))}^2} \frac{\partial\sigma(B_m)}{\partial B_m}[$]

and so on for higher derivatives.

3. Later in the previous exposition, I divided the main Taylor expansion into three groups G1, G2 nd G3. In this new note, Math for groups G1 and G3 remains almost identical (except that [$]\frac{\partial Z_t}{\partial B_m}[$] has changed) but I have re-done math for Group G2 that required more care with variable drift and variable volatility.
I want to emphasize that [$]\mu(B_m(t))[$] and [$]\sigma(B_m(t))[$] both are expressions with anywhere between 5 to 25 terms and each term is algebraic in [$]B_m[$] multiplied by some coefficient so both of these can be differentiated many times at the center and both [$]\mu(B_m(t))[$] and [$]\sigma(B_m(t))[$] can be easily integrated over the originating grid cell due to their algebraic(multiplied by a coefficient depending on transition time and SDE parameters) nature.

Please note that in the note below, I have used subdivision and grid cell interchangeably and they mean the same thing.

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+1.

Let us fix the notation first. We have an originating grid at time [$]t_1[$] and we indicate points in an arbitrary cell on this grid by [$]B_m[$](this is notation, we use when we integrate over the originating cell etc) 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_m/2[$] and [$]Z_m-\Delta Z_m/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_m/2[$] as [$]B_{m1}[$] and right boundary corresponding to [$]Z_m+\Delta Z_m/2[$] as [$]B_{m2}[$]

Please note that gaussian [$]Z_m[$] that underlies the originating bessel grid is different from the transition gaussian between any point on the originating grid cell [$]B_m[$] to some target point [$]B_n[$] on grid at time t+1 which is indicated with [$]Z_t[$]. By definition [$]Z_t=\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{\sigma(B_m(t))}[$]. Please notice that in previous program, the transition gaussian was not standard by definition but now we take it as standard Gaussian by definition.

Since in the following analytics, we have to calculate derivatives of gaussian [$]Z_t[$] with respect to both [$]B_n[$] and [$]B_m[$], we denote the derivatives with respect to [$]B_n[$] with a simplified notation using apostrophe on the original variable and derivatives with respect to [$]B_m[$] with a fraction notation. Please note that

since

[$]Z_t=\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{\sigma(B_m(t))}[$]

[$]\frac{\partial Z_t}{\partial B_n}=\frac{1}{\sigma(B_m(t))}[$]

[$]\frac{\partial Z_t}{\partial B_m}=\frac{-1-\frac{\partial\mu(B_m}{\partial B_m}}{\sigma(B_m(t))}-\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{{\sigma(B_m(t))}^2} \frac{\partial\sigma(B_m)}{\partial B_m}[$]

We will also need higher drivatives [$]\frac{\partial^2 Z_t}{\partial {B_m}^2}[$] and so on that we will easily find with repeated differentiation. We will need these derivatives expressions at the center of the originating grid cell and since expressions for [$]\mu(B_m(t))[$] and [$]\sigma(B_m(t))[$] both can be easily differentiated, this is not difficult.

We will calculate the derivatives of transition probability or transition CDF with respect to [$]B_m[$] by using the equations

[$]\frac{dP_{m,n}}{dB_m}=\frac{dP_{m,n}}{dZ_t} \frac{\partial Z_t}{\partial B_m}[$]

[$]\frac{d^2P_{m,n}}{{dB_m}^2}=\frac{d^2P_{m,n}}{d{Z_t}^2} {(\frac{\partial Z_t}{\partial B_m})}^2+\frac{dP_{m,n}}{dZ_t} \frac{\partial^2 Z_t}{\partial {B_m}^2}[$]

and so on for higher order

The transition probability between arbitrary point on mth cell [$]B_{m}[$] to arbitary point [$]B_n[$] on next grid is denoted as [$] p_{m,n}=\frac{1}{\sigma(B_m)\sqrt{2 \pi}}Exp(-.5 {(\frac{(B_n(t+1)-B_m(t)-\mu(B_m(t)))}{\sigma(B_m(t))})}^2)[$].Its CDF is defined by uppercase of PDF as integral of above probability density of transition Zt from -infinity to Zt.
When we expand around center of mth grid cell [$]B_{m0}[$], the above quantities change notation and the transition probability between center of mth cell [$]B_{m0}[$] to arbitary point [$]B_n[$] on next grid is denoted as [$] p_{m0,n}=\frac{1}{\sigma(B_m)\sqrt{2 \pi}}Exp(-.5 {(\frac{(B_n(t+1)-B_m0(t)-\mu(B_m0(t)))}{\sigma(B_m0(t))})}^2)[$].Its CDF is denoted by uppercase of PDF as [$] P_{m0,n} [$]

Now we come to derivatives of PDF in the originating cell [$]B_m[$] with underlying unit gaussian [$]Z_m[$]. Please note that this gaussian [$]Z_m[$] underlying originating cell [$]B_m[$] is different from transition gaussian [$]Z_t[$] between originating cell [$]B_m(t)[$] and target cell [$]B_n(t+1)[$]. In the following paragraph, we are only concerned with probability distribution and its derivatives in the originating cell [$]B_m[$]. In general, This probability distribution in originating cell is not a gaussian probability but this PDF and its derivatives in mth cell are related to underlying gaussian probability [$]Z_m[$] by change of variable formula for distributions.

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. (I calculated these derivatives in the main body of the program since they have to be calculated just once for every transition to target points and calculating them in the function that calculates CDF and its derivatives from originating cell to target points was not appropriate. They are simply input to the function that calculates CDF and its derivatives to target points at t+1 from a particular cell at time t.)

Again please note the difference that [$]P_{m0}[$] and [$]p_{m0}[$] are CDF and its PDF in the originating cell(calculated at its center) while [$]P_{m0,n}[$] and [$]p_{m0,n}[$] are transition CDF and PDF respectively between center of originating cell and the target point [$]B_n[$] on next time grid.

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_m) P_{m,n}(B_m) dB_m[$] Equation(A)

Here [$]p_m(B_m)[$] indicates the probability density function at arbitrary point [$]B_m[$] within the mth subdivision(not necessarily at center).

Here [$]P_{m,n}(B_m)[$] indicates the CDF of transition probability from an arbitrary point [$]B_m[$] within the mth subdivision(not necessarily from its center) to arbitrary point [$]B_n[$] on next (time level) 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 (in the following equation B indicates [$]B_m[$]. I have omitted m in subscript in the expansion after the first line.)

[$]\int_{B_{m1}}^{B_{m2}} p_m(B_m) P_{m,n} dB_m[$]

[$]=\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} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]

[$]+\int_{B_{m1}}^{B_{m2}} (\frac{d^2p_{m0}}{d{B}^2} P_{m0,n}+2 \frac{dp_{m0}}{dB} \frac{dP_{m0,n}}{dB}+ p_{m0} \frac{d^2P_{m0,n}}{dB^2}) \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} \frac{dP_{m0,n}}{dB}+3 \frac{dp_{m0}}{dB} \frac{d^2P_{m0,n}}{dB^2}+ p_{m0} \frac{d^3P_{m0,n}}{dB^3}) \frac{{(B - B_{m0})}^3}{6} dB[$] Equation(B)

[$]+ higher order terms in Taylor series[$]

We rearrange the above equation in further three groups of terms. Then We find analytic approximation to two of those three groups and add simply add the third group as such.

[$]\int_{B_{m1}}^{B_{m2}} p_m(B_m) P_{m,n}(B_m) dB_m[$]

[$]=G1 + G2 + G3[$]

where

[$]G1=\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[$]

and

[$]G2=\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]

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

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

[$]G3=\int_{B_{m1}}^{B_{m2}} (2 \frac{dp_{m0}}{dB} \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} \frac{dP_{m0,n}}{dB}+3 \frac{dp_{m0}}{dB} \frac{d^2P_{m0,n}}{dB^2}) \frac{{(B - B_{m0})}^3}{6} dB[$] Equation(C)

We notice that after taking [$]P_{m0,n}[$] common out of terms in G1, the terms in G1(and their taylor series continuation) are integral of a simple Taylor expansion of probability density function in the originating mth Grid at time t and hence represent the integral of probability mass in mth subdivision. Hence we can write this group G1 as

[$]G1= P_{m0,n} \int_{B_{m1}}^{B_{m2}} p_{m}(B_m) dB_m= P_{m0,n} \Delta P_m[$]

where [$]\Delta P_m[$] indicates total (integral of) probability mass in the mth subdivision.

Now we take the second group of terms.

[$]G2=\int_{B_{m1}}^{B_{m2}} ( p_{m0} \frac{dP_{m0,n}}{dB}) (B - B_{m0}) dB[$]

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

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

We notice that if we take [$] p_{m0}[$] common, this group of terms (and their taylor series continuation) would equal the (non-probability weighted) integral of transition CDF over the mth subdivision but in order to do that we will have to add one first term in the series that would equal [$]T2= \int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB[$]. Only if we could add T2 to this G2 group of terms, it would be equal to (non-probability weighted) integral of transition CDF over mth subdivision. So we can write

[$]G2= p_{m0} \int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m-\int_{B_{m1}}^{B_{m2}} p_{m0} P_{m0,n} dB [$]

The term T2 we added over G2 is very simple but we would still have to see how to analytically solve (non-probability weighted) integral of CDF over mth subdivision given as [$]\int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m[$]. In order to solve this integral by parts as

[$]\int_{B_{m1}}^{B_{m2}} P_{m,n}(B_m) dB_m=\int_{B_{m1}}^{B_{m2}} d[B_m P_{m,n}(B_m)] - \int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]

[$]=B_{m2} P_{m2,n} - B_{m1} P_{m1,n} - \int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]

****In order to solve the second term in last line, we write as

[$]\int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m [$]

We notice that [$]B_m=B_n-\sigma(B_m) Z_t -\mu(B_m)[$]

so the integral in question becomes

[$]\int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m =\int_{B_{m1}}^{B_{m2}} (B_n-\sigma(B_m) Z_t -\mu(B_m) \frac{d P_{m,n}}{dB_m} dB_m[$]

so we have

[$]\int_{B_{m1}}^{B_{m2}} B_m \frac{d P_{m,n}}{dB_m} dB_m = I_1 + I_2 +I_3[$]

where

[$]I_1=\int_{B_{m1}}^{B_{m2}} B_n \frac{d P_{m,n}}{dB_m} dB_m[$]

[$]I_2=-\int_{B_{m1}}^{B_{m2}} \sigma(B_m) Z_t \frac{d P_{m,n}}{dB_m} dB_m[$]

[$]I_3=-\int_{B_{m1}}^{B_{m2}} \mu(B_m) \frac{d P_{m,n}}{dB_m} dB_m[$]

integral [$]I_1[$] is straightforward since [$]B_n[$] is a constant so it can be solved as

[$]I_1=B_n (P_{m2,n}-P_{m1,n})[$]

We want to carefully solve integral [$]I_2[$]

[$]I_2=-\int_{B_{m1}}^{B_{m2}} \sigma(B_m) Z_t \frac{d P_{m,n}}{dB_m} dB_m[$]

[$]=\int_{B_{m1}}^{B_{m2}} \sigma(B_m) \frac{d p_{m,n}}{dB_m} dB_m[$]

where I have used [$]Z_t \frac{d P_{m,n}}{dB_m}=Z_t \frac{d P_{m,n}}{dZ_t}\frac{d Z_t}{dB_m}[$]

[$]=-\frac{d p_{m,n}}{dZ_t}\frac{d Z_t}{dB_m}=\frac{d p_{m,n}}{dB_m}[$]

so

[$]I_2=\int_{B_{m1}}^{B_{m2}} \sigma(B_m) \frac{d p_{m,n}}{dB_m} dB_m[$]

Just like the main original integral, we expand this integral again into Taylor series at the centre of grid cell and then among the expansion, we identify two Taylor expansions at the ends both of which can be solved analytically just like we have done for the main integral.

So we have

[$]I_2=\int_{B_{m1}}^{B_{m2}} \sigma(B_m) \frac{d p_{m,n}}{dB_m} dB_m[$]

[$]=\sigma(B_{m0}) \frac{d p_{m0,n}}{dB_m} (B_{m1}-B_{m2})[$]

[$]+(\frac{d\sigma(B_{m0})}{dB_m} \frac{d p_{m0,n}}{dB_m} + \sigma(B_{m0}) \frac{d^2 p_{m0,n}}{d{B_m}^2}) \int_{B_{m1}}^{B_{m2}} (B-B_{m0}) dB_m[$]

[$]+(\frac{d^2\sigma(B_{m0})}{d{B_m}^2} \frac{d p_{m0,n}}{dB_m} +2 \frac{d\sigma(B_{m0})}{dB_m} \frac{d^2 p_{m0,n}}{d{B_m}^2}+ \sigma(B_{m0}) \frac{d^3 p_{m0,n}}{d{B_m}^3}) \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^2}{2} dB_m[$]

We again divide the above expansion into three groups of terms [$]H1[$], [$]H2[$], [$]H3[$] where

[$]H1=\sigma(B_{m0}) \frac{d p_{m0,n}}{dB_m} (B_{m1}-B_{m2})[$]

[$]+(\frac{d\sigma(B_{m0})}{dB_m} \frac{d p_{m0,n}}{dB_m}) \int_{B_{m1}}^{B_{m2}} (B-B_{m0}) dB_m[$]

[$]+(\frac{d^2\sigma(B_{m0})}{d{B_m}^2} \frac{d p_{m0,n}}{dB_m}) \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^2}{2} dB_m[$]

[$]+(\frac{d^3\sigma(B_{m0})}{d{B_m}^3} \frac{d p_{m0,n}}{dB_m}) \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^3}{6} dB_m[$]

We immediately identify that

[$]H1=\frac{d p_{m0,n}}{dB_m} \int_{B_{m1}}^{B_{m2}} \sigma(B_m)dB_m[$]

where [$]\frac{d p_{m0,n}}{dB_m}[$] is a constant term and we had already indicated that in our set up [$]\sigma(B_m)[$] means a set of anywhere between 5 to 20 algebraic terms(with constant coefficient multipliers) that can easily be integrated over the grid cell analytically.

Now we identify the second group of terms

[$]H2=( \sigma(B_{m0}) \frac{d^2 p_{m0,n}}{d{B_m}^2}) \int_{B_{m1}}^{B_{m2}} (B-B_{m0}) dB_m[$]

[$]+ \sigma(B_{m0}) \frac{d^3 p_{m0,n}}{d{B_m}^3} \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^2}{2} dB_m[$]

[$]+ \sigma(B_{m0}) \frac{d^4 p_{m0,n}}{d{B_m}^4} \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^3}{6} dB_m[$]

We notice that if we can take [$]\sigma(B_m0)[$] constant and add a first order term in taylor expansion, we can again identify it as a Taylor series of [$]\frac{d p_{m,n}}{dB_m}[$] which can again be analytically integrated since its values(of its integral) at both ends are known. So we have

[$]H2=\sigma(B_{m0}) \int_{B_{m1}}^{B_{m2}} \frac{d p_{m,n}}{dB_m} dB_m- ( \sigma(B_{m0}) \frac{d p_{m0,n}}{dB_m}) \int_{B_{m1}}^{B_{m2}} dB_m[$]

[$]=\sigma(B_{m0}) (p_{m2,n}-p_{m1,n})- ( \sigma(B_{m0}) \frac{d p_{m0,n}}{dB_m}) (B_{m2}- B_{m1})[$]

We identify the third group of terms as

[$]H3=2 \frac{d\sigma(B_{m0})}{dB_m} \frac{d^2 p_{m0,n}}{d{B_m}^2} \int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^2}{2} dB_m[$]

[$]+(3 \frac{d^2\sigma(B_{m0})}{d{B_m}^2} \frac{d^2 p_{m0,n}}{d{B_m}^2} +3 \frac{d\sigma(B_{m0})}{dB_m} \frac{d^3 p_{m0,n}}{d{B_m}^3})\int_{B_{m1}}^{B_{m2}} \frac{(B-B_{m0})^3}{6} dB_m[$]

+Taylor continuation.

We would have to add this group of terms as such one by one and there is no analytic solution for them so We add the terms in group [$]H3[$] individually without any special analytical treatment.

So we can easily solve integral [$]I_2=H1+H2+H3[$]
The integral [$]I_3[$] is very similar to integral [$]I_2[$] so I am not writing how to solve it again.
Coming back to group G3 in the main Taylor expansion, We add the terms in group [$]G3[$] individually without any special analytical treatment.

In a few days, I will be modifying the program to add the above changes.