Proving Strong Order One for Scalar Milstein

From a Local Correction to a Global Error Bound

The scalar Milstein update adds a centered Brownian-increment-square term to Euler–Maruyama. Its purpose is not merely to improve a numerical path empirically: under smoothness and growth assumptions, it cancels the leading local error that prevents Euler–Maruyama from having strong order one. The proof has two parts. First, expand the exact solution over a single time step and bound the remainder. Then show that the centered random remainders accumulate like a martingale, while their conditional means accumulate more slowly.

Work on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_t)_{t\ge0},\mathbb{P})$ carrying a standard Brownian motion $B_t$. Consider the scalar SDE $dX_t=a(X_t)\,dt+b(X_t)\,dB_t$, with deterministic initial value $X_0=x_0$. Every stochastic integral below is an Ito integral with a progressively measurable integrand in $\mathcal{L}^2$, meaning $\mathbb{E}\int_0^T H_s^2\,ds<\infty$. Assume $a,b\in C^2(\mathbb{R})$, and that $a,b$ and their first and second derivatives are bounded. These assumptions imply global Lipschitz and linear growth bounds; they also ensure a unique strong solution and the moment bounds used below.

For a grid $t_k=kh$, $k=0,\ldots,N$, where $h=T/N$, write $\Delta B_k=B_{t_{k+1}}-B_{t_k}$. The scalar Milstein approximation is defined recursively by the following update.

Definition: With $Y_0=x_0$, the scalar Milstein scheme is $$ Y_{k+1}=Y_k+a(Y_k)h+b(Y_k)\Delta B_k+\frac12 b(Y_k)b'(Y_k)\big((\Delta B_k)^2-h\big). $$ Its strong order-one error bound means that, for fixed $T$, there is a constant $C$, independent of $N$, such that $\big(\mathbb{E}\max_{0\le k\le N}|X_{t_k}-Y_k|^2\big)^{1/2}\le Ch$.

The word “strong” refers to comparing the exact and approximate paths using the same Brownian motion. The bound here is for the maximum over grid points in mean square; it does not assert that every individual path has a deterministic error bound.

The One-Step Expansion

Set $c(x)=\tfrac12 b(x)b'(x)$. Starting an exact solution at time $t$ from $x=X_t$, its integral equation gives $X_{t+h}=x+\int_t^{t+h}a(X_s)\,ds+\int_t^{t+h}b(X_s)\,dB_s$. The key term comes from replacing $b(X_s)$ by its first-order expansion around $x$. In that expansion, the first Brownian displacement is $b(x)(B_s-B_t)$. Its integral is evaluated using Ito’s formula for $f(z)=z^2$: $f'(z)=2z$, $f''(z)=2$, and the Ito correction is $\tfrac12f''(z)\,ds=ds$. Thus $d((B_s-B_t)^2)=2(B_s-B_t)\,dB_s+ds$, and integration yields $\int_t^{t+h}(B_s-B_t)\,dB_s=\tfrac12((\Delta B)^2-h)$.

Lemma: Under the stated assumptions, for every grid interval, $$ X_{t_{k+1}}=X_{t_k}+a(X_{t_k})h+b(X_{t_k})\Delta B_k +c(X_{t_k})\big((\Delta B_k)^2-h\big)+\rho_k, $$ where $\mathbb{E}[\rho_k\mid\mathcal{F}_{t_k}]=q_k$ and $\mathbb{E}[\rho_k^2\mid\mathcal{F}_{t_k}]\le Ch^3$, $|q_k|\le Ch^2$. The constant may depend on $T$ and the coefficient bounds, but not on $k$ or $h$.

