Serving the Quantitative Finance Community

 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 11th, 2021, 10:07 am

In order to speed up the program, it will be better to just project a trial point from the centre of grid cells at time t, calculate CDF at the trial point and then use CDF defect to calculate dB that would be added to the trial point so that new point is at right CDF. This can be done in parallel for the whole grid.
Secondly, you could turn CDF from all points farther than 6/6.5 SDs to one(or zero depending upon sign of transition Zt) and turn all their derivatives to zero(when two points are .6/6.5 SDs apart).
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 17th, 2021, 8:46 am

For my next project, I want to solve for the densities of SDEs in PDE framework and also for the densities of dt- and dz-integrals also in PDE framework both with an underlying Z-grid. I have rough ideas about how to do this and will keep working out details in next few days.
There is another thing I want to do simultaneously in the meantime and that is to write an algorithm for monte carlo simulation of SDEs with arbitrary number of drift terms and arbitrary number of diffusion terms. I am sure this is not very difficult for friends to calculate coefficients of an SDE with more than two drift terms and multiple diffusion terms but just that calculating all the coefficients by hand and keeping track of everything is very laborious so I thought that I will write an automated algorithm that will do all the calculations (up to fourth order or eighth order) especially since I have already written an algorithm to do this for SDEs with two(or three including QDs) drift coefficients so I already have good experience in this domain.
I hope my detention ends soon and I start working with these new projects.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 19th, 2021, 7:34 am

Friends, In the program that I posted about a week ago, I noticed that drift is not very properly accounted for and therefore with difficult drifts, you have to take an extremely small step size to keep the results accurate. I have re-done my analytics to include the precise effect of variable drift and variable volatility over the originating grid cell in the calculation of CDF and its derivatives on next time level. When I release the next version with these effects properly accounted for, I am sure that friends would be able to take a large step size of 1/16 or 1/8 for even the diffusions with difficult drifts and we would be able to do quantizer-like simulations with great precision and a large step size. I have done the analytics for it and I hope to make changes to previous program right after I am released from detention.
For variable vols in SDEs with original coordinates, I will also add the higher order effect of second hermite polynomials so that our simulations with this method remain very accurate.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 20th, 2021, 12:51 pm

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.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 21st, 2021, 3:43 pm

Today weather was very pleasant and I did not want to do programming so I thought I will play with equations and try to see if I can improve my previous work on transition probabilities where I had calculated transition normals with hermite polynomial weighted linear equations and it was working well for CEV type noises without drift.
I wrote some equations that I want to share with friends though it will take some weeks to develop the method fully.
Let me state the problem.
Suppose we have a SDE in bessel/Lamperti coordinates and we call it [$]B_1(t)[$] or simply [$]B_1[$]. Its density is closely related to underlying normal [$]Z_1[$]'s density through change of variable for densities as

[$]p(B_1)=p(Z_1)\frac{dZ_1}{dB_1}[$]

suppose we add some gaussian noise (possibly [$]Z_{12}=Z\sqrt{dt}[$]) to this bessel desnity and we get a new density [$]p(B_2)[$] which can again be written in terms of another underlying normal as

[$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$]

suppose standard deviation of [$]Z_1[$] is given by [$]\sigma_1[$] and that of [$]Z_2[$] is given by [$]\sigma_2[$] while standard deviation of [$]Z_{12}[$] is given by [$]\sigma_{12}[$]

Let me recall some facts about addition of two independent gaussians through convolution. We are adding as

[$]Z_2= Z_{12} + Z_{1}[$]
and convolution equation for that would be (read here: https://en.wikipedia.org/wiki/Sum_of_no ... _variables)

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2 - Z_1)}{\sigma_{12}})^2] \frac{1}{\sigma_{1}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_1}{\sigma_{1}})^2]  dZ_1 [$]

After several manipulations we have an equation of the form

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2] \int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1 -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] dZ_1[$]

The part under integration in above equation integrates to unity and we get the formula for normal density [$]Z_2[$] with standard deviation [$]\sigma_2[$] as

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2][$]


Ok now we come back to actual business. And I repeat the introduction again as

