SDE discretization and weak convergence, financial engineering for quants lesson

When a price or a variance follows an SDE whose transition law has no closed form, a simulation has to step the equation forward on a time grid. How much does the grid distort a price, and does it matter whether the question concerns each simulated path or only the distribution of outcomes? This lesson derives the Euler-Maruyama step from the integral form of the SDE and the Milstein correction from the next stochastic Taylor term, states the strong and weak orders they reach under explicit hypotheses, shows how kinks, jumps and barriers in a payoff change the weak error, and gives a step of a square-root variance model that goes negative.

Key terms

Euler-Maruyama scheme

Definition Consider the scalar SDE $dX_t = b(t, X_t)\, dt + \sigma(t, X_t)\, dW_t$ on $[0, T]$ with a deterministic initial value $X_0$ , where $W$ is a standard Brownian motion under the measure in which the equation is stated. In integral form, over a step from $t_n$ to $t_{n+1} = t_n + h$ ,

\begin{aligned} &X_{t_{n+1}} - X_{t_n} \\[0.4em] &\quad = \int_{t_n}^{t_{n+1}} b(s, X_s)\, ds \\[0.4em] &\qquad + \int_{t_n}^{t_{n+1}} \sigma(s, X_s)\, dW_s . \end{aligned}

On the uniform grid $t_n = nh$ with $h = T/M$ , the Euler-Maruyama scheme freezes both coefficients at the start of each step. With $X^h_0 = X_0$ ,

\begin{aligned} X^h_{n+1} &= X^h_n + b_n h \\[0.4em] &\quad + \sigma_n \Delta W_n , \end{aligned}

where $b_n = b(t_n, X^h_n)$ , $\sigma_n = \sigma(t_n, X^h_n)$ and $\Delta W_n = W_{t_{n+1}} - W_{t_n}$ . The increments are independent and normal with mean $0$ and variance $h$ , and they are simulated as $\sqrt{h}\, Z_n$ with independent standard normal $Z_n$ . The coefficients are fixed before the increment is seen, as in the Itô integral. Evaluating the diffusion coefficient at the end of the step, after the increment is known, defines a different method, which in general converges to a different limit. Two separate limits are involved in a simulation. Refining $h$ is meant to bring the law of $X^h_T$ closer to the law of $X_T$ , while more sample paths estimate expectations under the law of $X^h_T$ more precisely and cannot remove the discretization bias of a fixed step. Whether refinement converges at all depends on the regularity of the coefficients, the moments of the solution and any constraint on its state.

In plain words Over a short step, Euler pretends that the drift and the volatility keep the values they had at the start of the step. It is the stochastic version of the forward Euler method for ordinary differential equations, with a Gaussian shock whose variance equals the step length.

Example Take constructed geometric Brownian motion inputs, $dS_t = \mu S_t\, dt + \sigma S_t\, dW_t$ with $S_0 = 100$ , $\mu = 0.05$ , $\sigma = 0.2$ and $h = 0.25$ , and one standard normal draw $Z = -1.3$ , so that $\Delta W = \sqrt{0.25} \cdot (-1.3) = -0.65$ . The drift contributes $\mu S_0 h = 1.25$ and the diffusion $\sigma S_0 \Delta W = -13$ , so the Euler step gives $S^h_1 = 100 + 1.25 - 13 = 88.25$ . For the same increment the exact transition gives

\begin{aligned} S_h &= 100\, e^{(\mu - \sigma^2/2) h + \sigma \Delta W} \\[0.4em] &= 100\, e^{0.0075 - 0.13} \\[0.4em] &\approx 88.47059 . \end{aligned}

The drift $\mu$ here is constructed and not a pricing drift. For pricing under $\mathbb{Q}$ the same step would use $r - q$ in place of $\mu$ .

Milstein correction

