Multilevel Monte Carlo method

Remarks

The telescoping sum $\mathbb E P_L = \sum_{\ell=0}^{L} \mathbb E[P_\ell - P_{\ell-1}]$ shows that $Y$ estimates $\mathbb E P_L$ without bias, so its mean squared error is the squared bias $(\mathbb E[P_L - P])^2$ plus $\sum_\ell \operatorname{Var} Y_\ell$. One chooses $L$ by clause 1 so that the bias is below $\varepsilon/\sqrt2$. The $N_\ell$ then minimise the cost under the constraint $\sum_\ell \operatorname{Var} Y_\ell \le \varepsilon^2/2$. A Lagrange multiplier gives $N_\ell$ proportional to $\sqrt{V_\ell / c_\ell}$, where $V_\ell$ and $c_\ell$ are the variance and cost of one sample on level $\ell$. The variance decay in clause 3 comes from coupling: $P_\ell$ and $P_{\ell-1}$ are computed from the same random input, so their difference is small. Example: consider $\mathbb E\, \varphi(S_T)$, where $S$ solves a stochastic differential equation and $\varphi$ is Lipschitz. Under standard conditions, the Euler–Maruyama scheme with $2^\ell$ time steps gives $\alpha = 1$, $\beta = 1$ and $\gamma = 1$. So the cost is $O(\varepsilon^{-2} (\ln \varepsilon)^2)$, against $O(\varepsilon^{-3})$ for plain Monte Carlo with the same scheme. In the case $\beta > \gamma$ the cost has the same order $\varepsilon^{-2}$ as plain Monte Carlo with exact samples of $P$ at unit cost.

Going deeper

For a student

The method, introduced by Giles for stochastic differential equations, is related to earlier multilevel work of Heinrich on parametric integration. It is a recursive variance reduction: each level acts as a control variate for the next. The three regimes of the theorem say whether the cost is dominated by the coarsest levels ($\beta > \gamma$), spread evenly ($\beta = \gamma$), or dominated by the finest level ($\beta < \gamma$). The Milstein scheme improves the strong convergence and gives $\beta = 2$ in many scalar problems, which moves the Euler–Maruyama case $\beta = \gamma$ to the optimal regime $\beta > \gamma$. The same framework applies to partial differential equations with random coefficients, where the levels are finite element meshes.