Suppose we have a SDE in bessel/Lamperti coordinates and we call it [$]B_1(t)[$] or simply [$]B_1[$]. Its density is closely related to underlying normal [$]Z_1[$]'s density through change of variable for densities as

[$]p(B_1)=p(Z_1)\frac{dZ_1}{dB_1}[$]

suppose we add some gaussian noise (possibly [$]Z_{12}=Z\sqrt{dt}[$]) to this bessel desnity and we get a new density [$]p(B_2)[$] which can again be written in terms of another underlying normal as

[$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$]


suppose we know [$]p(B_1)[$] and [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that

[$]B_2= Z_{12} + B_{1}[$]

Again both [$]B_1[$] and [$]B_2[$] are non-normal Bessel densities and we are adding a gaussian with certain variance to [$]B_1[$] in order to get density [$]B_2[$]. Again we know [$]p(B_1)[$] and related [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that [$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$] and we get [$]B_2[$] by adding [$]Z_{12}[$] to [$]B_1[$].

Our convolution equation now will be

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2 - Z_1(B_1))}{\sigma_{12}})^2] \frac{1}{\sigma_{1}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_1(B_1)}{\sigma_{1}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

We do exactly the same manipulations that we did for gaussian convolution and we have an equation of the form

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2] \int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

So we get a gaussian density [$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2)}{\sigma_{2}})^2][$] multiplied by 

[$]\int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

In pure Gaussian convolution, the above integral was becoming unity but now with weighted integration, this integral will no longer go to unity and will take a different value which should be equal to  [$]\frac{dZ_2}{dB_2}[$]
so we can write 
[$]\frac{dZ_2}{dB_2}(Z_2)=\int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]
where [$]\frac{dZ_2}{dB_2}(Z_2)[$] means  [$]\frac{dZ_2}{dB_2}[$] calculated at a particular value of [$]Z_2[$]
The above formula could be used in transition probabilities framework to caclulate the Bessel SDE variable at time t+1 when a gaussian is added to already existing bessel SDE variable especially when variance is changing all the time due to non-linear drift. 
We would have to repeat the integration for every point [$]Z_2[$] on the new grid therefore any practical algorithm would change integrations to summations over grid cells and some real work would have to be done before the idea could be successfully numerically implemented but I decided to share it with friends. Please pardon any errors.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 22nd, 2021, 9:26 am

I am writing this post to express my concern that mind control agencies want to put me on antipsychotic tablets so that mind control chemicals could be adjusted everyday in the tablets according to my changing neurotransmitters and I am sure they would move the psychiatrist to do precisely the same thing. If I have to take medication, I would much rather prefer injections since they cannot be altered everyday. There is absolutely no guarantee that I will ever get proper antipsychotic tablets, as they always had mind control chemicals in them and my mother kept the medicine on her so that it could be altered whenever mind control agencies wished. Though psychiatrist only mildly insisted on tablets, they continue to keep me detained and I am very afraid that finally they went to say that we have to prescribe tablets for you. Please stop mind control agencies from their ugly inhuman ploys to retard good-natured human beings because of their loathing of my religion at birth.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 22nd, 2021, 9:46 am

Even though psychiatrist today told that I would be getting my second half antipsychotic injection shot today, I am very afraid that eventually they want to put me on tablets after blaming some possible problem with injections and that is why they want to keep detaining me even after this second injection shot despite that seventeen days have passed since I was detained. Please stop mind control agency from their machinations.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 22nd, 2021, 10:47 am

My sister has tested positive for corona. She continued to come in contact with me several times when she arrived here to give me food. I have continued to request the psychiatrist since this morning to take my corona test so before they give me any antipsychotic injection (that can have a huge negative effect on immunity), they must take my corona test and they results would arrive in a few hours and once I am clear from corona I will take the injection but doctor has not listened yet. Let us see what the doctor says. 
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 23rd, 2021, 3:21 am

Today weather was very pleasant and I did not want to do programming so I thought I will play with equations and try to see if I can improve my previous work on transition probabilities where I had calculated transition normals with hermite polynomial weighted linear equations and it was working well for CEV type noises without drift.
I wrote some equations that I want to share with friends though it will take some weeks to develop the method fully.
Let me state the problem.
Suppose we have a SDE in bessel/Lamperti coordinates and we call it [$]B_1(t)[$] or simply [$]B_1[$]. Its density is closely related to underlying normal [$]Z_1[$]'s density through change of variable for densities as

