Multilevel Monte Carlo For Estimating the Expected Value of Stochastic Differential Equations

Motivation

A few days ago, I replied to a tweet about stochastic differential equations (SDEs).
(Insert screenshot of tweet + reply)

Quant Beckman’s response was a bit snarky—especially considering it was his own post—but it did get me thinking seriously about the computational challenges behind SDEs. Problems like computing expectations, option payoffs, or tail probabilities can easily become expensive, especially when the SDE has no closed-form solution.

But computationally expensive doesn’t mean hopeless.
It just means we need better algorithms.

This led me down the path of Multilevel Monte Carlo (MLMC), a method designed to dramatically speed up the estimation of expectations of SDEs by combining simulations at different time resolutions. Along the way, I also examined how approximate random variables, single vs double precision, and compensated summation (Kahan) affect accuracy.

Background: Multilevel Monte Carlo (MLMC)

1. The Problem

Consider the Ornstein–Uhlenbeck (OU) SDE:

d x_{t} = θ (μ - x_{t}) d t + σ d W_{t}

What we can say from this equation is that if we simulated the SDE from time $0$ to $T$ , the value $X_{T}$ is the random state of the process at the terminal time $T$ . Due to the stochastic nature of stochastic differential equations, $X_{T}$ is non determinstic, therefore it is oftentimes a hard problem to figure out how to estimate the expectation of the distribution

𝔼 [P (X_{T})]

where $P (X_{T})$ for us is just the distribution of $X_{T}$ . (Oftentimes, for instance, in options pricing $P (x) = m a x (x - K, 0)$ .

For many SDEs (including OU with general payoffs), computing this expectation analytically is difficult or impossible. So we turn to numerical simulation.

A simple and widely used scheme is the Euler–Maruyama method.

1.1 Discretizing the OU SDE with Euler–Maruyama

Given

d x_{t} = θ (μ - x_{t}) d t + σ d W_{t}, x (0) = x_{0}

integrate over the interval ([t_n, t_{n+1}]) with step size (\Delta t):

x_{n + 1} - x_{n} = \int_{t_{n}}^{t_{n + 1}} θ (μ - x_{s}) d s + \int_{t_{n}}^{t_{n + 1}} σ d W_{s} .

Euler–Maruyama approximates these integrals using:

drift evaluated at the left endpoint:

$θ (μ - x_{n}) Δ t \$
Brownian increment:

$Δ W_{n} ~ 𝒩 (0, Δ t)$

giving the update

x_{n + 1} = x_{n} + θ (μ - x_{n}) Δ t + σ \sqrt{Δ t} Z_{n}, Z_{n} ~ 𝒩 (0, 1) .

This produces a discrete approximation $P_{ℓ}$ with step size $Δ t_{ℓ} = T / 2^{ℓ}$ .

But refining the time step improves accuracy only at the cost of drastically more computation.

This is where MLMC enters.

2. The MLMC Idea

Instead of running a huge Monte Carlo simulation at the finest step size, MLMC decomposes the expectation using a telescoping sum:

𝔼 [P_{L}] = 𝔼 [P_{0}] + \sum_{ℓ = 1}^{L} 𝔼 [P_{ℓ} - P_{ℓ - 1}],

where:

$(P_{ℓ})$ is the payoff computed using time step $(Δ t_{ℓ} = T / 2^{ℓ})$ ,
$(P_{ℓ} - P_{ℓ - 1})$ is a level correction, computed by coupling coarse and fine simulations using the same random normal increments.

The key insight:

Coarse levels are cheap but inaccurate.
Fine levels are expensive but accurate.
MLMC mixes many cheap coarse samples with few expensive fine samples—
reducing total variance at minimal cost.

3. What Quantities MLMC Needs

To carry out the MLMC optimization, we need just three quantities per level $(ℓ)$ :

Variance

$V_{ℓ} = V a r (P_{ℓ} - P_{ℓ - 1})$

→ Determines how many samples we should take at each level.

Cost

$C_{ℓ} \propto 2^{ℓ}$

→ Number of time steps per sample.

Bias estimate

B_{L} = | 𝔼 [P_{L} - P_{L - 1}] |

→ Determines how deep the hierarchy must extend.

These three numbers feed directly into the MLMC formula for the optimal number of samples per level.