The Scalar Milstein Update and Brownian Increment Squares

From the Diffusion Variation to a Simulation Step

The Euler–Maruyama update uses the drift and diffusion coefficients at the start of each time step. The scalar Milstein update adds one term: it accounts for how the diffusion coefficient changes as the state moves during that step. That term involves an iterated Brownian integral, but in one dimension the integral can be computed from the same Brownian increment already used by Euler–Maruyama. The key identity is that the iterated integral equals half the squared increment minus half the step length.

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 Itô SDE $dX_t=\mu(t,X_t)\,dt+\sigma(t,X_t)\,dB_t$, with initial value $X_0$. Here $\mu$ is the drift coefficient and $\sigma$ is the diffusion coefficient. Assume, for example, that the coefficients are continuous in time, globally Lipschitz in the state, and of linear growth: for some constant $C$, $|\mu(t,x)-\mu(t,y)|+|\sigma(t,x)-\sigma(t,y)|\le C|x-y|$ and $|\mu(t,x)|+|\sigma(t,x)|\le C(1+|x|)$. Let $X_0$ be square-integrable and adapted. These standard conditions ensure a unique strong solution with finite second moments on bounded time intervals. In particular, the stochastic integral in the SDE is an Itô integral with progressively measurable integrand $\sigma(s,X_s)$ satisfying $\mathbb{E}\int_0^T \sigma(s,X_s)^2\,ds<\infty$.

Definition: On a grid $0=t_0<t_1<\cdots<t_N=T$, write $h_n=t_{n+1}-t_n$, $\Delta B_n=B_{t_{n+1}}-B_{t_n}$, and $X_n$ for the numerical approximation at $t_n$. The scalar Milstein update is $$ X_{n+1}=X_n+\mu(t_n,X_n)h_n+\sigma(t_n,X_n)\Delta B_n +\frac12\sigma(t_n,X_n)\,\partial_x\sigma(t_n,X_n)\big((\Delta B_n)^2-h_n\big). $$

The update is for an Itô SDE. In Stratonovich notation the drift is different: in one dimension, $dX_t=a(t,X_t)\,dt+\sigma(t,X_t)\circ dB_t$ converts to the Itô drift $\mu=a+\tfrac12\sigma\,\partial_x\sigma$ under the usual differentiability conditions. The Milstein correction above belongs to the Itô formulation; it is not a license to use the same drift in both conventions.

The Iterated Integral Identity

Fix a step from $t$ to $t+h$, and define $$ I_{t,t+h}=\int_t^{t+h}\left(\int_t^s dB_u\right)dB_s =\int_t^{t+h}(B_s-B_t)\,dB_s. $$

These are Itô integrals. The inner integrand is the constant process $1$, which is progressively measurable and has $\mathbb{E}\int_t^{t+h}1^2\,du=h<\infty$. For the outer integral, $B_s-B_t$ is adapted and progressively measurable, and $\mathbb{E}\int_t^{t+h}(B_s-B_t)^2\,ds=\int_t^{t+h}(s-t)\,ds=h^2/2<\infty$. Thus both integrals are defined in the square-integrable Itô sense.

Theorem: For every $h\ge0$, $$ I_{t,t+h}=\frac12\big((B_{t+h}-B_t)^2-h\big)\quad\text{almost surely}. $$

Proof. Apply Itô's formula to $f(s,x)=(x-B_t)^2$, with $x=B_s$ for $t\le s\le t+h$. Its derivatives are $f_s=0$, $f_x=2(x-B_t)$, and $f_{xx}=2$. The Itô correction is explicitly $\tfrac12 f_{xx}\,ds=ds$, since $dB_s\,dB_s=ds$. Therefore $$ (B_{t+h}-B_t)^2 =2\int_t^{t+h}(B_s-B_t)\,dB_s+\int_t^{t+h}ds =2I_{t,t+h}+h. $$ Rearranging gives the claimed identity. If $h=0$, both sides are zero, so the degenerate step is included. $\square$

The subtraction of $h$ is the effect of Brownian quadratic variation. In differential notation, $dB_t\,dB_t=dt$, while $dB_t\,dt=0$ and $dt\,dt=0$. The ordinary chain rule would omit the $\tfrac12 f_{xx}\,dt$ term and therefore would give the wrong iterated integral.

