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:
with initial condition, where denotes the Wiener process, and suppose that we wish to solve this SDE on some interval of time . Then the Milstein approximation to the true solution is the Markov chain defined as follows:

Partition the interval into equal subintervals of width :
Set
Recursively define for by: where denotes the derivative of with respect to and: are independent and identically distributed normal random variables with expected value zero and variance. Then will approximate for, and increasing will yield a better approximation.

Note that when this method is equivalent to the Euler–Maruyama method.
The Milstein scheme has both weak and strong order of convergence which is superior to the Euler–Maruyama method, which in turn has the same weak order of convergence but inferior strong order of convergence.

Intuitive derivation

For this derivation, we will only look at geometric Brownian motion, the stochastic differential equation of which is given by:
with real constants and. Using Itō's lemma we get:
Thus, the solution to the GBM SDE is:
where
The numerical solution is presented in the graphic for three different trajectories.

Computer implementation

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

-*- 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:
"""Random sample normal distribution."""
return np.random.normal
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 / N
TS = np.arange
Y_INIT = 1
# Vectors to fill
ys = np.zeros
ys = Y_INIT
for i in range:
t = * DT
y = ys
dw = dW
# Sum up terms as in the Milstein method
ys = y + \
Model.mu * y * DT + \
Model.sigma * y * dw + \
* y *
return TS, ys
def plot_simulations:
"""Plot several simulations in one image."""
for _ in range:
plt.plot)
plt.xlabel
plt.ylabel
plt.grid
plt.show
if __name__ "__main__":
NUM_SIMS = 2
plot_simulations