Milstein method

Remarks

The scheme is the Itô–Taylor expansion of the solution over one step, cut after the terms of order $h$ in mean square. Applying Itô's lemma to $b(X_t)$ shows that the leading neglected term of the Euler–Maruyama method is $b b' \int_{t_n}^{t_{n+1}} (W_s - W_{t_n})\, dW_s$, and this double integral equals $\frac12((\Delta W_n)^2 - h)$. The Euler–Maruyama method omits this term and has strong order only $1/2$ when $b' \ne 0$; when $b$ is constant, both methods coincide. Example: for geometric Brownian motion $dX = \mu X\, dt + \sigma X\, dW$, the method reads $Y_{n+1} = Y_n \left( 1 + \mu h + \sigma \Delta W_n + \frac{\sigma^2}{2} ((\Delta W_n)^2 - h) \right)$. For systems driven by several Wiener processes, the analogous scheme needs the iterated integrals $\int\!\!\int dW^i\, dW^j$, which are not functions of the increments unless a commutativity condition holds.

Going deeper

For a student

Strong order measures the error of individual paths and matters for multilevel Monte Carlo and for pathwise simulation. Weak order measures the error of expectations, and there Euler–Maruyama and Milstein both have order $1$. The extra term costs one derivative of $b$; derivative-free variants replace $b'$ by a difference quotient, as in stochastic Runge–Kutta methods.