Choosing Euler-Maruyama or Milstein on Test SDEs

What the Choice Depends On

Euler–Maruyama and Milstein use the same Brownian increments, but they do not use the same information about the diffusion coefficient. Euler–Maruyama evaluates the drift and diffusion at the current state. Scalar Milstein also uses the derivative of the diffusion coefficient, adding a term that accounts for variation of the diffusion across a time step. That extra calculation can improve strong accuracy, but it is unnecessary when the diffusion is constant and may be inconvenient when its derivative is unavailable or costly to evaluate.

Let $(\Omega,\mathcal{F},(\mathcal{F}_t)_{t\ge0},\mathbb{P})$ carry a Brownian motion $B_t$. Consider the scalar Ito SDE $dX_t=\mu(X_t)\,dt+\sigma(X_t)\,dB_t$, with initial value $X_0=x_0$. The drift is $\mu$, and the diffusion coefficient is $\sigma$. The stochastic integral is the Ito integral: on each finite interval $[0,T]$, its integrand $\sigma(X_t)$ is assumed progressively measurable and in $\mathcal{L}^2$, meaning $\mathbb{E}\int_0^T\sigma(X_t)^2\,dt<\infty$. We assume coefficients for which the SDE has a unique strong solution; in the examples below, global Lipschitz and linear growth conditions suffice. A strong solution is built with the given Brownian motion on the given probability space. A weak solution only requires a probability space and Brownian motion on which the equation holds. Pathwise uniqueness means two solutions driven by the same Brownian motion and initial value agree almost surely; uniqueness in law means their distributions agree.

Definition: On a uniform grid $t_n=nh$, with increments $\Delta B_n=B_{t_{n+1}}-B_{t_n}$, Euler–Maruyama and scalar Milstein use the updates $$ X_{n+1}^{\mathrm{EM}}=X_n^{\mathrm{EM}}+\mu(X_n^{\mathrm{EM}})h+\sigma(X_n^{\mathrm{EM}})\Delta B_n, \qquad X_{n+1}^{\mathrm{Mil}}=X_n^{\mathrm{Mil}}+\mu(X_n^{\mathrm{Mil}})h+\sigma(X_n^{\mathrm{Mil}})\Delta B_n+\frac12\sigma(X_n^{\mathrm{Mil}})\sigma'(X_n^{\mathrm{Mil}})\big((\Delta B_n)^2-h\big). $$ The Milstein update assumes that $\sigma$ is differentiable where it is evaluated. Both methods use the same increments when compared strongly. The identity behind the correction is the Ito multiplication rule $dB_t\,dB_t=dt$, together with $dB_t\,dt=0$ and $dt\,dt=0$; these are not ordinary differential multiplication rules.

For a fixed time interval, strong error compares approximations driven by the same Brownian path, for example through $\big(\mathbb{E}|X_T-X_N|^2\big)^{1/2}$. Weak error instead compares expectations of functions of the endpoint. The earlier material on strong and weak convergence establishes the general distinction. Here, the comparison is about the strong error and the cost of producing each update.

First Test: Additive Noise

Worked Example: When the Correction Is Zero

Consider $dX_t=2\,dt+3\,dB_t$, with $X_0=1$. Use $h=\tfrac14$ and increments $0.1,-0.2,0.05,0.3$ to compute the four steps of both methods.

Here $\mu(x)=2$, $\sigma(x)=3$, and $\sigma'(x)=0$. Thus the Milstein correction is zero at every step. The common update is $X_{n+1}=X_n+2(\tfrac14)+3\Delta B_n=X_n+\tfrac12+3\Delta B_n$. The successive calculations are $X_1=1+\tfrac12+3(0.1)=1.8$, $X_2=1.8+\tfrac12+3(-0.2)=1.7$, $X_3=1.7+\tfrac12+3(0.05)=2.35$, and $X_4=2.35+\tfrac12+3(0.3)=3.75$. Both schemes give $3.75$.

Check: Integrating the constant-coefficient SDE gives $X_1=1+2(1)+3B_1$. The supplied increments sum to $0.1-0.2+0.05+0.3=0.25$, so this expression gives $1+2+3(0.25)=3.75$, agreeing with the step-by-step updates. In this case, there is no discretization error at the grid points.

This example illustrates an exact structural criterion: if $\sigma$ is constant, then $\sigma'=0$ and the scalar Milstein correction vanishes. It does not say that Euler–Maruyama is exact for every SDE with constant diffusion; a nonconstant drift can still require time discretization.

Second Test: State-Dependent Diffusion

Worked Example: One Geometric Brownian Motion Step

For $dS_t=0.1S_t\,dt+0.4S_t\,dB_t$, take $S_0=2$, $h=\tfrac14$, and $\Delta B=0.3$. Compute both updates and compare them with the exact solution over this step.

