Milstein method

Consider the autonomous Itō stochastic differential equation:

$$\backslash mathrm\{d\}\; X\_t\; =\; a(X\_t)\; \backslash ,\; \backslash mathrm\{d\}\; t\; +\; b(X\_t)\; \backslash ,\; \backslash mathrm\{d\}\; W\_t$$

with initial condition $X\_\{0\}\; =\; x\_\{0\}$, where $W\_\{t\}$ denotes 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 $\backslash Delta\; t>0$:
• Set $Y\_0\; =\; x\_0;$
• Recursively define $Y\_n$ for $1\; \backslash leq\; n\; \backslash leq\; N$ by: $$Y\_\{n\; +\; 1\}\; =\; Y\_n\; +\; a(Y\_n)\; \backslash Delta\; t\; +\; b(Y\_n)\; \backslash Delta\; W\_n\; +\; \backslash frac\{1\}\{2\}\; b(Y\_n)\; b\text{'}(Y\_n)\; \backslash left(\; (\backslash Delta\; W\_n)^2\; -\; \backslash Delta\; t\; \backslash right)$$ where $b\text{'}$ denotes the derivative of $b(x)$ with respect to $x$ and: $$\backslash Delta\; W\_n\; =\; W\_\{\backslash tau\_\{n\; +\; 1\}\}\; -\; W\_\{\backslash tau\_n\}$$ are independent and identically distributed normal random variables with expected value zero and variance $\backslash Delta\; t$. Then $Y\_n$ will approximate $X\_\{\backslash tau\_n\}$ for $0\; \backslash leq\; n\; \backslash leq\; N$, and increasing $N$ will yield a better approximation.

Note that when $b\text{'}(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 $\backslash Delta\; t$ which is superior to the Euler–Maruyama method, which in turn has the same weak order of convergence $\backslash Delta\; t$ but inferior strong order of convergence $\backslash sqrt\{\backslash Delta\; t\}$.V. Mackevičius, Introduction to Stochastic Analysis, Wiley 2011

For this derivation, we will only look at geometric Brownian motion (GBM), the stochastic differential equation of which is given by:

$$\backslash mathrm\{d\}\; X\_t\; =\; \backslash mu\; X\; \backslash mathrm\{d\}\; t\; +\; \backslash sigma\; X\; d\; W\_t$$

with real constants $\backslash mu$ and $\backslash sigma$. Using Itō's lemma we get:

$$\backslash mathrm\{d\}\backslash ln\; X\_t=\; \backslash left(\backslash mu\; -\; \backslash frac\{1\}\{2\}\; \backslash sigma^2\backslash right)\backslash mathrm\{d\}t+\backslash sigma\backslash 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)\; =\; \backslash mu\; x,\; ~b(x)\; =\; \backslash 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/

File:SDE 2.jpg

The following Python code implements the Milstein method and uses it to solve the SDE describing geometric Brownian motion defined by

$$\backslash begin\{cases\}$$

dY_t= \mu Y \, {\mathrm d}t + \sigma Y \, {\mathrm d}W_t \\

Y_0=Y_\text{init}

\end{cases}

1. -*- coding: utf-8 -*-
2. 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[0] = Y_INIT

for i in range(1, TS.size):

t = (i - 1) * DT

y = ys[i - 1]

dw = dW(DT)

# Sum up terms as in the Milstein method

ys[i] = 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)

• Euler–Maruyama method

{{Reflist}}

• {{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 }}

Category:Numerical differential equations

Category:Stochastic differential equations