How the Identity Produces the Milstein Correction

Over one step, freeze the coefficients at the left endpoint $(t_n,X_n)$. The leading state change due to the noise is $\sigma(t_n,X_n)(B_s-B_{t_n})$. Expanding the diffusion coefficient to first order in that state change gives a contribution $\partial_x\sigma(t_n,X_n)\sigma(t_n,X_n)(B_s-B_{t_n})$ inside the stochastic integral. Integrating this contribution against $dB_s$ produces $\sigma\,\partial_x\sigma\,I_{t_n,t_{n+1}}$. The iterated-integral identity then gives the correction in the definition.

This is the scalar Milstein truncation of the Itô–Taylor expansion: terms of higher order in the step are omitted. It is not an exact solution formula. The one-dimensional identity is especially convenient because it reduces the needed iterated integral to $\Delta B_n$ alone. No additional random variable is required. For time-dependent coefficients, the displayed step still uses their left-endpoint values; the omitted time-variation terms are part of the higher-order remainder under the usual smoothness assumptions.

Worked Example: A Geometric-Type Step

Take $\mu(x)=\tfrac14x$, $\sigma(x)=\tfrac12x$, $X_n=2$, $h_n=\tfrac14$, and $\Delta B_n=\tfrac14$. Compute both the Euler–Maruyama and Milstein updates.

The drift contribution is $\mu(2)h_n=(\tfrac14\cdot2)(\tfrac14)=\tfrac18$. The diffusion contribution is $\sigma(2)\Delta B_n=(\tfrac12\cdot2)(\tfrac14)=\tfrac14$. Thus Euler–Maruyama gives $2+\tfrac18+\tfrac14=\tfrac{19}{8}$.

For the correction, $\partial_x\sigma=\tfrac12$, so $\tfrac12\sigma(2)\partial_x\sigma=\tfrac12\cdot1\cdot\tfrac12=\tfrac14$. Also $(\Delta B_n)^2-h_n=\tfrac1{16}-\tfrac14=-\tfrac3{16}$. The correction is therefore $-\tfrac3{64}$, and the Milstein value is $2+\tfrac18+\tfrac14-\tfrac3{64}=\tfrac{149}{64}$.

Check: the increment-square identity gives $I=\tfrac12(\tfrac1{16}-\tfrac14)=-\tfrac3{32}$. Multiplying by $\sigma(2)\partial_x\sigma=1\cdot\tfrac12=\tfrac12$ gives $-\tfrac3{64}$, exactly the correction used above.

Worked Example: Nonlinear Diffusion

Let $\mu(x)=0$, $\sigma(x)=1+x$, and take $X_n=1$, $h_n=\tfrac12$, and $\Delta B_n=\tfrac12$. Find the Milstein update and isolate the part that Euler–Maruyama misses.

At $X_n=1$, $\sigma(1)=2$ and $\partial_x\sigma(1)=1$. The Euler diffusion term is $2\cdot\tfrac12=1$, so Euler–Maruyama gives $1+1=2$. The squared-increment adjustment is $(\Delta B_n)^2-h_n=\tfrac14-\tfrac12=-\tfrac14$. The Milstein coefficient is $\tfrac12\cdot2\cdot1=1$, so the correction is $-\tfrac14$, and the update is $2-\tfrac14=\tfrac74$.

Check: $I=\tfrac12(-\tfrac14)=-\tfrac18$. The diffusion variation multiplies this integral by $\sigma(1)\partial_x\sigma(1)=2$, giving $2(-\tfrac18)=-\tfrac14$, as in the update.

Sampling and Conditional Checks

For simulation, take independent standard normal random variables $Z_n$ and set $\Delta B_n=\sqrt{h_n}Z_n$. This has the required Brownian increment distribution $N(0,h_n)$. Substitution yields the directly implementable form $$ X_{n+1}=X_n+\mu(t_n,X_n)h_n+\sigma(t_n,X_n)\sqrt{h_n}Z_n +\frac12\sigma(t_n,X_n)\partial_x\sigma(t_n,X_n)h_n(Z_n^2-1). $$ The same $Z_n$ must be used in both the diffusion term and the squared-increment correction: both arise from the same Brownian increment.