[$]p(B_1)=p(Z_1)\frac{dZ_1}{dB_1}[$]

suppose we add some gaussian noise (possibly [$]Z_{12}=Z\sqrt{dt}[$]) to this bessel desnity and we get a new density [$]p(B_2)[$] which can again be written in terms of another underlying normal as

[$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$]

suppose standard deviation of [$]Z_1[$] is given by [$]\sigma_1[$] and that of [$]Z_2[$] is given by [$]\sigma_2[$] while standard deviation of [$]Z_{12}[$] is given by [$]\sigma_{12}[$]

Let me recall some facts about addition of two independent gaussians through convolution. We are adding as

[$]Z_2= Z_{12} + Z_{1}[$]
and convolution equation for that would be (read here: https://en.wikipedia.org/wiki/Sum_of_no ... _variables)

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2 - Z_1)}{\sigma_{12}})^2] \frac{1}{\sigma_{1}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_1}{\sigma_{1}})^2]  dZ_1 [$]

After several manipulations we have an equation of the form

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2] \int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1 -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] dZ_1[$]

The part under integration in above equation integrates to unity and we get the formula for normal density [$]Z_2[$] with standard deviation [$]\sigma_2[$] as

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2][$]


Ok now we come back to actual business. And I repeat the introduction again as

Suppose we have a SDE in bessel/Lamperti coordinates and we call it [$]B_1(t)[$] or simply [$]B_1[$]. Its density is closely related to underlying normal [$]Z_1[$]'s density through change of variable for densities as

[$]p(B_1)=p(Z_1)\frac{dZ_1}{dB_1}[$]

suppose we add some gaussian noise (possibly [$]Z_{12}=Z\sqrt{dt}[$]) to this bessel desnity and we get a new density [$]p(B_2)[$] which can again be written in terms of another underlying normal as

[$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$]


suppose we know [$]p(B_1)[$] and [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that

[$]B_2= Z_{12} + B_{1}[$]

Again both [$]B_1[$] and [$]B_2[$] are non-normal Bessel densities and we are adding a gaussian with certain variance to [$]B_1[$] in order to get density [$]B_2[$]. Again we know [$]p(B_1)[$] and related [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that [$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$] and we get [$]B_2[$] by adding [$]Z_{12}[$] to [$]B_1[$].

Our convolution equation now will be

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2 - Z_1(B_1))}{\sigma_{12}})^2] \frac{1}{\sigma_{1}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_1(B_1)}{\sigma_{1}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

We do exactly the same manipulations that we did for gaussian convolution and we have an equation of the form

[$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{Z_2}{\sigma_{2}})^2] \int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

So we get a gaussian density [$]\frac{1}{\sigma_{2}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2)}{\sigma_{2}})^2][$] multiplied by 

[$]\int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

In pure Gaussian convolution, the above integral was becoming unity but now with weighted integration, this integral will no longer go to unity and will take a different value which should be equal to  [$]\frac{dZ_2}{dB_2}[$]
so we can write 
[$]\frac{dZ_2}{dB_2}(Z_2)=\int_{-\infty}^{+\infty} \frac{1}{\frac{\sigma_{12} \sigma_1}{\sigma_2} \sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1) -\frac{{\sigma_1}^2}{{\sigma_2}^2} Z_2 )}{\frac{\sigma_{12} \sigma_1}{\sigma_2}})^2] \frac{dZ_1}{dB_1} dB_1 [$]
where [$]\frac{dZ_2}{dB_2}(Z_2)[$] means  [$]\frac{dZ_2}{dB_2}[$] calculated at a particular value of [$]Z_2[$]
The above formula could be used in transition probabilities framework to caclulate the Bessel SDE variable at time t+1 when a gaussian is added to already existing bessel SDE variable especially when variance is changing all the time due to non-linear drift. 
We would have to repeat the integration for every point [$]Z_2[$] on the new grid therefore any practical algorithm would change integrations to summations over grid cells and some real work would have to be done before the idea could be successfully numerically implemented but I decided to share it with friends. Please pardon any errors.
Sorry friends, this post has serious errors. I was not being very careful. Convolution with [$]B_1[$] would require
[$]p(B_2)=\int_{-\infty}^{+\infty} p(B_2-B_1)p(B_1) dB_1[$]
and [$]p(Z(B_1))[$] is not given as an explicit function of [$]B_1[$] and  [$]Z_2 - Z(B_1)[$] is different from [$]Z_2 - B_1[$]
Sorry friends for this error. I will try to find a way around this though as I have some ideas about it but I cannot be sure.
Again please pardon this mistake.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 23rd, 2021, 6:28 am