Definition Keep the scalar SDE and a scalar Brownian motion, and suppose $\sigma$ is smooth enough to apply Itô’s formula to $\sigma(s, X_s)$ , with $\sigma_x = \partial \sigma / \partial x$ its derivative in the state variable. On $[t_n, t_{n+1}]$ the leading random part of $\sigma(s, X_s) - \sigma(t_n, X_{t_n})$ is $\sigma \sigma_x (t_n, X_{t_n}) \, (W_s - W_{t_n})$ , and the other terms contribute only at higher order in $h$ . Itô’s formula applied to $(W_s - W_{t_n})^2$ gives the iterated integral

\begin{aligned} &\int_{t_n}^{t_{n+1}} \bigl(W_s - W_{t_n}\bigr)\, dW_s \\[0.4em] &\quad = \tfrac12\bigl((\Delta W_n)^2 - h\bigr). \end{aligned}

Adding the resulting term to the Euler step gives the scalar Milstein scheme

\begin{aligned} &X^h_{n+1} = X^h_n + b_n h \\[0.4em] &\quad + \sigma_n \Delta W_n \\[0.4em] &\quad + \tfrac12 \sigma_n \sigma'_n \\[0.4em] &\qquad \times \bigl((\Delta W_n)^2 - h\bigr), \end{aligned}

where $\sigma'_n = \sigma_x(t_n, X^h_n)$ . Given $X^h_n$ the correction has mean zero, because $\mathbb{E}[(\Delta W_n)^2] = h$ . When $\sigma$ does not depend on $x$ , the correction vanishes and Milstein coincides with Euler. With several Brownian motions whose diffusion vector fields do not commute, the corresponding expansion also needs the iterated integrals of different components, whose antisymmetric parts, the Lévy areas, are not functions of the increments. Adding a scalar correction to each coordinate separately then does not in general give strong order $1$ .

In plain words Euler ignores the fact that the volatility itself moves during the step in response to the same shock. Milstein adds the first-order effect of that movement. The correction involves the squared shock minus its mean, so it changes each path without shifting the conditional mean of the step.

Example For geometric Brownian motion the diffusion coefficient is $\sigma(t, x) = \sigma x$ with constant volatility $\sigma$ , so $\sigma_x = \sigma$ and the correction is $\tfrac12 \sigma^2 S^h_n \bigl((\Delta W_n)^2 - h\bigr)$ . With the constructed inputs, $(\Delta W)^2 - h = 0.4225 - 0.25 = 0.1725$ and $\tfrac12 \sigma^2 S_0 = 2$ , so the correction is $0.345$ and the Milstein step gives $88.595$ . Against the exact $88.47059$ , the Euler error is about $-0.2206$ and the Milstein error about $0.1244$ . Milstein is closer for this draw, but one draw is not a convergence test.

Strong order of convergence

Definition A scheme has strong order $\gamma > 0$ at time $T$ if there are constants $C$ and $h_0 > 0$ , independent of $h$ , such that

\begin{aligned} &\Bigl(\mathbb{E}\bigl[\lvert X_T - X^h_T \rvert^2\bigr]\Bigr)^{1/2} \\[0.4em] &\quad \le C h^{\gamma} \end{aligned}

for every uniform step $h \le h_0$ , where $X^h_T = X^h_M$ is computed from the increments of the same Brownian path that drives $X$ . Some texts use $\mathbb{E}\lvert X_T - X^h_T \rvert$ or a maximum over grid times instead of this root-mean-square terminal error. The course states the following rates without proof. Let the SDE be scalar and autonomous, $dX_t = b(X_t)\, dt + \sigma(X_t)\, dW_t$ on a fixed horizon $[0, T]$ , with deterministic $X_0$ , and let $b$ and $\sigma$ be six times continuously differentiable with the functions and all their derivatives through order six bounded. These hypotheses make the coefficients globally Lipschitz and provide every moment needed. Then the Euler-Maruyama scheme has strong order $\tfrac12$ and the Milstein scheme strong order $1$ . In particular, Euler has strong order $1$ when $\sigma$ is constant, since it then coincides with Milstein. The hypotheses are sufficient, not necessary. The linear coefficients of geometric Brownian motion are unbounded, yet the same rates hold there by separate growth and moment estimates. Strong order also controls weak error for a test function $f$ with Lipschitz constant $L$ , because $\lvert \mathbb{E}[f(X_T)] - \mathbb{E}[f(X^h_T)] \rvert \le L\, \mathbb{E}\lvert X_T - X^h_T \rvert \le L C h^{\gamma}$ .