Proof. Fix $t=t_k$, let $x=X_t$, and write $Y_s^*=X_s-x$ for $t\le s\le t+h$. Taylor’s theorem gives $b(X_s)=b(x)+b'(x)Y_s^*+R_s$, with $|R_s|\le \tfrac12\|b''\|_\infty|Y_s^*|^2$. Also write $Y_s^*=b(x)(B_s-B_t)+A_s$, where $A_s=\int_t^s a(X_r)\,dr+\int_t^s(b(X_r)-b(x))\,dB_r$. Substituting both identities into the integral equation gives the displayed expansion, with $$ \rho_k=\int_t^{t+h}(a(X_s)-a(x))\,ds +b'(x)\int_t^{t+h}A_s\,dB_s +\int_t^{t+h}R_s\,dB_s. $$ The leading stochastic term is $b'(x)b(x)\int_t^{t+h}(B_s-B_t)\,dB_s =c(x)((\Delta B)^2-h)$, by the Ito calculation above.

Here are the estimates for the remainder. Bounded coefficients, the Ito isometry, and the standard Burkholder–Davis–Gundy moment inequality give $\mathbb{E}[|Y_s^*|^2\mid\mathcal{F}_t]\le C(s-t)$ and $\mathbb{E}[|Y_s^*|^4\mid\mathcal{F}_t]\le C(s-t)^2$. The Burkholder–Davis–Gundy inequality is a standard martingale inequality; its proof is beyond this tutorial. Since $a$ and $b$ are Lipschitz, the squared $L^2$ norm of the first, drift-integral term in $\rho_k$ is at most $h\int_t^{t+h}C(s-t)\,ds\le Ch^3$. The Ito isometry and $\mathbb{E}[|A_s|^2\mid\mathcal{F}_t]\le C(s-t)^2$ bound the squared $L^2$ norm of the second term by $C\int_t^{t+h}(s-t)^2\,ds\le Ch^3$. For the third term, the Taylor bound and the fourth-moment estimate give a squared $L^2$ norm at most $C\int_t^{t+h}(s-t)^2\,ds\le Ch^3$. The inequality $|u+v+w|^2\le3(|u|^2+|v|^2+|w|^2)$ therefore gives the stated conditional second-moment bound.

Both stochastic integrals in $\rho_k$ have conditional mean zero. For the remaining drift term, Taylor-expand $a(x+Y_s^*)$ to first order; its remainder is bounded by $\tfrac12\|a''\|_\infty|Y_s^*|^2$. Moreover, $\mathbb{E}[Y_s^*\mid\mathcal{F}_t]=\mathbb{E}[\int_t^s a(X_r)\,dr\mid\mathcal{F}_t]$, since the Brownian integral has conditional mean zero. Boundedness of $a$ makes this conditional expectation at most $C(s-t)$ in absolute value. The Taylor expansion and the second-moment estimate now give $|\mathbb{E}[a(X_s)-a(x)\mid\mathcal{F}_t]|\le C(s-t)$. Integrating over $s\in[t,t+h]$ proves $|q_k|\le Ch^2$. $\square$

Why the Local Remainder Still Gives Order One

A bound of order $h^{3/2}$ on each local remainder would not suffice if all remainders were added using the triangle inequality: $N$ such bounds would only give order $h^{1/2}$. The improvement comes from centering. Write $\rho_k=q_k+m_k$, where $q_k=\mathbb{E}[\rho_k\mid\mathcal{F}_{t_k}]$ and $m_k$ has conditional mean zero. The $q_k$ terms have total size at most $NCh^2= CTh$. The $m_k$ terms are martingale differences: their accumulated mean-square size is controlled by the sum of their conditional variances, which is at most $NCh^3=CTh^2$.

Theorem: Under the assumptions above, there is a constant $C_T$, independent of $N$, such that the scalar Milstein scheme satisfies $$ \mathbb{E}\left[\max_{0\le k\le N}|X_{t_k}-Y_k|^2\right]\le C_T h^2. $$ In particular, its grid-point strong error is of order one in root mean square.

Proof. Put $e_k=X_{t_k}-Y_k$. Subtract the Milstein update from the one-step expansion in the lemma. With $a_k=a(X_{t_k})-a(Y_k)$, $b_k=b(X_{t_k})-b(Y_k)$, and $c_k=c(X_{t_k})-c(Y_k)$, this gives $$ e_{k+1}=e_k+h a_k+b_k\Delta B_k+c_k\big((\Delta B_k)^2-h\big)+q_k+m_k. $$ The functions $a,b,c$ are globally Lipschitz: this holds for $a,b$ by assumption, and for $c=\tfrac12bb'$ because $b,b',b''$ are bounded. Consequently, $|a_k|+|b_k|+|c_k|\le C|e_k|$.

