Ito-Taylor Expansion Behind the Milstein Correction

Why an Ito–Taylor Expansion Adds a Term

For a scalar stochastic differential equation, Euler–Maruyama freezes both coefficients at the start of a time step. That captures the initial drift and diffusion, but not the change in the diffusion coefficient as the solution moves during the step. The first-order effect of that change is itself driven by Brownian motion. Integrating it against Brownian motion produces an iterated stochastic integral, and its value is the source of the scalar Milstein correction.

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 autonomous scalar SDE $dX_t=\mu(X_t)\,dt+\sigma(X_t)\,dB_t$, with $X_t=x$ at the start of the step. Assume for now that a strong solution exists: $X$ is adapted to the given filtration and satisfies the integral form of the SDE almost surely. We make no uniqueness claim here. On every finite time interval, assume $\mu,\sigma$ are globally Lipschitz and bounded, and that $\sigma$ has bounded first and second derivatives. These assumptions ensure the estimates below apply. Each stochastic integral is an Ito integral with a progressively measurable integrand $H_s$ in $\mathcal{L}^2$, meaning $\mathbb{E}\int_0^T H_s^2\,ds<\infty$. The iterated integrals below also have square-integrable integrands.

Write $h>0$, $\Delta B=B_{t+h}-B_t$, and use the multiplication rules $dB_t\,dB_t=dt$, $dB_t\,dt=0$, and $dt\,dt=0$. These rules encode quadratic variation, $[B]_t=t$; they are not the ordinary chain rule. The integral equation for the solution is

X_{t+h}=x+\int_t^{t+h}\mu(X_s)\,ds+\int_t^{t+h}\sigma(X_s)\,dB_s.

Expanding the Coefficients at the Start of the Step

To expose the first omitted term, apply Ito’s formula to the function $f(s,y)=\sigma(y)$ along $X_s$. Its partial derivatives are $f_s=0$, $f_y=\sigma'(y)$, and $f_{yy}=\sigma''(y)$. The formula, including its quadratic-variation correction, gives

