SME integration

In order to integrate a stochastic equation:

\[dX=a(X,t)dt+b(X,t)dW\]

if one wants to be more sophisticated than Euler integration, one can use Milstein integration:

\[X_{i+1}=X_i+a(X_i,t_i)\Delta t_i+b(X_i,t_i)\Delta W_i+ \frac{1}{2}b(X_i,t_i)\frac{\partial}{\partial X}b(X_i,t_i)\left( (\Delta W_i)^2-\Delta t_i\right)\]

Vector Milstein

What if we are interested in a vector-valued equation:

\[d\vec{\rho}=\vec{a}(\vec{\rho},t)dt+\vec{b}(\vec{\rho},t)dW\]

The way to generalize the Milstein scheme (while still restricting ourselves to a scalar-valued Wiener process) is

\[\rho^\mu_{i+1}=\rho^\mu_i+a^\mu_i\Delta t_i+b^\mu_i\Delta W_i+ \frac{1}{2}b^\nu_i\partial_\nu b^\mu_i\left((\Delta W_i)^2 -\Delta t_i\right),\]

where I have adopted an index convention for vectors such that

\[\begin{split}\begin{align} \vec{\rho}&=\rho^\mu\hat{e}_\mu \\ a^\mu_i&=a^\mu(\vec{\rho}_i,t_i) \\ \partial_\nu&=\frac{\partial}{\partial\rho^\nu}, \end{align}\end{split}\]

and indices that appear in both upper and lower positions in the same term are implicitly summer over.

For \(b^\mu=G^\mu_\nu\rho^\nu+k_\nu\rho^\nu\rho^\mu\) as defined in Vectorization we can write:

\[\begin{split}\begin{align} b^\nu\partial_\nu b^\mu&=\left(k_\nu G^\nu_\sigma\rho^\mu+ G^\mu_\nu G^\nu_\sigma+2k_\nu\rho^\nu(G^\mu_\sigma +k_\sigma\rho^\mu)\right)\rho^\sigma \\ b^\nu\partial_\nu b^\mu\hat{e}_\mu&=\left( \left(\vec{k}^\mathsf{T}G\vec{\rho}\right) +G^2+2(\vec{k}\cdot\vec{\rho})\left(G+\vec{k}\cdot \vec{\rho}\right)\right)\vec{\rho} \end{align}\end{split}\]

Order 1.5 Taylor scheme

To get higher order convergence in time, we can use a more complicated update formula (restricting ourselves to \(\vec{a}\) and \(\vec{b}\) with no explicit time dependence, as we have in our problem):

\[\begin{split}\begin{align} \rho^\mu_{i+1}&=\rho^\mu_i+a^\mu_i\Delta t_i+ b^\mu_i\Delta W_i+\frac{1}{2}b^\nu_i\partial_\nu b^\mu_i\left( (\Delta W_i)^2-\Delta t_i\right)+ \\ &\quad b^\nu_i\partial_\nu a^\mu_i\Delta Z_i +\left(a^\nu_i\partial_\nu +\frac{1}{2}b^\nu_ib^\sigma_i\partial_\nu\partial_\sigma\right) b^\mu_i\left(\Delta W_i\Delta t_i-\Delta Z_i\right)+ \\ &\quad\frac{1}{2}\left(a^\nu_i\partial_\nu +\frac{1}{2}b^\nu_ib^\sigma_i\partial_\nu\partial_\sigma\right) a^\mu_i\Delta t_i^2 +\frac{1}{2}b^\nu_i\partial_\nu b^\sigma_i\partial_\sigma b^\mu_i\left( \frac{1}{3}(\Delta W_i)^2-\Delta t_i\right)\Delta W_i \end{align}\end{split}\]

Recall from Vectorization that:

\[\begin{split}\begin{align} \vec{a}(\vec{\rho})&=Q\vec{\rho} \\ Q&:=(N+1)D(\vec{c})+ND(\vec{c}^*)+E(M,\vec{c})+F(\vec{h}) \end{align}\end{split}\]

\(\Delta Z\) is a new random variable related to \(\Delta W\):

\[\begin{split}\begin{align} \Delta W_i&=U_{1,i}\sqrt{\Delta t_i} \\ \Delta Z_i&=\frac{1}{2}\left(U_{1,i}+\frac{1}{\sqrt{3}}U_{2,i}\right) \Delta t_i^{3/2} \end{align}\end{split}\]

where \(U_1\), \(U_2\) are normally distributed random variables with mean 0 and variance 1.

The new terms in the higher-order update formula are given below:

\[\begin{split}\begin{align} b^\nu\partial_\nu a^\mu\hat{e}_\mu&=QG\vec{\rho} +(\vec{k}\cdot\vec{\rho})Q\vec{\rho} \\ a^\nu\partial_\nu b^\mu\hat{e}_\mu&=GQ\vec{\rho}+ (\vec{k}\cdot\vec{\rho})Q\vec{\rho}+\left( \vec{k}^\mathsf{T}Q\vec{\rho}\right)\vec{\rho} \\ a^\nu\partial_\nu a^\mu\hat{e}_\mu&=Q^2\vec{\rho} \\ b^\nu\partial_\nu b^\sigma\partial_\sigma b^\mu\hat{e}_\mu&=G^3\vec{\rho} +3(\vec{k}\cdot\vec{\rho})G^2\vec{\rho}+ 3\left(\vec{k}^\mathsf{T}G\vec{\rho}+ 2(\vec{k}\cdot\vec{\rho})^2\right)G\vec{\rho}+ \\ &\quad\left(\vec{k}^\mathsf{T}G^2\vec{\rho}+6(\vec{k}\cdot\vec{\rho}) \vec{k}^\mathsf{T}G\vec{\rho} +6(\vec{k}\cdot\vec{\rho})^3\right)\vec{\rho} \\ b^\nu b^\sigma\partial_\nu\partial_\sigma b^\mu\hat{e}_\mu&=2\left( \vec{k}^\mathsf{T}G\vec{\rho}+(\vec{k}\cdot\vec{\rho})^2\right)\left( G\vec{\rho}+(\vec{k}\cdot\vec{\rho})\vec{\rho}\right) \\ b^\nu b^\sigma\partial_\nu\partial_\sigma a^\mu\hat{e}_\mu&=0 \end{align}\end{split}\]

We explore testing the convergence rates in Testing.