I had written some equations that I already shared with friends. Here is the corrected version
Let me state the problem.

Suppose we have a SDE in bessel/Lamperti coordinates and we call it [$]B_1(t)[$] or simply [$]B_1[$]. Its density is closely related to underlying normal [$]Z_1[$]'s density through change of variable for densities as

[$]p(B_1)=p(Z_1)\frac{dZ_1}{dB_1}[$]

suppose we add some gaussian noise (possibly [$]Z_{12}=Z\sqrt{dt}[$]) to this bessel desnity and we get a new density [$]p(B_2)[$] which can again be written in terms of another underlying normal as

[$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$]


suppose we know [$]p(B_1)[$] and [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that

[$]B_2= Z_{12} + B_{1}[$]

Again both [$]B_1[$] and [$]B_2[$] are non-normal Bessel densities and we are adding a gaussian with certain variance to [$]B_1[$] in order to get density [$]B_2[$]. Again we know [$]p(B_1)[$] and related [$]p(Z_1)[$] and we have calculated [$]\frac{dZ_1}{dB_1}[$], now suppose we are also given [$]p(Z_2)[$] and we want to calculate [$]\frac{dZ_2}{dB_2}[$] analytically given that [$]p(B_2)=p(Z_2)\frac{dZ_2}{dB_2}[$] and we get [$]B_2[$] by adding [$]Z_{12}[$] to [$]B_1[$].

Our convolution equation now will be (I have given the resultant variable name [$]Z_2[$] but it is equivalent to [$]B_2[$] as name does not matter for analytics(I just wanted to recognize the gaussian so gave the name [$]Z_2[$]).Please note that in below [$]Z_1(B_1)[$] means that [$]Z_1[$] is given as a function of [$]B_1[$]. Brackets might be confusing so I thought I will give some orientation.

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2 - B_1)}{\sigma_{12}})^2] \frac{1}{\sigma_{1}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_1(B_1))}{\sigma_{1}})^2] \frac{dZ_1}{dB_1} dB_1 [$]

It is straightforward to do some manipulations to get an equation of the form

[$]\frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2)}{\sigma_{12}})^2] \int_{-\infty}^{+\infty} \frac{1}{\sigma_1 \sqrt{(2 \pi)}} \exp[-.5 (\frac{(B_1^2 -2 B_1 Z_2)}{{\sigma_{12}}^2})] \exp[-.5 (\frac{Z_1(B_1)^2 }{{\sigma_1}^2})] \frac{dZ_1}{dB_1} dB_1 [$]

So we get a gaussian density [$]\frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2)}{\sigma_{12}})^2][$] multiplied by 

[$]\int_{-\infty}^{+\infty} \frac{1}{\sigma_1 \sqrt{(2 \pi)}} \exp[-.5 (\frac{(B_1^2 -2 B_1 Z_2)}{{\sigma_{12}}^2})] \exp[-.5 (\frac{Z_1(B_1)^2 }{{\sigma_1}^2})] \frac{dZ_1}{dB_1} dB_1 [$]


In pure Gaussian convolution, the above integral was becoming unity but now with weighted integration, this integral will no longer go to unity and will take a different value which should be equal to  [$]\frac{dZ_2}{dB_2}[$]
so we can write 
[$]\frac{dZ_2}{dB_2}=\int_{-\infty}^{+\infty} \frac{1}{\sigma_1 \sqrt{(2 \pi)}} \exp[-.5 (\frac{(B_1^2 -2 B_1 Z_2)}{{\sigma_{12}}^2})] \exp[-.5 (\frac{Z_1(B_1)^2 }{{\sigma_1}^2})] \frac{dZ_1}{dB_1} dB_1 [$]


