the Milstein scheme
/ MIL-shtine /
The Milstein scheme is the next step up from Euler-Maruyama: a numerical method for SDEs that adds one correction term to recover full strong order 1, matching the deterministic Euler method's pathwise accuracy. It answers the question Euler-Maruyama leaves open: where did the lost half-order of pathwise accuracy go, and how do I get it back? The answer is that Euler ignores how the diffusion coefficient sigma changes within a time step, and Milstein puts exactly that contribution back in via Ito's formula.
For the scalar SDE dX = b(X) dt + sigma(X) dB, the Milstein recursion is X_{n+1} = X_n + b(X_n) h + sigma(X_n) Delta B_n + (1/2) sigma(X_n) sigma'(X_n) ( (Delta B_n)^2 - h ), where sigma' is the derivative of the diffusion coefficient and Delta B_n ~ N(0, h). The new term comes from Ito-expanding sigma(X_s) over the step: sigma(X_s) is approximately sigma(X_n) + sigma'(X_n)(X_s - X_n), and integrating the leading stochastic part of X_s - X_n against dB produces the iterated Ito integral integral integral dB dB = ((Delta B)^2 - h)/2 — and the appearance of the -h is the (dB)^2 = dt rule again. Because this captures the within-step variation of sigma, the strong order improves to 1: E[ |X_T - X^h_T| ] = O(h). The weak order remains 1, the same as Euler-Maruyama, so the benefit of Milstein is entirely in pathwise accuracy; if you only need expectations, Euler is already weak order 1 and Milstein gives no asymptotic improvement there.
Milstein is the method of choice when you need accurate sample paths — pathwise functionals, path-dependent options, sensitivity along trajectories. Two honest caveats stand out. First, Milstein needs the derivative sigma' of the diffusion coefficient, so it requires more smoothness than Euler and is awkward when sigma is only Lipschitz or given as a black box (then Euler, or a derivative-free Runge-Kutta-type variant, is used). Second, and importantly, the simple closed form above is special to one dimension and to commutative noise. In the general multidimensional case the correction involves iterated stochastic integrals integral integral dB^i dB^j (Levy areas) that are NOT simple functions of the increments and are expensive to simulate unless a commutativity condition on the diffusion columns holds — this is the real obstruction to high strong order in multiple dimensions.
For geometric Brownian motion sigma(x) = sigma x so sigma'(x) = sigma, and the Milstein step is X_{n+1} = X_n + mu X_n h + sigma X_n Delta B_n + (1/2) sigma^2 X_n ((Delta B_n)^2 - h). The extra term raises the strong order from 1/2 (Euler-Maruyama) to 1, halving pathwise error when you halve h, while leaving the weak order at 1.
Milstein for GBM: the (Delta B)^2 - h correction lifts strong order from 1/2 to 1.
Milstein improves only the STRONG (pathwise) order to 1; the weak order stays 1, so it gives no advantage for plain expectations over Euler-Maruyama. The clean formula is one-dimensional; in higher dimensions it needs iterated integrals (Levy areas) that are hard to simulate unless the noise commutes, and it requires the derivative sigma'.