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. [1] [2]
Consider the autonomous Itō stochastic differential equation: with initial condition , where stands for 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:
Note that when , i.e. the diffusion term does not depend on , 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, . [3]
For this derivation, we will only look at geometric Brownian motion (GBM), 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
See numerical solution is presented above for three different trajectories. [4]
The following Python code implements the Milstein method and uses it to solve the SDE describing the Geometric Brownian Motion defined by
# -*- coding: utf-8 -*-
# Milstein Method
import numpy as np
import matplotlib.pyplot as plt
class Model:
"""Stochastic model constants."""
μ = 3
σ = 1
def dW(Δt):
"""Random sample normal distribution."""
return np.random.normal(loc=0.0, scale=np.sqrt(Δt))
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)
ys0 = Y_INIT
for i in range(1, TS.size):
t = (i - 1) * DT
y = ysi - 1
dw = dW(DT)
# Sum up terms as in the Milstein method
ysi = y + \
Model.μ * y * DT + \
Model.σ * y * dw + \
(Model.σ**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)
{{
cite book}}
: CS1 maint: multiple names: authors list (
link)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. [1] [2]
Consider the autonomous Itō stochastic differential equation: with initial condition , where stands for 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:
Note that when , i.e. the diffusion term does not depend on , 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, . [3]
For this derivation, we will only look at geometric Brownian motion (GBM), 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
See numerical solution is presented above for three different trajectories. [4]
The following Python code implements the Milstein method and uses it to solve the SDE describing the Geometric Brownian Motion defined by
# -*- coding: utf-8 -*-
# Milstein Method
import numpy as np
import matplotlib.pyplot as plt
class Model:
"""Stochastic model constants."""
μ = 3
σ = 1
def dW(Δt):
"""Random sample normal distribution."""
return np.random.normal(loc=0.0, scale=np.sqrt(Δt))
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)
ys0 = Y_INIT
for i in range(1, TS.size):
t = (i - 1) * DT
y = ysi - 1
dw = dW(DT)
# Sum up terms as in the Milstein method
ysi = y + \
Model.μ * y * DT + \
Model.σ * y * dw + \
(Model.σ**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)
{{
cite book}}
: CS1 maint: multiple names: authors list (
link)