In plain words Strong error asks how closely the simulated path follows the true path driven by the same random shocks. It is the relevant error when individual paths matter, for example when coarse and fine simulations of the same path are compared, or when a quantity is computed along a path.

Example For the constructed geometric Brownian motion with $T = 1$ , a Monte Carlo estimate with $200000$ paths compares each scheme with the exact solution on the same Brownian path. Each entry below is a statistical estimate of the root-mean-square terminal error.

Step $h$	Euler error	Milstein error
$1/8$	$1.08$	$0.147$
$1/16$	$0.76$	$0.074$
$1/32$	$0.54$	$0.037$
$1/64$	$0.38$	$0.018$

Halving $h$ divides the Euler error by about $\sqrt{2} \approx 1.41$ and the Milstein error by about $2$ , as orders $\tfrac12$ and $1$ predict.

Weak order of convergence

Definition A scheme has weak order $\beta > 0$ for a class of test functions if for every $f$ in the class there are constants $C_f$ and $h_0 > 0$ with

\begin{aligned} &\bigl\lvert \mathbb{E}[f(X_T)] - \mathbb{E}[f(X^h_T)] \bigr\rvert \\[0.4em] &\quad \le C_f\, h^{\beta} \end{aligned}

for every uniform step $h \le h_0$ . Only the laws of $X_T$ and $X^h_T$ enter, so no coupling of paths is needed. The course states without proof that, under the hypotheses on $b$ and $\sigma$ of the previous card and for a test function $f$ that is six times continuously differentiable with $f$ and its derivatives through order six bounded, the Euler-Maruyama scheme has weak order $1$ . Under hypotheses of the same kind the Milstein scheme also has weak order $1$ , so its gain is in strong error. Because only laws matter, the normal increments may be replaced by other variables with matching low moments, for example $\pm\sqrt{h}$ with probability $\tfrac12$ each, and Euler keeps weak order $1$ under smoothness conditions of the same type. Weak order $1$ does not mean that coupled paths differ by order $h$ , since the strong order of Euler is only $\tfrac12$ .

In plain words Weak error asks only whether the simulated distribution gives the right averages, which is what a price needs. Euler’s paths drift away from the true paths by an amount of order $\sqrt{h}$ , but for smooth payoffs these errors largely cancel in the average, leaving an error of order $h$ .

Example For the constructed geometric Brownian motion, $\mathbb{E}[S^h_{n+1} \mid S^h_n] = S^h_n (1 + \mu h)$ , so the Euler mean is exactly $\mathbb{E}[S^h_T] = S_0 (1 + \mu h)^{T/h}$ , while $\mathbb{E}[S_T] = S_0 e^{\mu T} \approx 105.12711$ . With $T = 1$ the error of the mean is about $0.03258$ for $h = 1/4$ , $0.01636$ for $h = 1/8$ and $0.00820$ for $h = 1/16$ . It halves with $h$ , close to the first-order term

\begin{aligned} &\tfrac12 S_0 e^{\mu T} \mu^2 T h \\[0.4em] &\quad \approx 0.1314\, h . \end{aligned}

Geometric Brownian motion and $f(x) = x$ lie outside the bounded hypotheses, but here the error can be computed exactly. At $h = 1/8$ the error of the mean is more than $60$ times smaller than the root-mean-square path error of about $1.08$ in the previous card.

Payoff regularity