The coefficients are $\mu(s)=0.1s$, $\sigma(s)=0.4s$, and $\sigma'(s)=0.4$. Euler–Maruyama gives $S_1^{\mathrm{EM}}=2+0.1(2)(0.25)+0.4(2)(0.3)=2+0.05+0.24=2.29$.

For Milstein, the correction coefficient is $\tfrac12\sigma(2)\sigma'(2)=\tfrac12(0.8)(0.4)=0.16$, while $(\Delta B)^2-h=(0.3)^2-0.25=0.09-0.25=-0.16$. Thus the correction is $0.16(-0.16)=-0.0256$, and $S_1^{\mathrm{Mil}}=2.29-0.0256=2.2644$.

The exact solution is $S_t=2\exp((0.1-\tfrac12(0.4)^2)t+0.4B_t)$. To verify it, set $f(t,b)=2\exp((0.1-\tfrac12(0.4)^2)t+0.4b)$. Its derivatives are $f_t=(0.1-\tfrac12(0.4)^2)f=0.02f$, $f_b=0.4f$, and $f_{bb}=0.16f$. Ito’s formula, including its correction, gives $df=(f_t+\tfrac12f_{bb})dt+f_b\,dB_t=(0.02f+0.08f)dt+0.4f\,dB_t$, which is the stated SDE. At the supplied increment, the exact value is $2e^{0.02(0.25)+0.4(0.3)}=2e^{0.125}$.

Check: Substitution into the exact formula gives $2e^{0.125}\approx2.2663$. The absolute errors for this particular increment are $|2.29-2e^{0.125}|\approx0.0237$ for Euler–Maruyama and $|2.2644-2e^{0.125}|\approx0.0019$ for Milstein. This single increment is an illustration, not a convergence-rate estimate: strong rates concern expected errors over coupled paths as the step size decreases.

On this model, the correction requires one extra multiplication and a squared Brownian increment. For a general scalar diffusion it also requires evaluating $\sigma'$. In a vector SDE, Milstein corrections can involve additional noise interactions; the earlier material on multidimensional and noncommutative Milstein schemes addresses that complication. The scalar comparison here does not imply that this correction alone resolves multidimensional noise effects.

A Coupled Strong-Error Calculation

For geometric Brownian motion we can compare the endpoint mean-square errors exactly, rather than infer a rate from one sample path. Use $dS_t=S_t\,dB_t$, $S_0=1$, and $T=1$. This SDE has globally Lipschitz drift $\mu(s)=0$ and diffusion $\sigma(s)=s$, both of linear growth. Its exact solution is $S_1=e^{-1/2+B_1}$. For $N$ steps, let $h=1/N$ and use the same Brownian increments for the exact solution and for each numerical scheme.

Theorem: In this test, the endpoint mean-square errors are $$ \mathbb{E}|S_1-S_N^{\mathrm{EM}}|^2=e-(1+h)^N, \qquad \mathbb{E}|S_1-S_N^{\mathrm{Mil}}|^2=e-(1+h+\tfrac12h^2)^N. $$ Consequently, as $h=1/N$ tends to zero, the root mean-square endpoint errors have orders $h^{1/2}$ for Euler–Maruyama and $h$ for Milstein.

Proof. Write $Z=\Delta B_n$, which has distribution $N(0,h)$, and define the one-step exact multiplier $G=e^{-h/2+Z}$. The exact endpoint is the product of the $N$ multipliers $G$. The numerical one-step multipliers are $A_{\mathrm{EM}}=1+Z$ and $A_{\mathrm{Mil}}=1+Z+\tfrac12(Z^2-h)$. Independent increments make the corresponding multipliers independent across steps. The Gaussian moment generating function gives $\mathbb{E}G^2=e^h$, so $\mathbb{E}S_1^2=e^{Nh}=e$.

For Euler–Maruyama, $\mathbb{E}A_{\mathrm{EM}}^2=1+\mathbb{E}Z^2=1+h$. Also, completing the square in the normal density shows that weighting by $e^{-h/2+Z}$ shifts the mean of $Z$ from $0$ to $h$, while the resulting density still integrates to one. Therefore $\mathbb{E}[A_{\mathrm{EM}}G]=1+\mathbb{E}_{\mathrm{shifted}}Z=1+h$. It follows by independence that $\mathbb{E}(S_N^{\mathrm{EM}})^2=(1+h)^N$ and $\mathbb{E}[S_1S_N^{\mathrm{EM}}]=(1+h)^N$.