Proposition: If $h>0$, then $\mathbb{E}[I_{t,t+h}\mid\mathcal{F}_t]=0$ and $\mathbb{E}[I_{t,t+h}^2\mid\mathcal{F}_t]=h^2/2$. In particular, the iterated-integral increment is square-integrable and has conditional mean zero.

Proof. Brownian increments after $t$ are independent of $\mathcal{F}_t$, and $\Delta B=B_{t+h}-B_t$ has distribution $N(0,h)$. Thus $\mathbb{E}[\Delta B^2\mid\mathcal{F}_t]=h$, so the identity already proved gives conditional mean zero. For a centered normal variable with variance $h$, $\mathbb{E}[\Delta B^4]=3h^2$. Hence $$ \mathbb{E}[I_{t,t+h}^2\mid\mathcal{F}_t] =\frac14\mathbb{E}\big[(\Delta B^2-h)^2\mid\mathcal{F}_t\big] =\frac14(3h^2-2h^2+h^2)=\frac{h^2}{2}. $$ The conditional second moment is finite, proving square-integrability. At $h=0$, $I=0$, and both conditional moments are zero. $\square$

Worked Example: Checking the Correction's Scale

Suppose the coefficient $\sigma(t_n,X_n)\partial_x\sigma(t_n,X_n)$ equals $2$ and $h_n=\tfrac14$. Compute the conditional mean and variance of the Milstein correction.

The correction is $2I_{t_n,t_{n+1}}$. Its conditional mean is $2\cdot0=0$. Its conditional variance is its conditional second moment because its conditional mean is zero: $4\cdot(h_n^2/2)=4\cdot((\tfrac14)^2/2)=4\cdot\tfrac1{32}=\tfrac18$.

Check: using the normal representation, the correction is $\tfrac12\cdot2\cdot\tfrac14(Z_n^2-1)=\tfrac14(Z_n^2-1)$. Since $\mathbb{E}[Z_n^2]=1$ and $\mathbb{E}[Z_n^4]=3$, its mean is zero and its variance is $\tfrac1{16}(3-2+1)=\tfrac18$, in agreement.

Common Mistake: Dropping the Time-Step Subtraction

A tempting but incorrect replacement is $I_{t,t+h}=\tfrac12(\Delta B)^2$. It misses the Itô correction: the correct expression is $\tfrac12((\Delta B)^2-h)$. For instance, if $h=\tfrac12$ and $\Delta B=0$, the incorrect formula gives $0$, whereas the correct iterated integral value is $-\tfrac14$. A realized Brownian increment can be zero even though the quadratic-variation contribution over the interval is nonzero.

Another common error is to omit the factor $\tfrac12$ in the Milstein coefficient. The coefficient comes from the exact identity for the iterated integral, not just from differentiating $\sigma$. In code, a reliable check is to compare the correction computed as $\tfrac12\sigma\,\partial_x\sigma((\Delta B)^2-h)$ with $\sigma\,\partial_x\sigma I$, using $I=\tfrac12((\Delta B)^2-h)$.

The correction is not an independent noise source. It is a centered function of the same Brownian increment that drives the Euler diffusion term. The subtraction of $h$ centers the squared increment conditionally on the current information, while its size remains relevant to the stepwise approximation.

Check Your Understanding

Use the increment identity and the update formula; assume each Brownian increment is sampled with variance equal to its time-step length.

For $h=\tfrac13$ and $\Delta B=\tfrac13$, compute $I$ exactly.
For $\sigma(x)=3x$, $X_n=1$, $h=\tfrac14$, and $\Delta B=\tfrac12$, compute the Euler–Maruyama and Milstein diffusion contributions, including the correction.
Explain why the squared-increment term and the linear diffusion term must use the same sampled normal variable.
For $h>0$, use the conditional second-moment formula to find $\mathbb{E}[I_{t,t+h}^2\mid\mathcal{F}_t]$.
Identify the missing term in the incorrect formula $I=\tfrac12(\Delta B)^2$, and explain its origin using Itô's formula.