The above formula could be used in transition probabilities framework to calculate the Bessel SDE variable at time t+1 when a gaussian is added to already existing bessel SDE variable especially when variance is changing all the time due to non-linear drift. Any practical algorithm would change integrations to summations over grid cells and some real work would have to be done before the idea could be successfully numerically implemented.
The above integral can easily be solved by integrating over Taylor series and all the coefficients(except one term that can easily be handled) can be solved once on the grid and then summed for vatrious values of [$]Z_2[$] or [$]B_2[$] 

In the gaussian density we get[$]\frac{1}{\sigma_{12}\sqrt{(2 \pi)}} \exp[-.5 (\frac{(Z_2)}{\sigma_{12}})^2][$] we can easily convert it to unit variance by change of variable but same scaling will also have to be applied in the term which is being integrated.

Please pardon the mistakes I made in the previous post. Though I have done decent work on this forum, sometimes I cannot see what is otherwise totally obvious to people or may be the drugs are starting to have an effect on me again.

I had previously written a function that used hermite polynomials with linear system of equations to make sure that branching normal to next transition probabilities grid has exact mean, variance, skew and kurtosis and used that with CEV noises. But with above formula we can calculate resulting normal on SDE in bessel coordinates grid with drifts that change variance and then use that normal to create a branching normal even for SDEs with drift. The program for system of linear equations that ensured perfect normality for branching normals was a good one but it was not written with an eye to speed and I have been optimizing it for speed now so that it can be used in real time for SDEs with variance changing drift since pre-calculated normal grid cannot be used there. I want to use somewhere 24 to 48 steps per year depending upon accuracy required. I hope that things go well. 
 
User avatar
Cuchulainn
Posts: 64433
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 23rd, 2021, 12:06 pm

Learn LATEX and make it a pleasure to read.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 24th, 2021, 11:18 am

I feel a bit goofy writing this post. Friends recall this latest program in which I used gaussian CDF and its derivatives to find grids where the probability mass in each grid point corresponded to probability mass of underlying Z_grid after solving the quartic equation. 
But it seems that it is probably a very expensive and relatively inefficient way to do this. We can do this extremely inexpensively by proposing a "guess target grid" and finding only CDF at that grid(not any derivatives or any solution of quartic equation) and once we know the CDF on "guessed target grid", we will know the corresponding Z_values from inverting the CDF. And then we can simply interpolate for the desired CDF points only by interpolation(as opposed to optimization like quartic formulas that also include calculation of derivatives) of the Z-values that underlie the desired CDF. So we just have to solve only for CDF on a "guessed target grid" and the find out corresponding Z-grid from inverting CDF and then interpolate against appropriate Z-values that underlie the desired CDF.
I hope that probability transfer program might be useful for some other related purpose.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 24th, 2021, 3:01 pm

I had previously written a function that used hermite polynomials with linear system of equations to make sure that branching normal to next transition probabilities grid has exact mean, variance, skew and kurtosis and used that with CEV noises. But with above formula we can calculate resulting normal on SDE in bessel coordinates grid with drifts that change variance and then use that normal to create a branching normal even for SDEs with drift. The program for system of linear equations that ensured perfect normality for branching normals was a good one but it was not written with an eye to speed and I have been optimizing it for speed now so that it can be used in real time for SDEs with variance changing drift since pre-calculated normal grid cannot be used there. I want to use somewhere 24 to 48 steps per year depending upon accuracy required. I hope that things go well. 
I have decided to do write this transition probabilities program for SDEs where drift changes the variance with my Taylor based formula for calculation of CDF coming from an entire SDE grid cell. I will start working in another 3-4 days and I hope that I would complete it in 7-10 days. I will try to make it a good program that friends like.
I was earlier reluctant to follow this approach thinking it would be too inefficient.
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 26th, 2021, 1:43 am

