Numerical simulation for stochastic differential equations

stochastic calculus
A down-to-the-bone post, where few issues concerning the numerical simulation of stochastic differential equations are discussed with just a limited amount of technical detail.
Author

Angelo Maria Sabatini

Published

June 6, 2024

We consider here the class of stochastic processes known as diffusion processes, which are solutions to stochastic differential equation (SDE) of the form:

\[ dX(t)=b(t,X(t))\,dt+\sigma(t,X(t))\,dW(t) \tag{1}\]

with some initial condition \(X(0)\) which can be regarded as being either constant or random. The SDE is interpreted in the Itô sense:

\[ X(t)=X(0)+\int_0^tb(s,X(s))\,ds+\int_0^t\sigma(s,X(s))\,dW(s) \tag{2}\]

and \(\{W(t),t\in[0,T]\}\) is the so-called Brownian motion or Wiener process.

The Wiener process

The Wiener process \(\{W(t),t\in[0,T]\}\) has independent Gaussian increments and continuous sample paths, and comes endowed with the following properties:

  1. \(W(0)=0\) (with probability 1).

  2. For \(0\leq s\leq t\leq T\) the random variable given by the Gaussian increment \(W(t)-W(s)\) is distributed with mean zero and variance \(t-s\); equivalently, \(W(t)-W(s)\propto\sqrt{t-s}\,Z\), where \(Z\) denotes a standard Gaussian random variable, with a mean of zero and a variance of one.

  3. For \(0\leq s\leq t\leq u\leq v\leq T\) the increments \(W(t)-W(s)\) and \(W(v)-W(u)\) are independent (condition of independence on any two disjoint intervals).

The two deterministic functions \(b(\cdot,\cdot)\) and \(\sigma(\cdot,\cdot)\) are called, respectively, the drift and the diffusion coefficients of the SDE. Under a number of assumptions regarding their properties (i.e., global Lipschitz and linear growth, see (Iacus 2008) for mathematical details) the SDE has a unique solution with continuous sample paths such that

\[ E\left[\int_0^T\vert\,X(t)\,\vert^2\,dt\right]<\infty \tag{3}\]

For the sake of simpler notation, \(X(t)\rightarrow X_t\) is used in the following, therefore Equation 1 and Equation 2 can be written:

\[ \left\{ \begin{split} dX_t&=b(t,X_t)\,dt+\sigma(t,X_t)\,dW_t\\ X_t&=X_0+\int_0^tb(s,X_s)\,ds+\int_0^t\sigma(s,X_s)\,dW_s \end{split} \right. \tag{4}\]

The Itô formula is an important tool of stochastic calculus. It can be regarded as the stochastic version of the Taylor expansion of a function \(f(t,X_t)\) stopped at the second order, where \(X_t\) is a diffusion process. If the function \(f(\cdot,\cdot)\) is a twice differentiable function on both arguments, we have:

\[ df(t,X_t)=f_t(t,X_t)\,dt+f_x(t,X_t)\,dX_t+\frac{1}{2}f_{xx}(t,X_t)(dX_t)^2 \tag{5}\]

where:

\[ f_t(t,x)=\dfrac{\partial}{\partial t}f(t,x),\quad f_x(t,x)=\dfrac{\partial f(t,x)}{\partial x},\quad f_{xx}(t,x)=\dfrac{\partial^2 f(t,x)}{\partial x^2} \]

Using Equation 1 and the Itô lemma, which states that \((dW_t)^2=dt\), we have:

\[ (dX_t)^2=\left[b(t,X_t)\,dt+\sigma(t,X_t)\,dW_t\right]^2=\sigma^2(t,X_t)\,dt+\text{O}(dt^{3/2}) \tag{6}\]

Equation 5 can thus be written:

\[ df(t,X_t)=\left[f_t(t,X_t)+f_x(t,X_t)b(t,X_t)+\frac{1}{2}f_{xx}(t,X_t)\right]dt+f_x(t,X_t)\sigma(t,X_t)\,dW_t \tag{7}\]

Approximation methods

To apply a numerical method to the generic SDE of Equation 1 over the time interval \([0,T]\), the time interval needs to be discretized (Higham 2001). Let \(\Delta t=T/M\) for some positive integer \(M\), and \(t_n=n\Delta t\). The numerical approximation to \(X(t_n)\) will be denoted \(X_n\) (\(n=1,\cdots,M\)).

The Euler-Maruyama method takes the form:

\[ \left\{ \begin{split} X_n&=X_{n-1}+b(X_{n-1})\Delta t+\sigma(X_{n-1})\Delta W_{n-1}\\ X_0&=x_0 \end{split} \right. \tag{8}\]

It is worth noting that in the case when \(\sigma(X_t)=0\), then the problem becomes deterministic: this approximation method reduces to the standard Euler method for the ordinary differential equation \(\dot{X}_t=b(t,X_t)\).

The Milstein’s method adds a correction to the stochastic increment to produce the solution to the SDE. The correction arises because the Taylor expansion must be modified in the case of Itô calculus. Truncating the Taylor expansion to the second order at an appropriate point yields:

\[ \left\{ \begin{split} X_n&=X_{n-1}+b(t_{n-1},X_{n-1})\Delta t+\sigma(t_{n-1},X_{n-1})\Delta W_{n-1}\\ &+\dfrac{1}{2}\sigma(X_{n-1})\sigma_x(X_{n-1})\left[(\Delta W_{n-1})^2-\Delta t\right]\\ X_0&=x_0 \end{split} \right. \tag{9}\]

Strong and weak convergence

The methods of approximation of the continuous solution to an SDE are classified according to their different properties. Mainly two criteria of optimality are used in the literature: the strong and the weak (orders of) convergence.

A time-discretized approximation \(X_\delta\) of a continuous-time process \(X\), with \(\Delta t\) the time increment of the discretization, is said to be of strong order of convergence \(\gamma\) to \(X\) if for any fixed time horizon \(T\) it holds true that

\[ E[\,\vert X_{\Delta t}(\tau)-X(\tau)\vert\,]\leq C\Delta t^\gamma,\quad\forall\Delta t<\Delta t_0 \tag{10}\]

for any \(\tau\in[0,T]\), with \(\Delta t_0>0\) and \(C\) a constant that does not depend on \(\Delta t\).

Along with the strong convergence, the weak convergence can also be defined. In the same conditions as above, \(X_{\Delta t}\) is said to converge weakly of order \(\gamma\) to \(X\) if for any fixed horizon \(T\) and any \(2(\beta+1)\) continuous differentiable function \(g\) of polynomial growth, it holds true that

\[ \vert\,E_g[X_{\Delta t}(\tau)]-E_g[X(\tau)]\,\vert\leq C\Delta t^\beta \tag{11}\]

for any \(\tau\in[0,T]\), with \(\Delta t_0>0\) and \(C\) a constant that does not depend on \(\Delta t\).

Whereas the strong order of convergence measures the rate at which the mean of the error decays as \(\Delta t\rightarrow0\), the weak order of convergence measures the rate of decay of the error of the means.

Methods of approximation of some order that strongly converge usually have a higher order of weak convergence. This is the case with the Euler-Maruyama method, which is strongly convergent of order \(\gamma=1/2\) and weakly convergent of order \(\beta=1\) (under some smoothness conditions on the coefficients of the SDE). The Milstein’s method has both weak and strong order of convergence, \(\gamma=\beta=1\), which is superior to the Euler–Maruyama method.

It is interesting noting that, when the diffusion coefficient does not depend on the state variable of the process, the Euler-Maruyama and Milstein’s methods coincide. This is one important case, therefore, when the Euler-Maruyama method is strongly convergent of order \(\gamma=1\).

Numerical test

Geometric Brownian motion

The geometric Brownian motion (GBM) is the solution to the SDE:

\[ \left\{ \begin{split} dX_t&=\lambda\,X_t\,dt+\mu\,X_t\,dW_t\\ X_0&=x_0 \end{split} \right. \tag{12}\]

where the parameter \(\lambda\) is interpreted as the constant interest rate and \(\mu>0\) as the volatility. For a GBM, the drift and diffusion coefficients are both linearly related to the state variable \(X\), namely \(b(t,X_t)=\lambda\,X_t\) and \(\sigma(X_t)=\mu\,X_t\).

The explicit solution of the SDE Equation 12 can be found by operating the following change of variable:

\[ \begin{split} Y_t=\log X_t\quad\rightarrow\quad dY_t&=\frac{1}{X_t}\,dX_t-\frac{1}{2X_t^2}\,(dX_t)^2\\ &=\left(\lambda-\frac{1}{2}\mu^2\right)dt+\mu\,dW_t \end{split} \tag{13}\]

We obtain:

\[ X_t=x_0\exp\left[\left(\lambda-\frac{1}{2}\mu^2\right)t+\mu\,W_t\right] \tag{14}\]

For the GBM, the following stochastic difference equations are obtained for the Euler-Maruyama method:

\[ \left\{ \begin{split} X^E_n&=X^E_{n-1}(1+\lambda\Delta t+\mu\Delta W_{n-1})\\ X^E_0&=x_0 \end{split} \right. \tag{15}\]

and for the Milstein method:

\[ \left\{ \begin{split} X^M_n&=X^M_{n-1}+\lambda X^M_{n-1}\Delta t+\mu X^M_{n-1}\Delta W_{n-1}\\ &+\dfrac{1}{2}\mu^2 X^M_{n-1}\left[(\Delta W_{n-1})^2-\Delta t\right]\\ &=X^M_{n-1}\left[1+\left(\lambda-\dfrac{1}{2}\mu^2\right)\Delta t+\mu\Delta W_{n-1}+\dfrac{1}{2}\mu^2(\Delta W_{n-1})^2\right]\\ X^M_0&=x_0 \end{split} \right. \tag{16}\]

Recall that \(\Delta W_{n-1}=\sqrt{\Delta t}\,Z\), where \(Z\propto N(0,1)\).

Equation 15 reads:

\[ \left\{ \begin{split} X^E_n&=X^E_{n-1}(1+\lambda\Delta t+\mu\sqrt{\Delta t}\,Z)\\ X^E_0&=x_0 \end{split} \right. \tag{17}\]

Equation 16 reads:

\[ \left\{ \begin{split} X^M_n&=X^M_{n-1}\left[1+\left(\lambda+\dfrac{1}{2}\mu^2(Z^2-1)\right)\Delta t+\mu\sqrt{\Delta t}\,Z\right]\\ X^M_0&=x_0\\ \end{split} \right. \tag{18}\]

Comparing the exact solution Equation 14 with Equation 18 shows how the Milstein method makes the Taylor expansion exact up to order \(\text{O}(\Delta t)\).

The following code will test the strong convergence of both methods when they are applied to the linear SDE Equation 12, with \(\lambda=2,\mu=1,x_0=1\). One thousand sample paths are simulated over the time interval \([0,1]\), for step sizes being integer multiples of \(2^{-9}\;(8,4,2,1)\). The order of convergence is tested at the end point of the chosen time interval. The results of the simulations are shown in Figure 1.

Code 
library(tidyverse)

set.seed(1792)
lambda  <- 2 
mu      <- 1
Xzero   <- 1
T       <- 1
N       <- 2^9
dt      <- 1/N
M       <- 1000
Xerr_EM <- matrix(0, nrow = M, ncol = 5)
Xerr_M  <- matrix(0, nrow = M, ncol = 5)

for (i in seq(1, M, by = 1)) {
  dW    <- rnorm(N, mean = 0, sd = sqrt(dt))
  W     <- cumsum(dW)
  Xtrue <- Xzero*exp((lambda - 0.5*mu^2)*T + mu*W[N])
  for (j in 1:5) {
    R  <- 2^(j-1)    
    Dt <- R*dt 
    L  <- N/R
    Xtemp_EM <- Xzero
    Xtemp_M  <- Xzero
    for (k in seq(1, L, by = 1)) {
      Winc     <- sum(dW[(R*(k-1)+1):(R*k)])
      Xtemp_EM <- Xtemp_EM*(1 + Dt*lambda + mu*Winc)
      Xtemp_M  <- Xtemp_M*(1 + Dt*lambda + 0.5*mu^2*(Winc^2 - Dt) + mu*Winc)
    }
    Xerr_EM[i, j] <- abs(Xtemp_EM - Xtrue)
    Xerr_M[i, j]  <- abs(Xtemp_M - Xtrue)
  }
}
Dtvals <- c(1, 2, 4, 8)*dt
Xm_EM  <- array(0, dim = 4)
Xm_M   <- array(0, dim = 4)
for (i in 1:4) {
  Xm_EM[i] <- mean(Xerr_EM[, i])
  Xm_M[i]  <- mean(Xerr_M[, i])
}
t   <- rep(Dtvals, by = 2)
g   <- rep(c("Euler-Maruyama", "Milstein"), each = 4)
e   <- c(Xm_EM, Xm_M)
e_t <- c(sqrt(t), t)
df <- data.frame(t = t, method = g, error = e, error_true = e_t)
Figure 1: Strong error plot for the Euler-Maruyama and Milstein’s methods; dashed lines give the appropriate reference slope for each method, \(\gamma=0.5\) (Euler-Maruyama) and \(\gamma=1\) (Milstein).

Lamperti transform

An interesting application of the Itô formula Equation 7 is now discussed in connection with a formulation of Equation 4 where the diffusion coefficient is not an explicit function of time, but only depends on the state variable \(X\):

\[ dX_t=b(t,X_t)\,dt+\sigma(X_t)\,dW_t \tag{19}\]

The Lamperti transform is defined as follows:

\[ Y_t=F(X_t)=\int_z^{X_t}\frac{1}{\sigma(s)}\,ds \tag{20}\]

where \(z\) is any arbitrary value in the state space of \(X\). We assume that the function \(F(\cdot)\) defines a one-to-one mapping from the state space of \(X\) to \(\mathbb{R}\). We have:

\[ F_t(X_t)=0\quad\quad F_x(X_t)=\dfrac{1}{\sigma(X_t)}\quad\quad F_{xx}(X_t)=-\dfrac{\sigma_x(X_t)}{\sigma^2(X_t)} \tag{21}\]

Applying the Itô formula Equation 5, we obtain:

\[ \left\{ \begin{split} dY_t&=\left[\dfrac{b(t,X_t))}{\sigma(X_t)}-\dfrac{1}{2}\sigma_x(X_t)\right]dt+dW_t\\ X_t&=F^{-1}(Y_t) \end{split} \right. \tag{22}\]

The Lamperti transform changes the generic SDE Equation 4 into another SDE with unitary diffusion coefficient. As an example of application of the Lamperti transform, we consider here a result on transformations of SDEs that is relevant for the topic of interest in this post.

First, we discretize the Lamperti transform of Equation 22 using the Euler-Maruyama method (see Equation 8):

\[ Y_n=Y_{n-1}+\left[\dfrac{b(t_{n-1},X_{n-1})}{\sigma(X_{n-1})}-\dfrac{1}{2}\sigma_x(X_{n-1})\right]\Delta t+\sqrt{\Delta t}\,Z \tag{23}\]

Second, we apply the Taylor expansion to the inverse transform \(X_t=G(Y_t)\), noting that:

\[ \begin{split} G_y(Y_t)&=\sigma(G(Y_t))\\ G_{yy}(Y_t)&=\sigma(G(Y_t))\sigma_x(G(Y_t)) \end{split} \tag{24}\]

We obtain:

\[ G(Y_{n})=G(Y_{n-1})+G_y(Y_{n-1})(Y_n-Y_{n-1})+\frac{1}{2}G_{yy}(Y_{n-1})(Y_n-Y_{n-1})^2 \tag{25}\]

and hence:

\[ \begin{split} X_n&=X_{n-1}+\left[b(t_{n-1},X_{n-1})-\dfrac{1}{2}\sigma(X_{n-1})\sigma_x(X_{n-1})\right]\Delta t\\ &+\sigma(X_{n-1})\sqrt{\Delta t}\,Z+\frac{1}{2}\sigma(X_{n-1})\sigma_x(X_{n-1})\Delta t\,Z^2 \end{split} \tag{26}\]

In conclusion, the Milstein method on the original process and the Euler-Maruyama method on the transformed process are equal up to and including the order \(\text{O}(\Delta t)\).

In general if a transformation such as the Lamperti transform eliminates the interactions between the state of the process and the increments of the Wiener process, the performance of the Euler-Maruyama method on the transformed SDE are expected to improve.

The logarithmic transformation \(F(X_t)=\log X_t\) is an example of nonlinear transformation that is capable of eliminating the interaction between the state \(X_t\) of the GBM process and the increments of the Wiener process \(dW_t\), as shown in Equation 13. Therefore the Euler-Maruyama method can be confidently applied to the transformed process; the Milstein method is then obtained by simply taking the Taylor expansion of the inverse transform.

References

Higham, Desmond J. 2001. “An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations.” SIAM Review 43 (3): 525–46. https://doi.org/10.1137/s0036144500378302.
Iacus, Stefano M. 2008. Simulation and Inference for Stochastic Differential Equations. Springer New York. https://doi.org/10.1007/978-0-387-75839-8.