Nested Multilevel Monte Carlo Simulation

Python implementation of the numerical examples from the Master thesis:

Nested Multilevel Monte Carlo Simulation Nicolai Ehrhardt, TU Kaiserslautern, 2018 Supervised by Prof. Dr. Klaus Ritter, AG Computational Stochastics

Overview

This repository contains implementations and visualizations of Multilevel Monte Carlo (MLMC) methods applied to nested simulation problems arising in financial risk measurement, specifically the approximation of Value-at-Risk via nested expectations.

The underlying model is the Black-Scholes model with geometric Brownian motion. The code demonstrates:

• Brownian motion simulation and fine/coarse path coupling
• Euler-Maruyama and Milstein strong approximation schemes for SDEs
• Smoothing functions for indicator function approximation (Chapter 4)
• MLMC estimators for path-independent (European call) and path-dependent (lookback) options
• Nested simulation: outer real-world scenarios + inner risk-neutral pricing
• Convergence rate verification: MSE, level variance, and level mean

Mathematical Background

Nested Expectation Problem

The goal is to compute quantities of the form

$$a = E\big(f(h(X))\big),$$

where $f$ is a suitable mapping and $h$ is a function that depends on the expectation of another random variable $Y$:

$$h(x) = E\big(g(x, Y)\big).$$

In general, $h(x)$ cannot be evaluated exactly and must be approximated via Monte Carlo simulation. This gives rise to a nested simulation: the outer expectation over $X$ and the inner expectation over $Y$ are estimated iteratively.

Application: Value-at-Risk

In the financial setting, let $(S(x,t))_{t \in [T_0, T_1]}$ be an asset process (geometric Brownian motion) with initial value depending on a risk factor $x$. For a payoff function $\Phi$ (e.g. European call), the risk mapping is

$$h(x) = E\big(\Phi(S(x, \cdot))\big),$$

which gives the expected payoff at maturity $T_1$ conditional on the risk scenario $x$. The loss function is $L = C_0 - h(X)$, and the Value-at-Risk at level $\alpha$ is defined as

$$\text{VaR}_\alpha = \arg\min_k P({L \leq k}) \geq \alpha.$$

Thus, estimating VaR reduces to approximating the distribution function of $h(X)$:

$$F(s) = P({h(X) \leq s}) = E\big(\mathbb{1}_{(-\infty, s]}(h(X))\big).$$

This is exactly the nested expectation problem with $f = \mathbb{1}_{(-\infty, s]}$.

Black-Scholes Model

The numerical examples use the one-dimensional Black-Scholes model:

$$\mathrm{d}S(t) = \mu \cdot S(t),\mathrm{d}t + \sigma \cdot S(t),\mathrm{d}W(t), \quad S(0) = s_0,$$

with explicit solution $S(t) = s_0 \exp\big((\mu - \tfrac{\sigma^2}{2})t + \sigma W(t)\big)$ and parameters $s_0 = 100$, $\mu = 0.05$, $\sigma = 0.2$.

The Euler-Maruyama scheme approximates $S$ on a grid with $2^l$ time steps:

$$\hat{S}_{k+1} = \hat{S}_k + \mu,\hat{S}_k,\Delta_k + \sigma,\hat{S}_k,\Delta W_k,$$

while the Milstein method adds a correction term for higher-order strong convergence:

$$\check{S}_{k+1} = \check{S}_k + \mu,\check{S}_k,\Delta_k + \sigma,\check{S}_k,\Delta W_k + \tfrac{1}{2}\sigma^2,\check{S}_k\big((\Delta W_k)^2 - \Delta_k\big).$$

Multilevel Monte Carlo

The MLMC estimator is based on the telescoping decomposition

$$E(R^{(L)}) = E(R^{(0)}) + \sum_{l=1}^{L} E(R^{(l)} - R^{(l-1)}),$$

