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 Milstein who first published it in 1974.

Description

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

* Partition the interval ,T/math> into $N$ equal subintervals of width $\Delta t>0$: $0 = \tau_0 < \tau_1 < \dots < \tau_N = T\text\tau_n:=n\Delta t\text\Delta t = \frac$
* Set $Y_0 = x_0;$
* Recursively define $Y_n$ for $1 \leq n \leq N$ by: $Y_ = Y_n + a(Y_n) \Delta t + b(Y_n) \Delta W_n + \frac 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_ - W_$ are independent and identically distributed normal random variables with expected value zero and variance $\Delta t$. Then $Y_n$ will approximate $X_$ 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_$) 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$.

Intuitive derivation

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

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

The numerical solution is presented in the graphic for three different trajectories.Umberto Picchini, SDE Toolbox: simulation and estimation of stochastic differential equations with Matlab. http://sdetoolbox.sourceforge.net/

Computer implementation

The following Python code implements the Milstein method and uses it to solve the SDE describing geometric Brownian motion defined by
$\begin dY_t= \mu Y \, t + \sigma Y \, W_t \\ Y_0=Y_\text \end$

# -*- coding: utf-8 -*-
# Milstein Method

import numpy as np
import matplotlib.pyplot as plt

class Model:
"""Stochastic model constants."""
mu = 3
sigma = 1

def dW(dt):
"""Random sample normal distribution."""
return np.random.normal(loc=0.0, scale=np.sqrt(dt))

def run_simulation():
""" Return the result of one full simulation."""
# One second and thousand grid points
T_INIT = 0
T_END = 1
N = 1000 # Compute 1000 grid points
DT = float(T_END - T_INIT) / N
TS = np.arange(T_INIT, T_END + DT, DT)

Y_INIT = 1

# Vectors to fill
ys = np.zeros(N + 1)
ys = Y_INIT
for i in range(1, TS.size):
t = (i - 1) * DT
y = ys - 1 dw = dW(DT)

# Sum up terms as in the Milstein method
ys = y + \
Model.mu * y * DT + \
Model.sigma * y * dw + \
(Model.sigma**2 / 2) * y * (dw**2 - DT)

return TS, ys

def plot_simulations(num_sims: int):
"""Plot several simulations in one image."""
for _ in range(num_sims):
plt.plot(*run_simulation())

plt.xlabel("time (s)")
plt.ylabel("y")
plt.grid()
plt.show()

if __name__

"__main__":
NUM_SIMS = 2
plot_simulations(NUM_SIMS)

See also

* Euler–Maruyama method

References

Further reading

* {{cite book , author=Kloeden, P.E., & Platen, E. , title=Numerical Solution of Stochastic Differential Equations , publisher=Springer, Berlin , year=1999 , isbn=3-540-54062-8

Numerical differential equations
Stochastic differential equations