---
title: "Monte Carlo"
book: "Quantitative Methods"
subject: quant
language: en
chapter: 26
exercises: 8
source: https://one-course.com/books/quant/4/en/chapter/26-monte-carlo
---

# Chapter 26 — Monte Carlo

The overnight batch values 40 000 path-dependent trades with 100 000 paths each. The [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) on the book is 80 000 dollars, and the risk manager wants it under 20 000 by tomorrow. Monte Carlo error falls like one over the square root of the number of paths, so a factor of four in error is a factor of sixteen in paths, and sixteen times the machines is not on offer. The alternative is to make each path worth more. This chapter builds the [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) and its generators, then the variance-reduction techniques that change the constant ([antithetic variates](#def-qm-monte-carlo-vr), [control variates](#def-qm-monte-carlo-vr), stratification, [importance sampling](#def-qm-monte-carlo-is)), the low-discrepancy points that change the rate, and the discretisation and multilevel methods that decide what a path costs.

## 26.1 Estimators and generators

**Definition 26.1 (Monte Carlo method, variance reduction).**

The *Monte Carlo method* estimates $\mu = \E[f(X)]$ by the average $\hat\mu_n = \frac1n\sum_{i=1}^nf(X_i)$ of independent draws. *Variance reduction* is any change of [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) that keeps $\E[\hat\mu_n] = \mu$ (or makes the bias negligible) and lowers its variance at equal cost.

**Proposition 26.2 (Error and cost of the Monte Carlo estimator).**

If $\sigma^2 = \Var(f(X)) < \infty$, then $\sqrt n(\hat\mu_n - \mu) \Rightarrow \mathcal N(0, \sigma^2)$, so $\hat\mu_n \pm 1.96\,s_n/\sqrt n$ is an asymptotic 95% interval. If one draw costs $c$, the cost of reaching [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) $\epsilon$ is $c\sigma^2/\epsilon^2$: of two [unbiased estimators](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator), the better is the one with the smaller product of variance and cost per draw.

**Proof.** The first statement is the central limit theorem for independent identically distributed draws, with Slutsky’s lemma for the estimated $s_n$. The second: $\sigma/\sqrt n = \epsilon$ needs $n =
\sigma^2/\epsilon^2$ draws. ∎

The rate $n^{-1/2}$ does not depend on the dimension, which is why Monte Carlo prices 64-date averages and 30-asset baskets that no grid can reach; the constant $\sigma^2c$ is where the work is. On the batch’s representative trade, an arithmetic-average call on 64 fixings of a [geometric Brownian motion](https://one-course.com/books/quant/4/en/chapter/4-stochastic-differential-equations#def-qm-stochastic-differential-equations-gbm) ($S_0 = K = 100$, $r = 3\%$, $\sigma = 25\%$, one year), $2^{14}$ paths give 6.495 with a [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) of 0.077; the price is 6.4788.

**Definition 26.3 (Pseudo-random and counter-based generators, inverse transform sampling).**

A *pseudo-random number generator* is a deterministic recurrence $s_{k+1} = F(s_k)$, $u_k = G(s_k)$ whose outputs pass statistical tests of independent uniformity; the Mersenne Twister (Matsumoto and Nishimura, 1998) and PCG (O’Neill, 2014; NumPy’s default) are examples. A *counter-based generator* computes $u = G_k(c)$ directly from a counter $c$ and a key $k$ with a bijection $G_k$, so that any draw can be produced without its predecessors; Philox (Salmon, Moraes, Dror and Shaw, 2011) is one. *Inverse transform sampling* turns a uniform $U$ into a draw $F^{-1}(U)$ of the distribution function $F$.

A [counter-based generator](#def-qm-monte-carlo-prng) is what makes a distributed batch reproducible. The engine gives path $j$ of trade $i$ the counter $(b, j, i, 0)$ for its $b$-th block of four numbers, so the path’s numbers do not depend on how the paths are split among machines, on the order in which they run, or on how many there are: rerunning one path to debug a trade reproduces it bit for bit. The kit’s Philox4x64-10 reproduces NumPy’s implementation block for block, in Python, C++20 and Rust (the C++ version is listed in [Section 26.5](#sec-26-5)). Normals come from the inverse transform with Acklam’s rational approximation of $\Phi^{-1}$ (relative error below $1.15 \times 10^{-9}$), which, unlike Box–Muller, maps one uniform to one normal: the property low-discrepancy points need.

## 26.2 Variance reduction and importance sampling

**Definition 26.4 (Antithetic variates, control variate, stratified sampling).**

*Antithetic variates* pair each draw $Z$ with $-Z$ (each uniform $U$ with $1 - U$) and average $\frac12(f(Z) + f(-Z))$. A *control variate* is a variable $C$ with known mean, used in $\hat\mu_\beta = \frac1n\sum_i(Y_i - \beta(C_i - \E C))$. *Stratified sampling* partitions the space of an input into strata of known probability $p_s$ and draws a fixed number $n_s = np_s$ in each.

**Proposition 26.5 (Optimal control-variate coefficient).**

$\Var(Y - \beta(C - \E C))$ is minimised at $\beta^* = \Cov(Y, C)/\Var(C)$, where it equals $(1 - \rho^2)\Var(Y)$ with $\rho = \mathrm{Corr}(Y, C)$.

**Proof.** $\Var(Y - \beta C) = \Var Y - 2\beta\Cov(Y, C) + \beta^2\Var C$ is a quadratic in $\beta$ minimised at $\beta^*$; substituting gives $\Var Y - \Cov(Y, C)^2/\Var C = (1 - \rho^2)\Var Y$. ∎

For the average-price call, the geometric average $G = (\prod_kS_{t_k})^{1/64}$ is lognormal and its call has a closed form (Kemna and Vorst, 1990): 6.1793. Its payoff has correlation 0.99941 with the arithmetic payoff’s, so $1 - \rho^2 = 0.0012$ and, with $\beta = 1.04$ estimated from the same paths, the variance falls by a factor of 853 at the same number of paths. Antithetic pairs do much less: the two payoffs of a pair have correlation $-0.43$, a factor of 1.8 at equal payoff evaluations. Stratification is the idea behind the next section: it removes the variance between strata, and a low-discrepancy point set stratifies every coordinate, and many of their products, at once.

**Definition 26.6 (Importance sampling).**

*Importance sampling* draws $X$ from a density $g$ instead of $f$ and averages $h(X)f(X)/g(X)$, the payoff times the likelihood ratio (a [Radon–Nikodym derivative](https://one-course.com/books/quant/4/en/chapter/1-probability-at-speed#def-qm-probability-at-speed-change)); it is unbiased whenever $g > 0$ where $hf \ne 0$.

**Proposition 26.7 (Zero-variance importance density).**

For $h \ge 0$ with $\mu = \E_f[h(X)] > 0$, the density $g^* = hf/\mu$ makes $h(X)f(X)/g^*(X) = \mu$ for every draw: the [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) has zero variance.

**Proof.** Substitute: $hf/g^* = hf\mu/(hf) = \mu$ wherever $h > 0$, and $g^*$ puts no mass where $h = 0$. ∎

The optimal density needs the answer, but it says where to look: sample in proportion to where the payoff lives. For a digital paying if $S_T > 250$ (probability $1.21 \times 10^{-4}$, the normal driver above $z^* = 3.67$), 100 000 plain draws produce six hits and a relative [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) of 20%. Shifting the driver’s mean to $\theta$ and weighting each draw by $e^{-\theta Z +
\theta^2/2}$ brings it to 0.64% at $\theta = 3.8$, a variance factor of about 1 000 ([Figure 26.1](#fig-qm-monte-carlo-is)); shifting too far makes it worse again, as the weights of the few draws that matter explode.

![Importance sampling for a digital paying if S_T > 250 (S_0 = 100, = 25\%, one year; probability 1.21 × 10-4): relative standard error of 100 000 draws against the mean shift of the normal driver. Data: the chapter’s tutorial, seeded.](https://one-course.com/images/onecourse/chapters/quant-4/qm-monte-carlo/fig-743d9181728c.svg)

***Figure 26.1.** [Importance sampling](#def-qm-monte-carlo-is) for a digital paying if $S_T > 250$ ($S_0 = 100$, $\sigma = 25\%$, one year; probability $1.21 \times 10^{-4}$): relative [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) of 100 000 draws against the mean shift of the normal driver. Data: the chapter’s tutorial, seeded.*

## 26.3 Low-discrepancy sequences

**Definition 26.8 (Low-discrepancy sequence, Sobol sequence, quasi-Monte Carlo).**

The star discrepancy $D_n^*$ of points $x_1, \dots, x_n$ in $[0, 1)^d$ is the largest difference, over boxes $[0, a)$, between the fraction of points in the box and its volume. A *low-discrepancy sequence* has $D_n^* = O((\log n)^d/n)$. The *Sobol sequence* (Sobol, 1967) is a base-2 low-discrepancy sequence built coordinate by coordinate from primitive polynomials and direction numbers, by exclusive-or of binary digits. *Quasi-Monte Carlo* replaces the random draws of Monte Carlo by such points.

**Proposition 26.9 (Koksma–Hlawka inequality).**

For $f$ of bounded variation $V(f)$ in the sense of Hardy and Krause, $\bigl|\frac1n\sum_if(x_i) - \int_{[0,1)^d}f\bigr| \le V(f)\,D_n^*$.

The bound gives [quasi-Monte Carlo](#def-qm-monte-carlo-qmc) a rate close to $1/n$ for smooth integrands of moderate effective dimension, but it is a worst-case bound that cannot be estimated from the points. Randomisation restores an error estimate.

**Definition 26.10 (Randomised quasi-Monte Carlo, Brownian bridge construction).**

*Randomised quasi-Monte Carlo* randomises a low-discrepancy set so that each point is uniform while the set keeps its structure, and estimates the error from independent replicates. The kit uses Owen’s nested uniform scrambling (1995): each binary digit of each coordinate is flipped by a random bit attached to the node of the binary tree that its leading digits define. The *Brownian bridge construction* assigns the first coordinate to $W_T$, the next to $W_{T/2}$ given its ends, and so on down the levels, so that the first coordinates, where low-discrepancy points are best, carry most of the variance of the path.

[Figure 26.2](#fig-qm-monte-carlo-points) shows why: 256 scrambled Sobol points in two dimensions put exactly one point in each of 256 dyadic boxes of every shape, where independent draws leave gaps and clusters. On the average-price call with $2^{14}$ paths and 32 independent scrambles, scrambled Sobol points with the increments in time order reduce the variance by a factor of about 34; with the [Brownian bridge construction](#def-qm-monte-carlo-rqmc) by about 2 500; with the bridge and the [control variate](#def-qm-monte-carlo-vr) together by about 37 000, a [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) 190 times smaller at the same number of paths. The error also falls faster than $n^{-1/2}$: fitted slopes between $-0.7$ and $-0.8$ from $2^6$ to $2^{15}$ paths ([Figure 26.3](#fig-qm-monte-carlo-cost)).

![256 points in the unit square: independent Philox uniforms and the first two coordinates of Owen-scrambled Sobol points. Data: the chapter’s tutorial, seeded.](https://one-course.com/images/onecourse/chapters/quant-4/qm-monte-carlo/fig-bd0e31dd58f2.svg)

***Figure 26.2.** 256 points in the unit square: independent Philox uniforms and the first two coordinates of Owen-scrambled Sobol points. Data: the chapter’s tutorial, seeded.*

![Standard error of the average-price call (64 fixings) against the number of paths: plain Monte Carlo and the control variate (slope -1/2), and scrambled Sobol points (16 independent scrambles) with increments in time order, with the Brownian bridge construction, and with the bridge and the control variate. Data: the chapter’s tutorial, seeded.](https://one-course.com/images/onecourse/chapters/quant-4/qm-monte-carlo/fig-48f9999a56d7.svg)

***Figure 26.3.** [Standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) of the average-price call (64 fixings) against the number of paths: plain Monte Carlo and the [control variate](#def-qm-monte-carlo-vr) (slope $-\frac12$), and scrambled Sobol points (16 independent scrambles) with increments in time order, with the [Brownian bridge construction](#def-qm-monte-carlo-rqmc), and with the bridge and the [control variate](#def-qm-monte-carlo-vr). Data: the chapter’s tutorial, seeded.*

## 26.4 Discretisation and multilevel methods

**Definition 26.11 (Strong and weak orders of convergence, Milstein scheme).**

A scheme with step $h$ has *strong order of convergence* $\gamma$ if $\E|X_T^h - X_T| \le Ch^\gamma$, and *weak order of convergence* $\beta$ if $|\E g(X_T^h) - \E g(X_T)| \le Ch^\beta$ for smooth $g$. The *Milstein scheme* for $dX = a(X)\,dt + b(X)\,dW$ adds to the Euler–Maruyama step the term $\frac12b(X)b'(X)((\Delta W)^2 - h)$.

**Proposition 26.12 (Orders of Euler and Milstein).**

Under Lipschitz and growth conditions (and smoothness of $a$, $b$ for Milstein), the [Euler–Maruyama scheme](https://one-course.com/books/quant/4/en/chapter/4-stochastic-differential-equations#def-qm-stochastic-differential-equations-euler) has strong order $\frac12$ and weak order 1, and the [Milstein scheme](#def-qm-monte-carlo-order) strong order 1 (Kloeden and Platen, 1992).

**Idea.** The Euler step misses the double stochastic integral $\int\!\!\int dW\,dW = \frac12((\Delta W)^2 - h)$, of size $h$ with mean zero, times $bb'$; these errors add like a random walk over $1/h$ steps, $h\sqrt{1/h} = h^{1/2}$ pathwise, while their means cancel, leaving an $O(h)$ bias in law. Milstein keeps that term, and the next missing terms are of order $h^{3/2}$ per step. ∎

A price needs the weak order; a hedge simulated along the path, or a multilevel [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator), needs the strong one. On a [geometric Brownian motion](https://one-course.com/books/quant/4/en/chapter/4-stochastic-differential-equations#def-qm-stochastic-differential-equations-gbm) with $\sigma = 50\%$, driven by the same increments as the exact solution, the fitted strong slopes are 0.49 for Euler and 0.89 for Milstein, and the weak error of the mean is exactly $X_0|(1 + \mu h)^n - e^{\mu T}|$, slope 1.00 ([Figure 26.4](#fig-qm-monte-carlo-orders)).

![Strong error |X_Th - X_T| of the Euler and Milstein schemes (20 000 paths of a geometric Brownian motion, = 5\%, = 50\%, X_0 = 100, one year) and the exact weak error of the mean, against the time step. Data: the chapter’s tutorial, seeded.](https://one-course.com/images/onecourse/chapters/quant-4/qm-monte-carlo/fig-2514905cd724.svg)

***Figure 26.4.** Strong error $\E|X_T^h - X_T|$ of the Euler and [Milstein schemes](#def-qm-monte-carlo-order) (20 000 paths of a [geometric Brownian motion](https://one-course.com/books/quant/4/en/chapter/4-stochastic-differential-equations#def-qm-stochastic-differential-equations-gbm), $\mu = 5\%$, $\sigma = 50\%$, $X_0 = 100$, one year) and the exact weak error of the mean, against the time step. Data: the chapter’s tutorial, seeded.*

**Definition 26.13 (Multilevel Monte Carlo).**

*Multilevel Monte Carlo* (Giles, 2008) writes the finest-level expectation as a telescoping sum $\E P_L = \E P_0 + \sum_{l=1}^L\E[P_l - P_{l-1}]$, where $P_l$ is the payoff computed with step $h_l = 2^{-l}T$, and estimates each term independently, with $P_l$ and $P_{l-1}$ in each correction computed from the same Brownian path.

**Proposition 26.14 (Complexity of multilevel Monte Carlo).**

If the bias of $P_l$ is $O(2^{-\alpha l})$, the variance of the correction $V_l = O(2^{-\beta l})$ and its cost $C_l = O(2^{\gamma l})$, with $\alpha \ge \frac12\min(\beta, \gamma)$, then a mean square error $\epsilon^2$ costs $O(\epsilon^{-2})$ if $\beta > \gamma$, $O(\epsilon^{-2}(\log\epsilon)^2)$ if $\beta = \gamma$ and $O(\epsilon^{-2 - (\gamma - \beta)/\alpha})$ if $\beta < \gamma$; standard Monte Carlo costs $O(\epsilon^{-2 - \gamma/\alpha})$.

**Idea.** Minimising the cost $\sum_lN_lC_l$ subject to the variance constraint $\sum_lV_l/N_l = \epsilon^2/2$ gives $N_l \propto\sqrt{V_l/C_l}$ and a total cost

$$
2\epsilon^{-2}\Bigl(\sum_l\sqrt{V_lC_l}\Bigr)^2;
$$

the sum is dominated by the coarsest levels if $\beta > \gamma$, by the finest if $\beta < \gamma$, and grows like $L \sim \log\epsilon^{-1}$ if they are equal. The finest level $L$ is set by the bias, $2^{-\alpha L} \sim \epsilon$. ∎

The strong order sets $\beta$: Milstein on a Lipschitz payoff gives $\beta = 2 > \gamma = 1$, and the corrections’ variances in the tutorial fall by a factor of about four per level (0.23, 0.060, 0.016, 0.0045, 0.0012, 0.00030, 0.000079). For a European call ($\alpha = 1$), the driver listed in [Section 26.5](#sec-26-5) reaches a root mean square error of 0.005 with 7 levels at a cost of $3.2 \times 10^7$ time steps, where standard Monte Carlo at the same finest step would need $3.0 \times 10^9$: a factor of 94, which grows as $\epsilon$ falls (8, 20, 45, 94 for $\epsilon$ from 0.04 to 0.005; [Figure 26.5](#fig-qm-monte-carlo-mlmc)).

![Cost of pricing a European call (S_0 = K = 100, r = 3\%, = 25\%, one year) to a target root mean square error: multilevel Monte Carlo with Milstein corrections against standard Monte Carlo at the same finest step. Data: the chapter’s tutorial, seeded.](https://one-course.com/images/onecourse/chapters/quant-4/qm-monte-carlo/fig-1f976fd2f968.svg)

***Figure 26.5.** Cost of pricing a European call ($S_0 = K = 100$, $r = 3\%$, $\sigma = 25\%$, one year) to a target root mean square error: [multilevel Monte Carlo](#def-qm-monte-carlo-mlmc) with Milstein corrections against standard Monte Carlo at the same finest step. Data: the chapter’s tutorial, seeded.*

## 26.5 Tutorial: sixteen times fewer machines

**Goal.** Value the batch’s representative trade five ways and measure each [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator)’s variance at equal paths. **End state:** the factors 1.8, 853, 34, 2 500 and 37 000 and [Figure 26.3](#fig-qm-monte-carlo-cost).

1. **Counter-based uniforms.** Philox4x64-10 in C++20, reproducing NumPy’s `Philox` block for block. `inline Block philox4x64 (Block c, std::array<std::uint64_t , 2 > k) { constexpr std::uint64_t M0 = 0xD2E7470EE14C6C93ULL , M1 = 0xCA5A826395121157ULL ; constexpr std::uint64_t W0 = 0x9E3779B97F4A7C15ULL , W1 = 0xBB67AE8584CAA73BULL ; for (int round = 0 ; round < 10 ; ++round) { const unsigned __int128 p0 = static_cast <unsigned __int128>(M0) * c[0 ]; const unsigned __int128 p1 = static_cast <unsigned __int128>(M1) * c[2 ]; c = {static_cast <std::uint64_t >(p1 >> 64 ) ^ c[1 ] ^ k[0 ], static_cast <std::uint64_t >(p1), static_cast <std::uint64_t >(p0 >> 64 ) ^ c[3 ] ^ k[1 ], static_cast <std::uint64_t >(p0)}; k = {k[0 ] + W0, k[1 ] + W1}; } return c; }` **Listing 26.1.** Philox4x64-10: ten rounds of two wide multiplications (C++20). code/firm/mcengine/cpp/firm_mcengine.hpp
2. **Scrambled Sobol points.** Gray-code Sobol points from Joe and Kuo’s direction numbers, and Owen’s scrambling by hashed tree nodes. `def sobol_points (n: int , dim: int ) -> np.ndarray: """The first n Sobol points as 32-bit integers, shape (n, dim), in Gray-code order (point 0 is 0).""" V = sobol_directions(dim) i = np.arange(n, dtype=np.uint64) g = i ^ (i >> np.uint64(1 )) X = np.zeros((n, dim), dtype=np.uint64) for b in range (max (1 , int (n - 1 ).bit_length())): on = ((g >> np.uint64(b)) & np.uint64(1 )).astype(bool ) X[on] ^= V[:, b] return X def _mix64 (z): z = (z ^ (z >> np.uint64(30 ))) * np.uint64(0xBF58476D1CE4E5B9 ) z = (z ^ (z >> np.uint64(27 ))) * np.uint64(0x94D049BB133111EB ) return z ^ (z >> np.uint64(31 )) def owen_scramble (X: np.ndarray, seed: int ) -> np.ndarray: """Nested uniform scrambling: digit l of coordinate d is flipped by a random bit attached to the node (d, l, first l digits) of the binary tree, drawn by hashing. Returns uniforms (y + 1/2) 2^-32.""" Y = X.copy() for d in range (X.shape[1 ]): x = X[:, d] for lev in range (32 ): prefix = x >> np.uint64(32 - lev) if lev else np.zeros_like(x) node = np.uint64(d << 40 ) ^ np.uint64(lev << 32 ) ^ prefix bit = _mix64(np.uint64(seed & MASK) ^ _mix64(node)) >> np.uint64(63 ) Y[:, d] ^= bit << np.uint64(31 - lev) return (Y.astype(float ) + 0.5 ) * 2.0 **-32` **Listing 26.2.** Sobol points and nested uniform scrambling (Python). code/firm/mcengine/firm_mcengine.py
3. **[Estimators](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator).** Run `table()` in `qm_mc.py` : plain, antithetic, [control variate](#def-qm-monte-carlo-vr) , and scrambled Sobol with and without the bridge (32 scrambles of $2^{14}$ points); then `error_vs_cost()` .
4. **Multilevel.** Run `mlmc_study()`; the driver is below. `def mlmc (sampler, eps: float , L0: int = 2 , N0: int = 10_000 , alpha: float = 1.0 , max_level: int = 12 ) -> dict : """Giles's multilevel Monte Carlo. sampler(l, n) returns n samples of Y_l = P_l - P_{l-1} (P_0 at l = 0) and the cost of one sample; the sample sizes minimise cost for variance eps^2/2, and levels are added until the estimated bias |E Y_L| / (2^alpha - 1) is below eps / sqrt(2).""" L = L0 s1, s2, N, C = [0.0 ] * (L + 1 ), [0.0 ] * (L + 1 ), [0 ] * (L + 1 ), [0.0 ] * (L + 1 ) dN = [N0] * (L + 1 ) while True : for lev in range (L + 1 ): if dN[lev] > 0 : y, c = sampler(lev, dN[lev]) s1[lev] += float (y.sum()) s2[lev] += float ((y * y).sum()) N[lev] += dN[lev] C[lev] = c m = [s1[k] / N[k] for k in range (L + 1 )] V = [max (s2[k] / N[k] - m[k] ** 2 , 1e-300 ) for k in range (L + 1 )] tot = sum (math.sqrt(V[k] * C[k]) for k in range (L + 1 )) Nopt = [math.ceil(2 / eps**2 * math.sqrt(V[k] / C[k]) * tot) for k in range (L + 1 )] dN = [max (0 , Nopt[k] - N[k]) for k in range (L + 1 )] if any (dN): continue bias = max (abs (m[L]), abs (m[L - 1 ]) / 2 **alpha) / (2 **alpha - 1 ) if bias < eps / math.sqrt(2 ) or L == max_level: return {" est " : sum (m), " N " : N, " V " : V, " mean " : m, " C " : C, " L " : L, " cost " : sum (N[k] * C[k] for k in range (L + 1 ))} L += 1 s1.append(0.0 ), s2.append(0.0 ), N.append(0 ), C.append(0.0 ) dN = [0 ] * L + [N0]` **Listing 26.3.** Giles’s multilevel Monte Carlo driver (Python). code/firm/mcengine/firm_mcengine.py

**What to change next.** Put the strike at 130 and watch the [control variate](#def-qm-monte-carlo-vr)’s factor fall with the correlation; replace the bridge by a principal-component construction; price a barrier with the bridge crossing probability of chapter 2 inside the multilevel driver.

## 26.6 Build: the Monte Carlo engine, stage three

**Purpose.** Make each path worth more and every path reproducible: the generators, point sets and [estimators](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) that the firm’s pricing and risk batches call.

**Interface.** `philox4x64`, `philox_uniforms(n, dim, seed, stream, first_path)`; `sobol_points`, `owen_scramble`; `norm_ppf`; `bridge_from_normals`; `asian_payoffs`, `geometric_asian_price`, `control_variate`; `mlmc(sampler, eps)`. C++20 and Rust twins of `philox4x64`, `sobol_point` (eight dimensions) and `owen_scramble`.

**Rules.** Counters, not sequential state, identify paths; every estimate is reported with its [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) (from replicates for [randomised quasi-Monte Carlo](#def-qm-monte-carlo-rqmc), never from within one point set); unscrambled Sobol points are never used for error bars; direction numbers carry their licence.

**Acceptance tests.** `code/firm/mcengine/{tests,cpp,rust}`: Philox against NumPy’s; splitting paths reproduces them; Joe and Kuo’s published first points; identical scrambled integers in the three languages; the net property after scrambling; $\Phi^{-1}$ against the standard library to $10^{-8}$; bridge covariance; the geometric closed form against simulation and against Black–Scholes for one fixing; the multilevel driver on a sampler with a known answer.

**Stretch.** Joe and Kuo’s full 21 201 dimensions; a principal-component path construction; conditional Monte Carlo for barriers; adjoint Greeks along the paths (chapter 28).

Sources and further reading

- N. Metropolis and S. Ulam, “The Monte Carlo method”, *Journal of the American Statistical Association* 44, 1949; P. P. Boyle, “Options: a Monte Carlo approach”, *Journal of Financial Economics* 4, 1977.
- P. Glasserman, *Monte Carlo Methods in Financial Engineering* , Springer, 2003.
- J. K. Salmon, M. A. Moraes, R. O. Dror and D. E. Shaw, “Parallel random numbers: as easy as 1, 2, 3”, SC11, 2011; M. Matsumoto and T. Nishimura, *ACM TOMACS* 8, 1998; M. E. O’Neill, PCG, Harvey Mudd College report HMC-CS-2014-0905, 2014.
- I. M. Sobol’, *USSR Computational Mathematics and Mathematical Physics* 7, 1967; S. Joe and F. Y. Kuo, *SIAM Journal on Scientific Computing* 30, 2008 (direction numbers `new-joe-kuo-6.21201` ); A. B. Owen, “Randomly permuted $(t, m, s)$ -nets and $(t, s)$ -sequences”, 1995; B. Moskowitz and R. E. Caflisch, *Mathematical and Computer Modelling* 23, 1996.
- A. G. Z. Kemna and A. C. F. Vorst, *Journal of Banking and Finance* 14, 1990.
- M. B. Giles, “Multilevel Monte Carlo path simulation”, *Operations Research* 56, 2008; P. E. Kloeden and E. Platen, *Numerical Solution of Stochastic Differential Equations* , Springer, 1992; G. N. Milstein, *Theory of Probability and its Applications* 19, 1975.
- P. J. Acklam, “An algorithm for computing the inverse normal cumulative distribution function” (web note, 2004).

## 26.7 Exercises

**Exercise 26.1 ★.**

A batch has [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) 80 000 dollars with 100 000 paths. How many paths reach 20 000, and what variance-reduction factor does the same at 100 000 paths?

**Solution of Exercise 26.1.**

Error scales as $n^{-1/2}$: a factor of four needs $16 \times 100\,000 = 1.6$ million paths, or a variance-reduction factor of 16 at 100 000 paths.

**Exercise 26.2 ★.**

A [control variate](#def-qm-monte-carlo-vr) has correlation 0.9 with the payoff. What is the variance factor? And with 0.99?

**Solution of Exercise 26.2.**

$1/(1 - \rho^2)$: 5.3 at $\rho = 0.9$ and 50.3 at $\rho = 0.99$. High correlations are needed for large factors.

**Exercise 26.3 ★.**

Why do [antithetic variates](#def-qm-monte-carlo-vr) help a monotone payoff and not an even function of $Z$ such as $Z^2$?

**Solution of Exercise 26.3.**

$\Var(\frac12(f(Z) + f(-Z))) = \frac12\Var f(Z)(1 + \mathrm{Corr}(f(Z), f(-Z)))$ per pair, i.e. per two evaluations; it beats plain sampling when the correlation is negative, which a monotone $f$ guarantees. For an even function $f(Z) = f(-Z)$, the correlation is 1 and the pair is worth one draw at the cost of two.

**Exercise 26.4 ★★.**

Show that proportional stratification never increases the variance: with strata $s$ of probability $p_s$ and conditional variances $\sigma_s^2$, the stratified variance is $\frac1n\sum_sp_s\sigma_s^2 \le
\sigma^2/n$.

**Solution of Exercise 26.4.**

By the law of total variance, $\sigma^2 = \sum_sp_s\sigma_s^2 + \sum_sp_s(\mu_s - \mu)^2$. The stratified [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) $\sum_sp_s\bar Y_s$ with $n_s = np_s$ has variance $\sum_sp_s^2\sigma_s^2/n_s =
\frac1n\sum_sp_s\sigma_s^2$, which drops the nonnegative between-strata term.

**Exercise 26.5 ★★.**

For $P(Z > z^*)$ with $Z$ standard normal, compute the second moment of the importance-sampling [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) with the shift $\theta$, and show that it is minimised near $\theta = z^*$ for large $z^*$.

**Solution of Exercise 26.5.**

The [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) is $\mathbf 1\{Z > z^*\}e^{-\theta Z + \theta^2/2}$ with $Z \sim \mathcal N(\theta, 1)$; its second moment is

$$
\E_f\bigl[\mathbf 1\{Z > z^*\}e^{-\theta Z + \theta^2/2}\bigr] = e^{\theta^2}\bar\Phi(z^* + \theta).
$$

The derivative of its logarithm is $2\theta - \varphi(z^* + \theta)/\bar\Phi(z^* + \theta) \approx 2\theta - (z^* + \theta)$ for large arguments (Mills’s ratio), which vanishes at $\theta \approx z^*$. At $z^* = 3.67$ the exact minimiser is 3.80, where the simulation’s minimum is.

**Exercise 26.6 ★★.**

Derive the Milstein step for [geometric Brownian motion](https://one-course.com/books/quant/4/en/chapter/4-stochastic-differential-equations#def-qm-stochastic-differential-equations-gbm) and show that it is the second-order expansion of the exact step $X e^{(\mu - \sigma^2/2)h + \sigma\Delta W}$ in $\Delta W$.

**Solution of Exercise 26.6.**

With $a(x) = \mu x$ and $b(x) = \sigma x$, $bb' = \sigma^2x$: $X_{k+1} = X_k(1 + \mu h + \sigma\Delta W + \frac12\sigma^2(\Delta W^2 - h))$. Expanding the exact step, $e^{(\mu - \sigma^2/2)h + \sigma\Delta W} = 1 + \sigma\Delta W + (\mu - \frac12\sigma^2)h + \frac12\sigma^2\Delta W^2 + O(h^{3/2})$, the same terms.

**Exercise 26.7 ★★★.**

*Coding.* Compute the variance factor of the [control variate](#def-qm-monte-carlo-vr) at strikes 80, 100 and 130, and explain the trend.

**Solution of Exercise 26.7.**

With $2^{14}$ paths: correlation 0.99953, 0.99941 and 0.99425, variance factors 1 071, 853 and 87 at strikes 80, 100 and 130. Far out of the money the geometric average (always below the arithmetic one) finishes in the money less often than the arithmetic one, and the two payoffs decouple near the strike: the control is worth less where both are rare.

**Exercise 26.8 ★★★.**

*Find the flaw.* “We moved the batch to Sobol points, 100 000 per trade, and report the [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) $s/\sqrt n$ from the points as before; it fell by a factor of ten.”

**Solution of Exercise 26.8.**

Sobol points are not independent, so $s/\sqrt n$ is not the [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) of their average: it measures the spread of the payoffs, not of the estimate, and says nothing about the error, which may be smaller or, with a poor path construction, larger. The error must come from independent scrambles (replicates), and unscrambled points give no error estimate at all.

## 26.8 Problem: Sixteen Times Fewer Machines

**Problem 26.1.**

Weekend problem — a factor of four in error, without sixteen times the paths

The representative trade is an arithmetic-average call on 64 fixings ($S_0 = K = 100$, $r = 3\%$, $\sigma = 25\%$, one year), valued with $2^{14}$ paths; the batch needs a factor of four in [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator).

**Part I — The plain [estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator).**

1. What are the plain estimate, its [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) and the reference price?
2. What variance factor does the target require?
3. What does the cost-variance product say about comparing [estimators](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) ?
4. Why does the [counter-based generator](#def-qm-monte-carlo-prng) make the batch reproducible across machines?
5. Why is the inverse transform, rather than Box–Muller, used to feed low-discrepancy points?

**Part II — [Variance reduction](#def-qm-monte-carlo-mc).**

6. What factor do antithetic pairs give, and why so little?
7. What is the [control variate](#def-qm-monte-carlo-vr) , its closed-form value, its correlation and its factor?
8. What is the optimal coefficient, and what was estimated?
9. Where does [importance sampling](#def-qm-monte-carlo-is) help, and by how much on the deep out-of-the-money digital?
10. Why does too large a shift make it worse?

**Part III — Points and paths.**

11. What factor do scrambled Sobol points give with the increments in time order, and with the bridge?
12. Why does the bridge matter so much?
13. What do the two together give, and what is the error’s slope in the number of paths?
14. How is the error of [randomised quasi-Monte Carlo](#def-qm-monte-carlo-rqmc) estimated?
15. What do the strong and weak orders measure, and what are their fitted values here?

**Part IV — Judgement.**

16. What does [multilevel Monte Carlo](#def-qm-monte-carlo-mlmc) save on the European call at $\epsilon = 0.005$ , and why?
17. Would the factor of 37 000 hold across the batch?
18. What must the batch log so that one trade’s paths can be regenerated?
19. State the *named result* : the variance-reduction factor of the [control variate](#def-qm-monte-carlo-vr) plus [randomised quasi-Monte Carlo](#def-qm-monte-carlo-rqmc) on the representative trade, against the 16 required.
20. In one sentence: where does the work of a Monte Carlo engine go?

**Solution of Problem 26.1.**

**1.** 6.495 with [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) 0.077 from $2^{14}$ paths; reference 6.4788 ([control variate](#def-qm-monte-carlo-vr) plus scrambled Sobol, 32 replicates of $2^{14}$, [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) 0.00007). **2.** 16. **3.** At equal accuracy the cost is $c\sigma^2/\epsilon^2$, so [estimators](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) are compared by variance times cost per draw; all the factors here are at equal payoff evaluations. **4.** Each path’s numbers are a function of (seed, trade, path, block), not of a sequential state, so any split of the paths reproduces them. **5.** It maps one uniform coordinate to one normal monotonically, so the structure of the point set survives; Box–Muller mixes two coordinates through a cosine. **6.** 1.8: the two payoffs of a pair have correlation $-0.43$; the payoff is monotone in the path but convex, and its positive part discards the down paths. **7.** The geometric-average call, 6.1793 in closed form; correlation 0.99941; factor 853. **8.** $\beta^* = \Cov(Y, C)/\Var C$; estimated from the same paths as 1.04 (the bias this introduces is of order $1/n$). **9.** Where the payoff lives in a small region of the driver’s space: at $\theta = 3.8$ the relative [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) falls from 20% to 0.64%, a factor of about 1 000 in variance. **10.** The likelihood ratio $e^{-\theta Z + \theta^2/2}$ becomes very variable: the second moment $e^{\theta^2}\bar\Phi(z^* + \theta)$ grows again beyond the optimum. **11.** About 34 with the increments in time order and about 2 500 with the bridge. **12.** Low-discrepancy points are most uniform in their first coordinates; the bridge puts $W_T$ and the coarse shape of the path there, the directions that carry most of the payoff’s variance. **13.** About 37 000 ([standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) 190 times smaller); fitted slopes between $-0.7$ and $-0.8$ against $-0.5$ for Monte Carlo. **14.** From the spread of independent scrambles: $R$ replicates give $R$ independent unbiased estimates, whose sample standard deviation divided by $\sqrt R$ is the [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator). **15.** Strong: the pathwise error, slopes 0.49 (Euler) and 0.89 (Milstein); weak: the error in law, slope 1.00 for the mean. **16.** Cost $3.2 \times 10^7$ time steps against $3.0 \times 10^9$, a factor of 94: most samples are drawn on coarse, cheap levels, and the fine corrections have small variance because Milstein’s strong order 1 makes $V_l$ fall by four per level. **17.** No: it relies on a smooth payoff with a closed-form control and a low effective dimension. Barriers, digitals and early exercise break the smoothness, and most trades have no control as good; each trade class must be measured. **18.** Seed, trade and path counters, generator and version, scrambling seeds, the path construction, and the point-set size. **19.** Named result: *sixteen times fewer machines*: on the representative trade, the [control variate](#def-qm-monte-carlo-vr) alone reduces the variance by 853 and with scrambled Sobol points and the [Brownian bridge](https://one-course.com/books/quant/4/en/chapter/2-brownian-motion#def-qm-brownian-motion-bridge) by about 37 000, against the 16 required; the target is met at today’s number of paths, with the generation overhead to be measured in the production engine. **20.** Into the constant in front of $n^{-1/2}$: what each path is worth, and what it costs.

## 26.9 Interview questions

**Interview question 26.1 ★ researcher, developer.**

How does the Monte Carlo error scale with the number of paths and with the dimension?

**Solution of Interview question 26.1.**

The [standard error](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator) is $\sigma/\sqrt n$ whatever the dimension: the dimension enters only through $\sigma$ and the cost of a path, which is why Monte Carlo wins in high dimension and loses to grids in one or two.

*What the interviewer is looking for: $n^{-1/2}$, dimension-free, the constant.*

**Interview question 26.2 ★★ researcher.**

Name three variance-reduction techniques and when each works.

**Solution of Interview question 26.2.**

[Antithetic variates](#def-qm-monte-carlo-vr) (monotone payoffs, cheap, small gains); [control variates](#def-qm-monte-carlo-vr) (a correlated quantity with known mean: the geometric Asian, the underlying, a Black–Scholes price); [importance sampling](#def-qm-monte-carlo-is) (rare events: shift the drift towards the region that matters); stratification or [quasi-Monte Carlo](#def-qm-monte-carlo-qmc) (smooth payoffs in a few important directions).

*What the interviewer is looking for: When each works, not just names.*

**Interview question 26.3 ★★ developer.**

How do you make a distributed Monte Carlo batch reproducible regardless of the number of machines?

**Solution of Interview question 26.3.**

Use a [counter-based generator](#def-qm-monte-carlo-prng) (Philox, Threefry) keyed by the seed, with the counter built from trade, path and block indices, so that numbers are a function of their coordinates, not of the order of generation; reduce partial sums in a fixed order.

*What the interviewer is looking for: Counters over sequential streams, and a fixed reduction order.*

**Interview question 26.4 ★★ researcher.**

Why might [quasi-Monte Carlo](#def-qm-monte-carlo-qmc) do badly on a 250-step path, and what fixes it?

**Solution of Interview question 26.4.**

Its uniformity is good only in the first coordinates, and in time order the path’s variance is spread evenly over 250 of them. Fix it with a [Brownian bridge](https://one-course.com/books/quant/4/en/chapter/2-brownian-motion#def-qm-brownian-motion-bridge) or principal-component construction that puts the important directions first, and scramble to get error estimates.

*What the interviewer is looking for: Effective dimension and path construction.*

**Interview question 26.5 ★★ researcher, risk.**

How would you estimate the probability of a loss beyond a very high threshold by simulation?

**Solution of Interview question 26.5.**

[Importance sampling](#def-qm-monte-carlo-is): tilt the distribution (shift the mean of the drivers, or exponentially tilt the loss) so that the loss beyond the threshold is common, and weight by the likelihood ratio; choose the tilt near the most likely point of the rare region; check the weights’ spread. Splitting or conditional Monte Carlo are alternatives.

*What the interviewer is looking for: A [change of measure](https://one-course.com/books/quant/4/en/chapter/1-probability-at-speed#def-qm-probability-at-speed-change) with the likelihood ratio, and a sensible choice of tilt.*

**Interview question 26.6 ★★★ researcher.**

Explain [multilevel Monte Carlo](#def-qm-monte-carlo-mlmc) and why the strong order of the scheme matters for it.

**Solution of Interview question 26.6.**

Write the finest-level expectation as a telescoping sum of corrections between successive step sizes, simulate the fine and coarse paths of each correction from the same Brownian path, and allocate samples in proportion to $\sqrt{V_l/C_l}$. The correction’s variance $V_l$ is governed by the strong error, so a higher strong order makes the fine levels need very few samples; with $V_l$ falling faster than the cost rises, the total cost is $O(\epsilon^{-2})$, as for an [unbiased estimator](https://one-course.com/books/quant/4/en/chapter/11-estimation#def-qm-estimation-estimator).

*What the interviewer is looking for: Telescoping sum, coupling, and the variance of corrections from the strong order.*