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 the method in 1974.
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:
- partition the interval into equal subintervals of width :
- recursively define for by:
where denotes the derivative of with respect to and:
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, .
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:
See numerical solution is presented above for three different trajectories.
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 num_sims = 1 # One Example # 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 ## Initial Conditions y_init = 1 mu = 3 sigma = 1 # dw Random process def dW(delta_t): """" Random sample normal distribution""" return np.random.normal(loc=0.0, scale=np.sqrt(delta_t)) # vectors to fill ts = np.arange(t_init, t_end + dt, dt) ys = np.zeros(N + 1) ys = y_init # Loop for _ in range(num_sims): for i in range(1, ts.size): t = (i - 1) * dt y = ys[i - 1] # Milstein method ys[i] = y + mu * dt * y + sigma* y* dW(dt) + 0.5* sigma**2 * (dW(dt)**2 - dt) plt.plot(ts, ys) # Plot plt.xlabel("time (s)") plt.grid() h = plt.ylabel("y") h.set_rotation(0) plt.show()
- Mil'shtein, G. N. (1974). "Approximate integration of stochastic differential equations". Teoriya Veroyatnostei i ee Primeneniya (in Russian). 19 (3): 583–588.
- Mil’shtein, G. N. (1975). "Approximate Integration of Stochastic Differential Equations". Theory of Probability & Its Applications. 19 (3): 557–000. doi:10.1137/1119062.
- V. Mackevičius, Introduction to Stochastic Analysis, Wiley 2011
- Umberto Picchini, SDE Toolbox: simulation and estimation of stochastic differential equations with Matlab. http://sdetoolbox.sourceforge.net/