d\sigma(X_s) =\left(\mu(X_s)\sigma'(X_s)+\frac12\sigma(X_s)^2\sigma''(X_s)\right)ds +\sigma(X_s)\sigma'(X_s)\,dB_s.

Define $L^1\sigma(y)=\sigma(y)\sigma'(y)$. Integrate this identity from $t$ to $s$, then substitute it into the stochastic integral for $X_{t+h}$. The term involving $L^1\sigma(x)$ is integrated twice against Brownian motion. Define that iterated integral by

I_{11}(t,t+h)=\int_t^{t+h}\left(\int_t^s dB_u\right)dB_s.

The inner integral is adapted, and its squared expectation at time $s$ is $s-t$, so the outer integrand is in $\mathcal{L}^2$ on the step. The Ito integral is a square-integrable martingale: its increments have conditional expectation zero given the filtration at the starting time, as follows first for adapted step integrands from Brownian independent increments and then for $\mathcal{L}^2$ integrands by the Ito isometry and $L^2$ approximation. Thus these integrals are true martingales, not merely local martingales.

Theorem: Under the stated boundedness and smoothness assumptions, for $t\le T-h$, $$ X_{t+h}=x+\mu(x)h+\sigma(x)\Delta B+L^1\sigma(x)I_{11}(t,t+h)+R_{t,h}, $$ where $L^1\sigma(x)=\sigma(x)\sigma'(x)$ and $\mathbb{E}[R_{t,h}^2]\le C h^3$. The finite constant $C$ may depend on the coefficient bounds and $T$, but not on $t,h$, or the starting value $x$.

Proof. Let $\delta_s=X_s-x$, $M_s=\int_t^s\sigma(X_u)\,dB_u$, and let $M_\sigma=\sup_y|\sigma(y)|$, $M_\mu=\sup_y|\mu(y)|$, and $L_\sigma=\sup_y|\sigma'(y)|$. The integral equation and Ito isometry give, for $r=s-t\le T$, $\mathbb{E}[\delta_s^2]\le2M_\mu^2r^2+2M_\sigma^2r\le C_1r$. The drift part of the remainder is $D=\int_t^{t+h}(\mu(X_s)-\mu(x))\,ds$. If $L_\mu$ is a Lipschitz constant for $\mu$, the integral form of the Cauchy–Schwarz inequality gives $\mathbb{E}[D^2]\le h\int_t^{t+h}L_\mu^2\mathbb{E}[\delta_s^2]\,ds\le C_2h^3$.

For the diffusion part, define $U_s=\delta_s-\sigma(x)(B_s-B_t)$. Subtracting the frozen Brownian term from the integral equation yields $U_s=\int_t^s\mu(X_u)\,du+\int_t^s(\sigma(X_u)-\sigma(x))\,dB_u$. By Cauchy–Schwarz, Ito isometry, and the Lipschitz bound for $\sigma$, $\mathbb{E}[U_s^2]\le2M_\mu^2r^2+2L_\sigma^2\int_t^s\mathbb{E}[\delta_u^2]\,du\le C_3r^2$.

We also need a fourth-moment bound for $\delta_s$. Ito’s formula applied to $f(y)=y^4$, whose derivatives are $f'(y)=4y^3$ and $f''(y)=12y^2$, gives for $M_s$ the correction term $\tfrac12f''(M_s)\sigma(X_s)^2ds=6M_s^2\sigma(X_s)^2ds$. Taking expectations, the Ito integral term has mean zero, and Ito isometry gives $\mathbb{E}[M_u^2]\le M_\sigma^2(u-t)$. Consequently $\mathbb{E}[M_s^4]\le6M_\sigma^2\int_t^s\mathbb{E}[M_u^2]\,du\le3M_\sigma^4r^2$. The drift increment is bounded in absolute value by $M_\mu r$, and $|a+b|^4\le8(|a|^4+|b|^4)$; hence $\mathbb{E}[\delta_s^4]\le C_4r^2$.

Let $L_{\sigma,2}=\sup_y|\sigma''(y)|$. Taylor’s formula with bounded second derivative gives $|\sigma(X_s)-\sigma(x)-\sigma'(x)\delta_s|\le\tfrac12L_{\sigma,2}\delta_s^2$. Therefore $\sigma(X_s)-\sigma(x)-L^1\sigma(x)(B_s-B_t)$ is the sum of $\sigma'(x)U_s$ and that Taylor remainder. The bounds on $\mathbb{E}[U_s^2]$ and $\mathbb{E}[\delta_s^4]$ show that the expectation of its square is at most $C_5r^2$. By Ito isometry, the squared $L^2$ norm of the diffusion remainder $\int_t^{t+h}[\sigma(X_s)-\sigma(x)-L^1\sigma(x)(B_s-B_t)]\,dB_s$ is at most $C_5\int_0^h r^2\,dr=C_5h^3/3$. Combining this with the drift bound proves $\mathbb{E}[R_{t,h}^2]\le Ch^3$, after increasing $C$ if needed. $\square$

The Iterated Integral Is a Centered Square

Proposition: For every $h>0$, $$ I_{11}(t,t+h)=\frac12\big((\Delta B)^2-h\big) $$ almost surely.

Proof. Apply Ito’s formula to $f(y)=y^2$, with $f'(y)=2y$ and $f''(y)=2$, along $Y_s=B_s-B_t$. Since $dY_s=dB_s$, the explicit correction is $\tfrac12 f''(Y_s)d[Y]_s=ds$, because $[B]_s-[B]_t=s-t$. Integrating from $t$ to $t+h$ gives $(\Delta B)^2=2\int_t^{t+h}(B_s-B_t)\,dB_s+h$. The integral is $I_{11}(t,t+h)$ by its definition, so rearranging proves the identity. $\square$

Substitution into the theorem gives the expansion $X_{t+h}=x+\mu(x)h+\sigma(x)\Delta B+\tfrac12\sigma(x)\sigma'(x)((\Delta B)^2-h)+R_{t,h}$. Euler–Maruyama keeps the first three terms and omits the iterated-integral term. The correction is not an arbitrary adjustment: it is exactly the leading contribution from the diffusion coefficient changing within the step. Since $\mathbb{E}[(\Delta B)^2]=h$, this correction has conditional mean zero given $\mathcal{F}_t$, but it is generally nonzero on an individual step.

Worked Example: Constant Diffusion Has No Correction

Take $dX_s=3\,ds+2\,dB_s$, starting at $x=1$, over a step of length $h=\tfrac14$. Here $\mu(x)=3$, $\sigma(x)=2$, and $\sigma'(x)=0$, so $L^1\sigma(x)=2\cdot0=0$. The expansion through the iterated-integral term is therefore $X_{t+h}=1+3(\tfrac14)+2\Delta B+R_{t,h}$. The retained terms give $1+\tfrac34+2\Delta B=\tfrac74+2\Delta B$.

Check: The SDE has the exact solution $X_{t+h}=1+3h+2\Delta B$, obtained by integrating its constant coefficients. Substitution of $h=\tfrac14$ gives $1+\tfrac34+2\Delta B=\tfrac74+2\Delta B$, matching the expansion with zero remainder. Constant diffusion has no within-step diffusion variation to produce a correction.

Worked Example: Geometric Diffusion and the Iterated Term

Consider $dX_s=\tfrac14X_s\,ds+\tfrac12X_s\,dB_s$, with starting value $x=2$, step $h=\tfrac14$, and use the increment $\Delta B=\tfrac14$ to evaluate the algebraic terms. This increment is an input for the calculation, not a positive-probability conditioning event. The coefficients at $x=2$ are $\mu(x)=\tfrac12$, $\sigma(x)=1$, and $\sigma'(x)=\tfrac12$, so $L^1\sigma(x)=1\cdot\tfrac12=\tfrac12$. The Euler–Maruyama terms sum to $2+(\tfrac12)(\tfrac14)+1(\tfrac14)=2+\tfrac18+\tfrac14=\tfrac{19}{8}$. The iterated integral is $I_{11}=\tfrac12((\tfrac14)^2-\tfrac14)=\tfrac12(\tfrac1{16}-\tfrac4{16})=-\tfrac3{32}$. Thus the correction is $(\tfrac12)(-\tfrac3{32})=-\tfrac3{64}$, and the expansion through this term is $\tfrac{19}{8}-\tfrac3{64}=\tfrac{152}{64}-\tfrac3{64}=\tfrac{149}{64}$, apart from the remainder.

Check: Substituting the increment back into the centered-square expression gives $\tfrac12((\Delta B)^2-h)=\tfrac12(\tfrac1{16}-\tfrac4{16})=-\tfrac3{32}$, the same iterated integral. Subtracting the correction from the corrected value gives $\tfrac{149}{64}+\tfrac3{64}=\tfrac{152}{64}=\tfrac{19}{8}$, the Euler–Maruyama value.

Worked Example: A Nonlinear Diffusion Coefficient

Take $dX_s=(1+X_s)\,dB_s$ with $x=1$, $h=\tfrac12$, and algebraic increment $\Delta B=\tfrac12$. Here $\mu=0$, $\sigma(1)=2$, and $\sigma'(1)=1$, so $L^1\sigma(1)=2$. Euler–Maruyama gives $1+2(\tfrac12)=2$. The iterated integral is $I_{11}=\tfrac12((\tfrac12)^2-\tfrac12)=\tfrac12(\tfrac14-\tfrac24)=-\tfrac18$. The correction is $2(-\tfrac18)=-\tfrac14$, so the expansion through this term gives $2-\tfrac14=\tfrac74$.

Check: The coefficient of the centered square is $\tfrac12\sigma(1)\sigma'(1)=\tfrac12(2)(1)=1$. Directly, its contribution is $1((\tfrac12)^2-\tfrac12)=\tfrac14-\tfrac12=-\tfrac14$, agreeing with $L^1\sigma(1)I_{11}=2(-\tfrac18)=-\tfrac14$.

Interpretation and a Common Pitfall

The remainder estimate in the theorem says its root-mean-square size is at most of order $h^{3/2}$ under the stated assumptions. It is a local expansion, not by itself a proof of a global convergence rate for a time-stepping method. Such a rate requires controlling how local errors accumulate across steps. Nor does the estimate say that the remainder vanishes on each path or that every individual Brownian increment produces a large correction.

Common mistake: dropping the $-h$ in the centered square. A tempting but incorrect calculation is to replace $I_{11}$ by $\tfrac12(\Delta B)^2$, yielding a correction $\tfrac12\sigma(x)\sigma'(x)(\Delta B)^2$. The proposition shows that this leaves out the term $-\tfrac12h\sigma(x)\sigma'(x)$. For example, in the geometric-diffusion calculation above, that incorrect expression would give $\tfrac12(\tfrac12)(\tfrac1{16})=\tfrac1{64}$. The correct correction is $-\tfrac3{64}$, because the increment-square term must be centered by $h=\tfrac14$. The missing term also changes the conditional mean: the correct correction has conditional mean zero, whereas the uncentered expression has conditional mean $\tfrac12\sigma(x)\sigma'(x)h$.

The key distinction is between a single squared increment and its conditional mean. The identity $dB_t\,dB_t=dt$ records quadratic variation; it does not assert $(\Delta B)^2=h$ on every step. Retaining the centered difference is what preserves the actual random iterated integral.

Check Your Understanding

Use the definitions and the proved identities above to answer each question.

For $\sigma(x)=c$, where $c$ is constant, compute $L^1\sigma(x)$. What does this imply about the iterated-integral term?
For $\sigma(x)=x$ and starting value $x=3$, compute $L^1\sigma(x)$.
Show from the proposition and $\mathbb{E}[(\Delta B)^2\mid\mathcal{F}_t]=h$ that the iterated-integral correction has conditional mean zero.
For $h=\tfrac13$ and $\Delta B=\tfrac23$, compute $I_{11}$ exactly.
Why does an $L^2$ remainder bound of order $h^{3/2}$ not alone prove a global strong convergence rate for a numerical scheme?

Ito-Taylor Expansion Behind the Milstein Correction

What You'll Learn