I am sharing this vector(array) version of the program that calculates CDF to an array of target points on the next grid simultaneously in matlab. It is a vector version of the earlier program. It is just to give an idea to friends that fast programs are possible to do this because I have not carefully tested this program and it may have some silly programming mistake somewhere but I am sure you will be able to fix the error(if there is any) very easily if you know logic of the method since it is a simple program.
After I leave detention, I will write and share a vector version of this program(in the context of larger program that advances the SDE for several steps) with variable drift and variable volatility also adding second and more hermite polynomials for accuracy(I have already done the analytics for adding hermite polynomials with precision.)
Here is the vector program. I am very sure it would work well but you have to check it. It is just meant to give an idea.



function [Pmn] = CalculateCDFArraySignIndex(Bn,Mm,sigma0,Bm0,Bm1,Bm2,dPm,pm0,dpm0,d2pm0,d3pm0,d4pm0,d5pm0)

%This function calculates CDF Pmn, PDF pmn, first derivative of Pdf dpmn,
%second derivative of PDF d2pmn, and third derivative of PDF d3pmn at a
%grid of target points Bn on next time level grid. All these CDF, PDF and their
%derivatives are calculated from an entire subdivision Bm at time t grid.
%Bn is the array of target points at time t+1 and Mm is the fianl index of
%array Bn starting from one.
%Bm1 and Bm2 are boundaries of grid cell Bm and Bm0 is center of grid cell
%Bm. 
% sigma0 is the volatility associated with transition gaussian originating
% from the center of grid Bm given as Bm0.
%dPm is the integrate probability mass in mth originating subdivision.
%pm0 is the value of probability distribution at the center of mth cell at 
%Bm0.
%dpm0 is the value of first derivative of probability distribution at the
%center of mth subdivision at Bm0;
%and so on for d2pm0(second derivative of pdf at Bm0), d3pm0, d4pm0, d5pm0.
%pdf and all its derivatives at center of mth subdivsion Bm0 are input to
%the function.
%This program is not written with regard to efficiency at all and I will
%optimize the final version for efficiency.

%The following two loops assume that Bn array is arranged in an ascending
%order as it usually is. The lower two arrays find the lower and upper
%limit of indices such that target array index in between is withing -7 SD
%to + 7SD of the transition normal from the originating cell under question
%(change SDs appropriately if you wish).
mm1=1;
for mm=1:Mm
    if( (Bn(mm)-Bm0) <-7*sigma0)
        mm1=mm;
    end
end

mm2=Mm;
for mm=Mm:-1:1
    if( (Bn(mm)-Bm0) >7*sigma0)
        mm2=mm;
    end
end


Bt(mm1:mm2)=Bn(mm1:mm2)-Bm0;
Zt(mm1:mm2)=(Bn(mm1:mm2)-Bm0)/sigma0;
Zt2(mm1:mm2)=Zt(mm1:mm2).*Zt(mm1:mm2);
Zt3(mm1:mm2)=Zt2(mm1:mm2).*Zt(mm1:mm2);
Zt4(mm1:mm2)=Zt3(mm1:mm2).*Zt(mm1:mm2);


cc1=sigma0;
cc2=cc1*sigma0;
cc3=cc2*sigma0;
cc4=cc3*sigma0;

Pm0n(mm1:mm2)=normcdf(Bt(mm1:mm2),0,cc1);
dPm0ndZt(mm1:mm2)=normpdf(Bt(mm1:mm2),0,cc1);
pm0n(mm1:mm2)=normpdf(Bt(mm1:mm2),0,cc1);
d2Pm0ndZt2(mm1:mm2)=-pm0n(mm1:mm2)*(Zt(mm1:mm2))./cc1;
d3Pm0ndZt3(mm1:mm2)=pm0n(mm1:mm2)*(Zt2(mm1:mm2)-1)./cc2;
d4Pm0ndZt4(mm1:mm2)=-pm0n(mm1:mm2)*(Zt3(mm1:mm2)-3*Zt(mm1:mm2))./cc3;
d5Pm0ndZt5(mm1:mm2)=pm0n(mm1:mm2)*(Zt4(mm1:mm2)-6*Zt2(mm1:mm2)+3)./cc4;

Zm1(mm1:mm2)=Bn(mm1:mm2)-(Bm1);
Zm2(mm1:mm2)=Bn(mm1:mm2)-(Bm2);