Definition Payoff regularity is the smoothness, as a function of the terminal state or of the path, that a weak error estimate requires of the payoff. The weak order $1$ result above assumes a smooth test function with bounded derivatives. A call payoff $(x - K)^{+}$ is Lipschitz with constant $1$ but has a kink at $K$ , a digital payoff $\mathbf{1}_{\{x > K\}}$ has a jump, and a barrier condition depends on the whole path. For such payoffs the smooth result does not apply as stated, and the rate has to be established by a separate argument. For a Lipschitz payoff the strong bound already gives weak order at least $\tfrac12$ . A barrier that is monitored continuously in the contract but checked only at grid dates in the simulation misses crossings between dates, and for Brownian motion and geometric Brownian motion this monitoring error is of order $\sqrt{h}$ , a known result stated here without proof. Smoothing the payoff and conditional monitoring between dates are remedies, each with its own analysis. For a Brownian motion with constant volatility $s$ , possibly with drift, and a barrier level $H$ above two consecutive grid values $x$ and $y$ , the probability that the path between them crosses $H$ is

\exp\Bigl(-\frac{2(H - x)(H - y)}{s^2 h}\Bigr).

In plain words The argument behind weak order one relies on derivatives of the payoff. At a kink or a jump those derivatives do not exist, so the argument does not go through as it stands, even though a rate may still be proved in another way. A barrier adds a separate problem, because a path checked only at grid dates can cross the barrier and come back unnoticed.

Example Take a constructed log price $\ln(S_t/S_0) = 0.2 W_t$ with no drift and an up barrier at $S_0 e^{0.1}$ , about $110.52$ when $S_0 = 100$ . The barrier is reached by time $1$ exactly when the maximum of $W$ on $[0, 1]$ is at least $0.5$ , which by the reflection principle has probability $2(1 - \Phi(0.5)) \approx 0.6171$ . Checking only at $M$ equally spaced dates, with $W$ sampled exactly there, gives Monte Carlo estimates of about $0.441$ for $M = 4$ , $0.522$ for $M = 16$ , $0.568$ for $M = 64$ and $0.592$ for $M = 256$ , each with standard error below $0.0004$ . The gaps $0.176$ , $0.095$ , $0.049$ and $0.025$ roughly halve each time $M$ is multiplied by four, the signature of an error of order $\sqrt{h}$ , although the path has no error at the grid dates. Adding the bridge crossing probability between dates brings the estimate for $M = 4$ to about $0.617$ .

Square-root state positivity

Definition The Cox-Ingersoll-Ross square-root model for a variance is

\begin{aligned} dv_t &= \kappa(\theta - v_t)\, dt \\[0.4em] &\quad + \xi \sqrt{v_t}\, dW_t \end{aligned}

with $\kappa, \theta, \xi > 0$ and $v_0 \ge 0$ . It has a unique strong solution, and that solution stays nonnegative for all $t$ . If in addition $2\kappa\theta \ge \xi^2$ , the Feller condition, and $v_0 > 0$ , it stays strictly positive almost surely. These facts are stated here without proof, and the coefficient $\sqrt{v}$ is not Lipschitz at $0$ , so the global Lipschitz results of the earlier cards do not apply. Square-root state positivity is this property of the exact solution, never leaving $[0, \infty)$ . The Euler step

\begin{aligned} v^h_{n+1} &= v^h_n + \kappa(\theta - v^h_n)\, h \\[0.4em] &\quad + \xi \sqrt{v^h_n h}\, Z_n \end{aligned}

does not have it. For every $v^h_n > 0$ the step is negative when

Z_n < -\frac{v^h_n + \kappa(\theta - v^h_n) h}{\xi \sqrt{v^h_n h}},

an event of positive probability, and then the next square root is undefined. Truncating at zero, reflecting, sampling the exact noncentral chi-square transition and other positivity-preserving constructions are different numerical methods, each needing its own convergence and bias analysis.

In plain words The exact square-root process has a drift that pulls it up and a shock that shrinks to zero as the variance approaches zero, so it cannot cross below zero. The Euler step sizes its Gaussian shock at the start of the step, and a large enough negative shock carries the variance below zero in one step.