Sum the recursion from $k=0$ to $n-1$. The drift sum satisfies $\max_{n\le N}|\sum_{k<n}ha_k|^2\le T\sum_{k<N}h|a_k|^2\le C\sum_{k<N}h|e_k|^2$, by Cauchy–Schwarz. The sum of the $q_k$ has absolute value at most $NCh^2=CTh$. The remaining terms form a martingale sum: $\Delta B_k$, $(\Delta B_k)^2-h$, and $m_k$ are conditionally centered given $\mathcal{F}_{t_k}$. The increments from distinct time steps are orthogonal in $L^2$, by conditional expectation.

For completeness, $\mathbb{E}[(\Delta B_k)^2\mid\mathcal{F}_{t_k}]=h$ and $\mathbb{E}[((\Delta B_k)^2-h)^2\mid\mathcal{F}_{t_k}]=2h^2$, using the second and fourth moments of a centered normal variable with variance $h$. The conditional cross moment between $\Delta B_k$ and $(\Delta B_k)^2-h$ is zero by symmetry. Using $|u+v+w|^2\le3(u^2+v^2+w^2)$, the martingale increment at step $k$ has second moment at most $C(h\mathbb{E}|e_k|^2+h^3)$. Doob’s $L^2$ maximal inequality, which bounds the expected squared maximum of a martingale by four times its terminal second moment, therefore bounds the expected squared maximum of the accumulated martingale by $C\sum_{k<N}(h\mathbb{E}|e_k|^2+h^3)$.

Combining these bounds, and using $|u_1+\cdots+u_4|^2\le4\sum_i|u_i|^2$, yields, for each $n\le N$, $$ u_n:=\mathbb{E}\max_{0\le j\le n}|e_j|^2 \le C\sum_{k<n}h\,u_k+Ch^2. $$ Here $e_0=0$, and $Nh^3=Th^2$. The discrete Gronwall inequality gives $u_N\le Ch^2$: iterating the displayed inequality bounds $u_N$ by $Ch^2(1+Ch)^N\le Ch^2e^{CT}$. Taking square roots proves the claimed order-one bound. $\square$

Worked Examples

Worked Example: Additive Noise Is Captured Exactly

Take $dX_t=2\,dt+3\,dB_t$ with $X_0=1$, on a grid with step $h=\tfrac14$. Find one Milstein step and compare it with the exact solution.

Here $a(x)=2$, $b(x)=3$, and $b'(x)=0$, so $c(x)=\tfrac12b(x)b'(x)=0$. The update from $Y_k=1+2t_k+3B_{t_k}$ is $Y_{k+1}=Y_k+2(\tfrac14)+3\Delta B_k=Y_k+\tfrac12+3\Delta B_k$. The exact integral equation gives $X_{t_{k+1}}=X_{t_k}+\int_{t_k}^{t_{k+1}}2\,ds+\int_{t_k}^{t_{k+1}}3\,dB_s =X_{t_k}+\tfrac12+3\Delta B_k$. Thus, if $Y_k=X_{t_k}$, then $Y_{k+1}=X_{t_{k+1}}$ exactly. Starting from $Y_0=X_0=1$, induction gives equality at every grid point.

Check: the Milstein correction vanishes because the diffusion coefficient is constant. Substitution into the exact integral equation gives precisely the same drift increment $\tfrac12$ and Brownian increment $3\Delta B_k$, so the computed update has zero error.

Worked Example: The Geometric Brownian One-Step Correction

For $dX_t=X_t\,dB_t$, start a step at $X_t=2$ with $h=\tfrac14$, and denote the Brownian increment by $z$. Compare the Euler and Milstein updates with the exact one-step solution.