For Milstein, put $U=Z+\tfrac12(Z^2-h)$, so $A_{\mathrm{Mil}}=1+U$. The centered normal moments give $\mathbb{E}U=0$, $\mathbb{E}Z^2=h$, $\mathbb{E}(Z^2-h)^2=2h^2$, and $\mathbb{E}[Z(Z^2-h)]=0$. Hence $\mathbb{E}U^2=h+\tfrac12h^2$, and $\mathbb{E}A_{\mathrm{Mil}}^2=1+h+\tfrac12h^2$. Under the shifted normal law used above, $Z$ has mean $h$ and variance $h$, so $\mathbb{E}_{\mathrm{shifted}}(Z^2-h)=h^2$. Thus $\mathbb{E}[A_{\mathrm{Mil}}G]=1+h+\tfrac12h^2$. Independence now gives $\mathbb{E}(S_N^{\mathrm{Mil}})^2=\mathbb{E}[S_1S_N^{\mathrm{Mil}}]=(1+h+\tfrac12h^2)^N$.

For either method, expand the squared error as $\mathbb{E}S_1^2+\mathbb{E}S_N^2-2\mathbb{E}[S_1S_N]$ and substitute the moments just calculated. This gives both stated error formulas. To obtain the rates, use the Taylor expansion $\log(1+h)=h-\tfrac12h^2+O(h^3)$, so $e-(1+h)^{1/h}=\tfrac12e\,h+O(h^2)$. Also $\log(1+h+\tfrac12h^2)=h+O(h^3)$, so $e-(1+h+\tfrac12h^2)^{1/h}=O(h^2)$; its leading term is positive and of order $h^2$, as follows by expanding one further term, $\log(1+h+\tfrac12h^2)=h-\tfrac16h^3+O(h^4)$. Taking square roots yields orders $h^{1/2}$ and $h$, respectively. $\square$

This calculation uses the same increments for the numerical and exact paths, as a strong comparison requires. For example, at $N=1$, the formulas give mean-square errors $e-2$ for Euler–Maruyama and $e-\tfrac52$ for Milstein. Both are positive, and the Milstein error is smaller in this test. The asymptotic rate statement concerns refinement $N\to\infty$, not a guarantee that Milstein is closer on every individual increment or at every coarse step.

How to Choose and What Can Go Wrong

The standard scalar convergence results established in the earlier material on strong order for Euler–Maruyama and scalar Milstein give strong orders $1/2$ and $1$, respectively, under their stated regularity and moment hypotheses. The exact geometric Brownian calculation above demonstrates those rates in one test case. It does not replace those hypotheses: lack of smoothness, growth or moment control can invalidate an order claim. Nor does a higher strong order automatically mean a faster computation. If evaluating $\sigma'$ is expensive, derivatives are unavailable, or the simulation is dominated by another cost, the extra accuracy per step may not repay its cost.

Common mistake: Omitting the $-h$ in the Milstein correction changes the method. For the geometric Brownian step above, using $\tfrac12\sigma(S_0)\sigma'(S_0)(\Delta B)^2$ instead of $\tfrac12\sigma(S_0)\sigma'(S_0)((\Delta B)^2-h)$ gives a false correction $0.16(0.09)=0.0144$, and the incorrect update $2.3044$. The proper correction is $0.16(0.09-0.25)=-0.0256$, giving $2.2644$. The subtraction of $h$ is essential: it is the centered Brownian-increment-square term generated by the Ito expansion.

Test or consideration	Euler–Maruyama	Scalar Milstein
Constant diffusion	Correction-free update	Same update; correction vanishes
State-dependent scalar diffusion	Uses $\mu$, $\sigma$, and $\Delta B$	Also evaluates $\sigma'$ and $(\Delta B)^2-h$
Strong endpoint rate in the geometric Brownian test	$h^{1/2}$	$h$
Practical trade-off	Simpler update	Higher accuracy per step under suitable assumptions

A useful decision is therefore model-specific: check whether the diffusion depends on the state, whether its derivative can be evaluated reliably, and whether the target is strong path accuracy. If the diffusion is additive, the scalar correction provides no benefit. If it varies smoothly with the state and strong accuracy matters, Milstein is a natural candidate. Always compare methods under a shared Brownian coupling and assess errors over a sequence of step sizes; a favorable single path is not evidence of a convergence rate.

Check Your Understanding

Use the update formulas and distinguish a structural comparison from a conclusion based on one sample path.

For $dX_t=dt+2\,dB_t$, explain why the scalar Milstein correction is zero.
For $\sigma(x)=x^2$, write the diffusion-derivative factor $\tfrac12\sigma(x)\sigma'(x)$ in the Milstein correction.
Why must Euler–Maruyama and Milstein use the same Brownian increments for a direct strong-error comparison?
In the geometric Brownian test, what are the endpoint mean-square error formulas when $N=1$?
Why does a smaller error on one increment not establish that a method has a higher strong convergence order?

Choosing Euler-Maruyama or Milstein on Test SDEs

What You'll Learn