Milstein method

In mathematics, the Milstein method is a technique for the approximate numerical solution of a stochastic differential equation. It is named after Grigori N. Milstein who first published it in 1974.

Description
Consider the autonomous Itō stochastic differential equation: $$\mathrm{d} X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t$$ with initial condition $$X_{0} = x_{0}$$, where $$W_{t}$$ stands for the Wiener process, and suppose that we wish to solve this SDE on some interval of time $$[0,T]$$. Then the Milstein approximation to the true solution $$X$$ is the Markov chain $$Y$$ defined as follows:


 * partition the interval $$[0,T]$$ into $$N$$ equal subintervals of width $$\Delta t>0$$: $$0 = \tau_0 < \tau_1 < \dots < \tau_N = T\text{ with }\tau_n:=n\Delta t\text{ and }\Delta t = \frac{T}{N}$$
 * set $$Y_0 = x_0;$$
 * recursively define $$Y_n$$ for $$1 \leq n \leq N$$ by: $$Y_{n + 1} = Y_n + a(Y_n) \Delta t + b(Y_n) \Delta W_n + \frac{1}{2} b(Y_n) b'(Y_n) \left( (\Delta W_n)^2 - \Delta t \right)$$ where $$b'$$ denotes the derivative of $$b(x)$$ with respect to $$x$$ and: $$\Delta W_n = W_{\tau_{n + 1}} - W_{\tau_n}$$ are independent and identically distributed normal random variables with expected value zero and variance $$\Delta t$$. Then $$Y_n$$ will approximate $$X_{\tau_n}$$ for $$0 \leq n \leq N$$, and increasing $$N$$ will yield a better approximation.

Note that when $$ b'(Y_n) = 0 $$, i.e. the diffusion term does not depend on $$X_{t}$$, this method is equivalent to the Euler–Maruyama method.

The Milstein scheme has both weak and strong order of convergence, $$\Delta t$$, which is superior to the Euler–Maruyama method, which in turn has the same weak order of convergence, $$\Delta t$$, but inferior strong order of convergence, $$\sqrt{\Delta t}$$.

Intuitive derivation
For this derivation, we will only look at geometric Brownian motion (GBM), the stochastic differential equation of which is given by: $$\mathrm{d} X_t = \mu X \mathrm{d} t + \sigma X d W_t$$ with real constants $$\mu$$ and $$\sigma$$. Using Itō's lemma we get: $$\mathrm{d}\ln X_t= \left(\mu - \frac{1}{2} \sigma^2\right)\mathrm{d}t+\sigma\mathrm{d}W_t$$

Thus, the solution to the GBM SDE is: $$ \begin{align} X_{t+\Delta t}&=X_t\exp\left\{\int_t^{t+\Delta t}\left(\mu-\frac{1}{2}\sigma^2\right)\mathrm{d}t+\int_t^{t+\Delta t}\sigma\mathrm{d}W_u\right\} \\ &\approx X_t\left(1+\mu\Delta t-\frac{1}{2} \sigma^2\Delta t+\sigma\Delta W_t+\frac{1}{2}\sigma^2(\Delta W_t)^2\right) \\ &= X_t + a(X_t)\Delta t+b(X_t)\Delta W_t+\frac{1}{2}b(X_t)b'(X_t)((\Delta W_t)^2-\Delta t) \end{align}$$ where $$ a(x) = \mu x, ~b(x) = \sigma x $$

See numerical solution is presented above for three different trajectories.

Computer implementation
The following Python code implements the Milstein method and uses it to solve the SDE describing the Geometric Brownian Motion defined by $$\begin{cases} dY_t= \mu Y \, {\mathrm d}t + \sigma Y \, {\mathrm d}W_t \\ Y_0=Y_\text{init} \end{cases}$$