The coefficients are $a(x)=0$, $b(x)=x$, and $c(x)=\tfrac12x$. Hence Euler gives $2+2z$, whereas Milstein gives $2+2z+1\cdot(z^2-\tfrac14)=\tfrac74+2z+z^2$. The exact solution over this step is $2\exp(z-\tfrac18)$, obtained by Ito’s formula applied to the exponential solution. Expanding the exponent through its quadratic Brownian terms, $\exp(z-\tfrac18)=1+z+\tfrac12(z^2-\tfrac14)+$ terms of higher stochastic order. Multiplying by $2$ gives the Milstein terms $2+2z+(z^2-\tfrac14)$; Euler omits the last parenthesized term.

At the algebraic value $z=0$, the exact endpoint is $2e^{-1/8}$, Euler gives $2$, and Milstein gives $\tfrac74$. The exact Milstein discrepancy at this value is $2e^{-1/8}-\tfrac74$. This particular increment is an algebraic check, not an event of positive probability under Brownian motion.

Check: substituting $z=0$ into the Milstein formula gives $2+0+(0-\tfrac14)=\tfrac74$, as above. The correction has mean zero over a Brownian step because $\mathbb{E}[z^2]=h=\tfrac14$, despite being nonzero at this selected increment.

Worked Example: Size of the Leading Term Euler Omits

At $x=2$ for the same diffusion $b(x)=x$, take a step $h=\tfrac14$. Compute the variance and standard deviation of the correction term in the Milstein update.

The coefficient is $\tfrac12xb'(x)=\tfrac12(2)(1)=1$. If $Z=\Delta B\sim N(0,\tfrac14)$, the correction is $Z^2-\tfrac14$. Its mean is $\mathbb{E}[Z^2]-\tfrac14=\tfrac14-\tfrac14=0$. A centered normal variable of variance $h$ has fourth moment $3h^2$, so $\mathbb{E}[(Z^2-\tfrac14)^2]=3(\tfrac14)^2-2(\tfrac14)(\tfrac14)+(\tfrac14)^2=\tfrac3{16}-\tfrac18+\tfrac1{16}=\tfrac18$. The variance is therefore $\tfrac18$, and the standard deviation is $\sqrt{\tfrac18}=\tfrac1{2\sqrt2}$.

Check: the same variance follows from $\operatorname{Var}(Z^2)=2h^2=2(\tfrac14)^2=\tfrac18$. Centering changes the mean of the squared-increment term to zero; it does not remove its random fluctuation.

Common Pitfall: Forgetting to Center the Square

A tempting but incorrect update replaces $(\Delta B_k)^2-h$ by $(\Delta B_k)^2$. This leaves a nonzero conditional mean in the purported correction, so it is no longer the Milstein method used in the theorem. For the last example, its uncentered term has mean $1\cdot\mathbb{E}[Z^2]=\tfrac14$, whereas the correct term has mean zero. The missing subtraction is not cosmetic: it is the Ito correction in $\int_t^{t+h}(B_s-B_t)\,dB_s=\tfrac12((\Delta B)^2-h)$.

The order-one conclusion also relies on the smoothness and growth assumptions stated above. In particular, the proof uses bounded derivatives for the Taylor remainder and global Lipschitz bounds for stability. The theorem is a grid-point strong estimate for scalar noise; it does not establish a multidimensional Milstein result.

Check Your Understanding

Use the one-step expansion and the distinction between centered and uncentered terms in your answers.

For $dX_t=a\,dt+\sigma\,dB_t$ with constants $a,\sigma$, explain why the scalar Milstein scheme agrees with the exact solution at every grid point.
Why does an $O(h^{3/2})$ local remainder not by itself yield a global $O(h)$ error by the triangle inequality?
For a Brownian increment with variance $h$, compute $\mathbb{E}[(\Delta B)^2-h]$ and $\mathbb{E}[((\Delta B)^2-h)^2]$.
Identify where the bounded second derivative of $b$ enters the local expansion proof.
State the difference between the theorem’s mean-square maximum bound and a deterministic pathwise error bound.