dBm=Bm2-Bm1;
dBm2=1/3* ((Bm2 - Bm0)^3 - (Bm1 - Bm0)^3)/2;
dBm3=1/4* ((Bm2 - Bm0)^4 - (Bm1 - Bm0)^4)/6;
dBm4=1/5* ((Bm2 - Bm0)^5 - (Bm1 - Bm0)^5)/24;
dBm5=1/6* ((Bm2 - Bm0)^6 - (Bm1 - Bm0)^6)/120;
dBm6=1/7* ((Bm2 - Bm0)^7 - (Bm1 - Bm0)^7)/720;

%Below are calculations for CDF at the target point Bn on next grid with respect to entire mth grid cell.
Integral1(mm1:mm2)=(0.398942* exp(-((0.5* Zm1(mm1:mm2).^2)/sigma0^2)) - 0.398942* exp(-((0.5* Zm2(mm1:mm2).^2)/sigma0^2)))*sigma0;%Integral of Zp(Z)dZ term in by parts integration in G2.
Integral2(mm1:mm2)=-Bn(1:Mm)*(normcdf(Zm2(mm1:mm2),0,sigma0)-normcdf(Zm1(mm1:mm2),0,sigma0));    %Integral of p(Z)dZ term multiplied by -Bn in G2. 
Integral3(mm1:mm2)=Bm2*normcdf(Zm2(mm1:mm2),0,sigma0)-Bm1*normcdf(Zm1(mm1:mm2),0,sigma0);  %The boundaries integral in by parts integration in G2.

Integral0(mm1:mm2)=dPm.*Pm0n(mm1:mm2); % G1 group of terms in wilmott explanation.

IntegralCDF(mm1:mm2)=Integral0(mm1:mm2)+ ...  %Integral0 is Group1(G1) in the explanation on wilmott.
     pm0.*(+Integral1(mm1:mm2)+Integral2(mm1:mm2)+Integral3(mm1:mm2))-  Pm0n(mm1:mm2)*pm0*dBm+ ... %This line is group2 or G2 in explanation.
     + 2* dpm0 * (dPm0ndZt(mm1:mm2).*-1) *dBm2 + ...  %This line and lines below are G3 in explanation.
    + (3*d2pm0*(dPm0ndZt(mm1:mm2).*-1) + 3*dpm0*(d2Pm0ndZt2(mm1:mm2)))*dBm3 + ...
    + (4*d3pm0*(dPm0ndZt(mm1:mm2).*-1) + 6* d2pm0*(d2Pm0ndZt2(mm1:mm2)) + 4* dpm0 * (d3Pm0ndZt3(mm1:mm2)*-1))*dBm4 + ...
    + (5*d4pm0*(dPm0ndZt(mm1:mm2).*-1) + 10* d3pm0*(d2Pm0ndZt2(mm1:mm2))+ 10* d2pm0*(d3Pm0ndZt3(mm1:mm2).*-1) + 5* dpm0 * (d4Pm0ndZt4(mm1:mm2)))*dBm5 + ...
    + (6*d5pm0*(dPm0ndZt(mm1:mm2).*-1) + 15* d4pm0*(d2Pm0ndZt2(mm1:mm2))+ 20* d3pm0*(d3Pm0ndZt3(mm1:mm2).*-1) + 15* d2pm0 * (d4Pm0ndZt4(mm1:mm2))+ ...
    6* dpm0 * (d5Pm0ndZt5(mm1:mm2)*-1)) * dBm6;

Pmn(mm1:mm2)=IntegralCDF(mm1:mm2);
if(mm1>1)
Pmn(1:mm1-1)=0;
end
if (mm2<Mm)
Pmn(mm2+1:Mm)=1;
end
end
 
User avatar
Amin
Topic Author
Posts: 2728
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

March 26th, 2021, 3:14 am

It seems that there is no plan to end my detention anytime soon. Yesterday, I talked to my sister on phone and asked when they want to end my detention and she was trying to evade the questions saying things like, "God willing, you will come home soon." and ""you will return home eventually." and all sort of similar things.
Psychiatrists and everyone clearly knows who is right and wrong but they continue to play games and pretend as if they do not know anything. Good people everywhere know the truth and I am very sure they are making strong efforts to rescue me from this but forces of evil are not ready to yield. Please protest against my detention and try to do something to end it if you can.