Example Take constructed inputs $v_0 = 0.001$ , $\kappa = 1$ , $\theta = 0.04$ , $\xi = 0.5$ , $h = 0.1$ and $Z = -3$ . The drift adds $\kappa(\theta - v_0) h = 0.0039$ and the diffusion adds $\xi \sqrt{v_0 h}\, Z = 0.5 \cdot 0.01 \cdot (-3) = -0.015$ , so the Euler value is $0.001 + 0.0039 - 0.015 = -0.0101$ . Here $2\kappa\theta = 0.08$ is below $\xi^2 = 0.25$ , so the exact process can reach zero, but it never goes below it. From this state the threshold is $-0.0049/0.005 = -0.98$ , so the Euler step is negative with probability $\Phi(-0.98) \approx 0.164$ , about one step in six.

Check yourself

For constructed geometric Brownian motion inputs $S_0 = 50$ , $\mu = 0.04$ , $\sigma = 0.3$ , $h = 0.04$ and $Z = 1.5$ , compute one Euler step, one Milstein step and the exact transition with the same increment.

Show the answer

The increment is $\Delta W = \sqrt{0.04} \cdot 1.5 = 0.3$ . The drift term is $0.04 \cdot 50 \cdot 0.04 = 0.08$ and the diffusion term is $0.3 \cdot 50 \cdot 0.3 = 4.5$ , so the Euler value is $54.58$ . The Milstein correction is

\begin{aligned} &\tfrac12 \cdot 0.09 \cdot 50 \cdot (0.09 - 0.04) \\[0.4em] &\quad = 2.25 \cdot 0.05 = 0.1125, \end{aligned}

which gives $54.6925$ . The exact exponent is $(0.04 - 0.045) \cdot 0.04 + 0.3 \cdot 0.3 = 0.0898$ , so the exact value is $50\, e^{0.0898} \approx 54.69777$ . The Euler error is about $-0.1178$ and the Milstein error about $-0.0053$ . Both numbers come from one draw and say nothing by themselves about the orders of convergence.

A colleague prices a contract on the constructed geometric Brownian motion with Euler steps of $h = 1/64$ and a million paths. Knowing that the root-mean-square path error at this step is about $0.38$ , the colleague concludes that the estimated mean of $S_T$ is off by about $0.38$ and that ten times more paths will fix it. What is wrong with both conclusions?

Show the answer

The value $0.38$ is a strong error, the typical distance between coupled Euler and exact paths, and a mean is a weak quantity. The Euler mean is exactly $S_0(1 + \mu h)^{1/h}$ , and its error is

\begin{aligned} &100\bigl(e^{0.05} - (1 + 0.05/64)^{64}\bigr) \\[0.4em] &\quad \approx 0.00205, \end{aligned}

about $185$ times smaller than $0.38$ . More paths reduce only the sampling error of the estimate of $\mathbb{E}[S^h_T]$ . The discretization bias of about $0.00205$ stays the same for every number of paths, and only a finer grid, or here the exact transition, removes it.

Parameters $\kappa = 2$ , $\theta = 0.04$ and $\xi = 0.3$ satisfy the Feller condition, and a modeler argues that Euler steps of $h = 1/52$ therefore cannot produce a negative variance. Check the condition, then find the probability that one Euler step from $v^h_n = 0.01$ is negative.

Show the answer

The condition holds, since $2\kappa\theta = 0.16 \ge \xi^2 = 0.09$ . It is a statement about the exact process, which then stays strictly positive, and it says nothing about the scheme. From $v^h_n = 0.01$ the deterministic part is $0.01 + 2 \cdot 0.03/52 \approx 0.011154$ and the shock scale is $0.3\sqrt{0.01/52} \approx 0.0041603$ , so the step is negative when

\begin{aligned} Z_n &< -\frac{0.011154}{0.0041603} \\[0.4em] &\approx -2.681 , \end{aligned}

with probability $\Phi(-2.681) \approx 0.0037$ . That is small for one step, but a simulation makes many steps on many paths. In a run of $200000$ paths over five years of weekly steps from $v_0 = 0.04$ , continued with the positive part inside the coefficients, about $30$ percent of the paths produced at least one negative Euler value.

SDE discretization and weak convergence

Key terms

Euler-Maruyama scheme

Milstein correction

Strong order of convergence

Weak order of convergence

Payoff regularity

Square-root state positivity

Check yourself

How this lesson connects

Builds on

Used later in