where $R^{(l)} = \Phi(\hat{S}^{(2^l)}(T_1))$ uses a level-$l$ discretization. Each summand is estimated independently with $M_l$ samples, using coupled fine and coarse simulations driven by the same Brownian path. The variance of $R^{(l)} - R^{(l-1)}$ decays with rate $\beta$, which determines the overall cost of the algorithm (Theorem 3.1).

Smoothing Function

Since the indicator function $\mathbb{1}_{(-\infty, s]}$ is not Lipschitz-continuous, a polynomial smoothing function $\Psi$ is used to approximate it. The polynomial $p$ of degree $r+1$ is constructed to satisfy moment conditions, and $\Psi$ is defined as

$$\Psi(x) = \begin{cases} 1 & \text{if } x < -1, \ p(x) & \text{if } |x| \leq 1, \ 0 & \text{if } x > 1, \end{cases}$$

so that $\Psi\big(\tfrac{t-s}{\delta}\big) \approx \mathbb{1}_{(-\infty,s]}(t)$, where $\delta$ controls the smoothing width and $r$ the polynomial degree.

Highlights

Nested Simulation

Outer real-world scenarios (left of the vertical line) fan out into multiple inner risk-neutral pricing scenarios:

2000 Geometric Brownian Motion Paths

GBM paths with the log-normal terminal distribution emerging at maturity:

Fine and Coarse Path Coupling

The multilevel method exploits coupled fine/coarse discretizations of the same Brownian path:

Smoothing of the Indicator Function

Polynomial approximations of the indicator function for different smoothness parameters r and smoothing widths delta:

Convergence Rates

MSE of Euler-Maruyama (slope -1) vs Milstein (slope -2) confirms the theoretical strong convergence orders:

MLMC level variance and level mean decay at the predicted rates:

Empirical Distribution Function

Approximation of the distribution function of h(X) using the Black-Scholes pricing formula, converging with increasing sample size:

Structure

core.py                      # Core functions (simulation, MLMC, smoothing, Black-Scholes)
plot_brownian_paths.py        # Figures 6 & 7: Independent and fine/coarse Brownian paths
plot_gbm.py                  # Figure 8: GBM Euler vs Milstein approximation
plot_gbm_distribution.py     # Figures 4 & 16: GBM paths with log-normal distribution, empirical CDF
plot_nested_simulation.py    # Figure 9: Nested simulation (outer + inner scenarios)
plot_nested_arrow.py         # Figure 1: Schematic nested simulation diagram
plot_smoothing.py            # Figures 2 & 3: Smoothing of indicator function
plot_mse_convergence.py      # Figure 5: MSE convergence of Euler and Milstein
plot_level_variance.py       # Figures 10 & 11: Level variance and level mean
matlab/                      # Original MATLAB scripts
output/                      # Pre-generated plots

Getting started

pip install -r requirements.txt
python plot_brownian_paths.py
python plot_smoothing.py
python plot_nested_simulation.py
# ...

Key results

The numerical experiments confirm the theoretical convergence rates from the thesis:

Method	Quantity	Observed rate	Theoretical rate
Euler-Maruyama	MSE at terminal time	-1.02	-1 (Thm 5.2)
Milstein	MSE at terminal time	-1.99	-2 (Thm 5.5)
Euler MLMC	Level variance	-0.96	-1 (Rem 5.9)
Milstein MLMC	Level variance	-1.92	-2 (Rem 5.13)

References

• M. B. Giles. Multilevel Monte Carlo methods. Acta Numerica, 2015.
• M. B. Giles, T. Nagapetyan, K. Ritter. Multilevel Monte Carlo Approximation of Distribution Functions and Densities. SIAM/ASA J. Uncertainty Quantification, 2015.
• P. Kloeden, E. Platen. Numerical Solution of Stochastic Differential Equations. Springer, 1999.
• T. Mueller-Gronbach, E. Novak, K. Ritter. Monte Carlo-Algorithmen. Springer, 2012.

License

MIT