Monte Carlo Integration & Precision

Let $\theta$ denote the integral of a real-valued function $h(\mathbf{x})$ over a domain $\mathcal{D}$ with respect to a continuous probability density function $f(\mathbf{x})$:

If the original integral is an unweighted Lebesgue integral over a bounded hypercube $\mathcal{D} = [a_1, b_1] \times \dots \times [a_d, b_d]$ with total hypervolume $V = \prod_{j=1}^d (b_j - a_j)$:

where $\mathbf{U} \sim \mathcal{U}(\mathcal{D})$ is uniformly distributed across the domain.

1. The Monte Carlo Estimator

Let $\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_n$ be an independent and identically distributed (i.i.d.) random sample drawn from density $f(\mathbf{x})$. The Crude Monte Carlo Estimator of $\theta$ is the sample mean of the evaluated function values:

2. Theoretical Statistical Properties

1. Unbiasedness:

The estimator is strictly unbiased for any finite sample size $n \ge 1$.

2. Strong Consistency (Convergence Almost Surely):

Assuming the function is absolutely integrable ($\mathbb{E}[|h(\mathbf{X})|] < \infty$), by the Kolmogorov Strong Law of Large Numbers:

3. Finite-Sample Variance:

Assuming square integrability ($\sigma^2 = \mathrm{Var}_f(h(\mathbf{X})) = \int (h(\mathbf{x}) - \theta)^2 f(\mathbf{x})\,dx < \infty$):

The standard error of the Monte Carlo estimator decays at the universal rate:

By the Lindeberg-Lévy Central Limit Theorem, as the simulation sample size $n \to \infty$, the standardized estimation error converges in distribution to the standard Gaussian distribution:

In practice, the true population variance $\sigma^2$ is unknown and must be estimated consistently from the simulated trajectory using the sample variance:

By Slutsky's theorem:

Asymptotic Confidence Intervals

An asymptotic $(1 - \alpha) \times 100\%$ confidence interval for the true integral $\theta$ is:

where $z_{1 - \alpha/2} = \Phi^{-1}(1 - \alpha/2)$ is the upper standard normal critical quantile (e.g., $z_{0.975} \approx 1.960$ for a $95\%$ confidence interval).

To appreciate why Monte Carlo integration is irreplaceable in modern science, consider the mathematical convergence rates of classical deterministic numerical integration algorithms in $d$-dimensional space. ### 1. Deterministic Numerical Quadrature Error Rates To evaluate an integral over the unit hypercube $[0, 1]^d$ using a grid-based quadrature rule: - Let $k$ denote the number of grid subdivisions along each individual coordinate axis. - The total number of evaluation nodes across the hypercube is: $$n = k^d \implies k = n^{1/d}$$ For standard classical quadrature rules: 1. **Riemann Sum / Midpoint Rule**: The error scales as $\mathcal{O}(1 / k) = \mathcal{O}(n^{-1/d})$. 2. **Trapezoidal Rule**: For twice-differentiable integrands, the error scales as $\mathcal{O}(1 / k^2) = \mathcal{O}(n^{-2/d})$. 3. **Simpson's Rule**: For four-times differentiable integrands, the error scales as $\mathcal{O}(1 / k^4) = \mathcal{O}(n^{-4/d})$. ### 2. The Dimensionality Transition Point Notice the critical mathematical contrast: - **Deterministic Quadrature**: Error $= \mathcal{O}(n^{-r / d})$, where $r$ is the algebraic order of the quadrature method. The convergence rate **degrades exponentially as dimension $d$ grows**! - **Monte Carlo Integration**: Error $= \mathcal{O}(n^{-1/2})$, completely **independent of the dimension $d$**! We equate the convergence exponents: $$\frac{r}{d} = \frac{1}{2} \implies d = 2r$$ - For the Trapezoidal rule ($r = 2$): - In $d = 1$: Quadrature achieves $\mathcal{O}(n^{-2})$, far faster than Monte Carlo $\mathcal{O}(n^{-1/2})$. - In $d = 4$: Quadrature achieves $\mathcal{O}(n^{-2/4}) = \mathcal{O}(n^{-1/2})$, breaking even with Monte Carlo. - In $d = 10$: Quadrature achieves $\mathcal{O}(n^{-0.20})$, while Monte Carlo maintains $\mathcal{O}(n^{-0.50})$. - In $d = 100$: Quadrature achieves $\mathcal{O}(n^{-0.02})$, requiring $10^{100}$ evaluations for two digits of precision, whereas Monte Carlo requires $10^4$ evaluations for the exact same accuracy! Monte Carlo integration is the only computationally tractable method for multidimensional problems when $d \ge 8$.

A central question in experimental simulation design is: how many Monte Carlo samples $n$ must be drawn to guarantee that the estimate $\hat{\theta}_n$ lies within a prescribed margin of error $\epsilon$ with high probability?

1. Absolute Error Bound

Let the researcher demand that the absolute error does not exceed $\epsilon > 0$ with confidence $1 - \alpha$:

Under the asymptotic normal approximation:

Squaring both sides yields the minimum sample size:

Key Observation: To gain one additional decimal digit of precision (reducing $\epsilon$ by a factor of 10), the required sample size $n$ must increase by a factor of 100 ($10^2$).

2. Relative Error Bound

When the magnitude of $\theta$ is unknown or very small (e.g., simulating rare failure events in structural safety where $\theta = 10^{-6}$), an absolute error bound is meaningless. The researcher specifies a relative precision error $\epsilon_r \in (0, 1)$:

The relative standard error requirement is:

The ratio $CV = \sigma / \theta$ is the Coefficient of Variation of the summand $h(\mathbf{X})$.

For a rare-event probability $\theta = P(X \in \mathcal{E})$ estimated via indicator variables $h(\mathbf{X}) = I(\mathbf{X} \in \mathcal{E})$:

The required sample size is:

If $\theta = 10^{-5}$ and the researcher desires a $10\%$ relative error ($\epsilon_r = 0.10$) at $95\%$ confidence ($z \approx 1.96$):

Crude Monte Carlo becomes computationally inefficient for rare-event estimation, necessitating advanced variance reduction techniques (such as Importance Sampling).

While Crude Monte Carlo integration achieves an asymptotic standard error of $O(n^{-1/2})$ regardless of dimension $d$, this convergence rate is agonizingly slow: achieving four decimal places of accuracy requires $n \approx 10^8$ sample draws.

Quasi-Monte Carlo (QMC) replaces pseudo-random independent draws with deterministic, highly uniform points known as Low-Discrepancy Sequences. By deliberately preventing points from clustering or leaving empty voids in the unit hypercube $[0, 1]^d$, QMC achieves an asymptotic convergence rate approaching:

For moderate dimensions $d$, this represents an almost linear convergence rate $O(n^{-1})$, vastly outperforming pseudo-random Monte Carlo.

1. Discrepancy and the Koksma-Hlawka Inequality

Let $P_n = \{\mathbf{x}_1, \dots, \mathbf{x}_n\}$ be a deterministic point set in the unit hypercube $I^d = [0, 1]^d$. For any sub-box anchored at the origin $[\mathbf{0}, \mathbf{y}) = \prod_{j=1}^d [0, y_j)$, the empirical count of points falling in the box is:

The local discrepancy function measures the deviation between the empirical fraction of points and the theoretical volume $\lambda([\mathbf{0}, \mathbf{y})) = \prod_{j=1}^d y_j$:

The Star Discrepancy $D_n^*(P_n)$ is the supremum norm of this deviation over all origin-anchored boxes:

Theorem (The Koksma-Hlawka Inequality): If $f: [0, 1]^d \to \mathbb{R}$ has bounded variation in the sense of Hardy and Krause, denoted $V_{\text{HK}}(f) < \infty$, then the deterministic integration error is bounded by:

This fundamental theorem separates the integration error into two independent components:

1. $V_{\text{HK}}(f)$: depends purely on the smoothness and regularity of the integrand $f$.
2. $D_n^*(P_n)$: depends purely on the geometric uniformity of the point set $P_n$.

2. Constructions of Low-Discrepancy Sequences

A sequence $\{\mathbf{x}_n\}$ is called a low-discrepancy sequence if its first $n$ elements satisfy $D_n^*(P_n) \le C_d \frac{(\log n)^d}{n}$ for all $n$. Three major families dominate computational practice:

1. Halton Sequences: Generalize the 1D van der Corput sequence to $d$ dimensions. Let $p_1, \dots, p_d$ be the first $d$ prime numbers. The $i$-th point is:

where $\phi_p(i) = \sum_{k=0}^M a_k(i) p^{-(k+1)}$ is the radical inverse function expanding integer $i = \sum_{k=0}^M a_k(i) p^k$ in base $p$.

2. Sobol' Sequences: Constructed in base 2 using primitive polynomials over the Galois field $\mathbb{F}_2$. A set of direction integers $v_j$ defines the bitwise XOR transformations generating coordinates. Sobol' sequences satisfy the $(t, s)$-sequence criteria, guaranteeing optimal uniformity across dyadic elementary intervals.

3. Faure and Niederreiter Sequences: Exploit algebraic geometry and polynomial rings over general finite fields $\mathbb{F}_q$.

A naive inspection of the Koksma-Hlawka upper bound $O((\log n)^d / n)$ suggests that QMC integration should fail catastrophically in high dimensions: if $d = 100$, $(\log n)^{100}$ does not decrease below $n$ until $n$ exceeds the number of subatomic particles in the observable universe.

Yet, in financial engineering and computational physics, Sobol' sequences routinely beat Crude Monte Carlo in 100-dimensional and 360-dimensional mortgage-backed security simulations. Ian H. Sloan and Henryk Woźniakowski (1998), along with Russel E. Caflisch, William J. Morokoff, and Art B. Owen (1997), solved this paradox by introducing the concept of Effective Dimension.

1. ANOVA Decomposition of Functions

Any square-integrable function $f: [0, 1]^d \to \mathbb{R}$ admits an exact orthogonal Analysis of Variance (ANOVA) decomposition:

where $f_0 = \int f(\mathbf{x})\,d\mathbf{x}$, each term $f_u$ depends only on the coordinate subset $u \subseteq \{1, \dots, d\}$, and all terms are mutually orthogonal in $L^2([0, 1]^d)$:

The total variance $\sigma^2 = \mathrm{Var}(f(\mathbf{X}))$ decomposes into variance components $\sigma_u^2 = \mathrm{Var}(f_u(\mathbf{X}_u))$:

2. Truncation Dimension and Superposition Dimension

Caflisch, Morokoff, and Owen defined two operational metrics:

1. Effective Dimension in the Truncation Sense: The smallest integer $d_T$ such that the variance explained by the first $d_T$ coordinates constitutes at least $1 - \alpha$ (typically $99\%$) of the total variance:

If a 360-dimensional problem has truncation dimension $d_T = 3$, only the first 3 coordinates drive $99\%$ of the variation, while the remaining 357 coordinates act as negligible perturbations.

2. Effective Dimension in the Superposition Sense: The smallest integer $d_S$ such that terms depending on at most $d_S$ variables explain at least $1 - \alpha$ of the total variance:

In many quantitative models, even if all 100 coordinates are equally active, functions are dominated by low-order interactions ($|u| \le 2$, meaning single-variable and pairwise effects). Because low-discrepancy sequences exhibit exceptional bivariate and trivariate projections, QMC integrates these low-order ANOVA components with near $O(n^{-1})$ accuracy.

A major practical drawback of pure deterministic QMC is that the Koksma-Hlawka bound provides an uncomputable worst-case deterministic error bound: the variation $V_{\text{HK}}(f)$ is impossible to evaluate analytically for complex simulation models. Consequently, deterministic QMC cannot produce a valid statistical standard error or confidence interval.

Randomized Quasi-Monte Carlo (RQMC) resolves this dilemma by randomizing a low-discrepancy point set $P_n$ such that:

1. Each individual randomized point $\mathbf{X}_i$ is uniformly distributed over the hypercube: $\mathbf{X}_i \sim \mathcal{U}([0, 1]^d)$.
2. The collective low-discrepancy structure of the ensemble is strictly preserved.

1. Cranley-Patterson Rotation (Random Shift)

The simplest RQMC randomization, introduced by R. Cranley and T. N. Patterson (1976), generates a single uniform shift vector $\mathbf{U} \sim \mathcal{U}([0, 1]^d)$ and shifts all deterministic points modulo 1:

By repeating this procedure independently $M$ times with distinct shift vectors $\mathbf{U}_1, \dots, \mathbf{U}_M$, we obtain $M$ independent, unbiased estimators $\hat{\theta}^{(1)}, \dots, \hat{\theta}^{(M)}$, where:

The grand RQMC estimator is the sample average:

Its empirical standard error is calculated directly using elementary sample variance formulas:

2. Owen's Scrambled Nets

Art B. Owen (1995) developed a deeper randomization scheme: Nested Scrambling. For digital $(t, m, s)$-nets in base $b$, Owen's method applies independent random permutations to the digits of each coordinate expansion, conditioned on prefix history.

Theorem (Owen): For any function $f$ whose partial derivatives satisfy Lipschitz continuity, integration over an Owen-scrambled net of size $n$ achieves a root mean square error (RMSE) of:

This improves upon the deterministic rate $O(n^{-1})$ and destroys the classical Monte Carlo rate $O(n^{-1/2})$, establishing RQMC as the premier integration technology for smooth simulation models.

In computational finance, mathematical biology, and fluid dynamics, models are defined by Stochastic Differential Equations (SDEs):

To estimate the expected payoff $\mathbb{E}[g(\mathbf{X}_T)]$, the continuous SDE must be discretized into $M$ discrete timesteps of size $h = T/M$ using the Euler-Maruyama or Milstein scheme, producing an approximation $\hat{X}_{T, h}$.

The total mean squared error (MSE) consists of two competing errors:

To achieve an accuracy tolerance of $\epsilon$, one must choose $h = O(\epsilon)$ and $N = O(\epsilon^{-2})$. Because each simulated path requires $M = O(h^{-1}) = O(\epsilon^{-1})$ operations, the total computational complexity of standard Monte Carlo integration is:

Michael B. Giles (2008) revolutionized computational simulation by introducing Multilevel Monte Carlo (MLMC), slashing the complexity to $O(\epsilon^{-2})$.

1. Telescoping Sum Decomposition

Let $L$ denote the finest discretization level, with geometric timestep refinement $h_l = h_0 \cdot M^{-l}$ for levels $l = 0, 1, \dots, L$ (typically with $M = 2$ or $4$). Let $P_l$ denote the simulation approximation of payoff $g(\mathbf{X}_T)$ computed with timestep $h_l$. By linearity of expectation, the expectation at the finest level $P_L$ can be expressed as a telescoping sum:

Instead of evaluating $P_L$ directly on a massive sample size, MLMC estimates each correction term $\mathbb{E}[P_l - P_{l-1}]$ independently using $N_l$ simulation paths:

The Crucial Coupling Trick: To evaluate $P_l^{(i)} - P_{l-1}^{(i)}$, the fine path $P_l$ and the coarse path $P_{l-1}$ are simulated using the exact same Brownian motion path $W_t$. The fine path evaluates increments $\Delta W$ at steps of size $h_l$, while the coarse path simply sums adjacent Brownian increments $\Delta W_{\text{coarse}} = \Delta W_1 + \dots + \Delta W_M$.

Because both paths are driven by identical Brownian realizations, $P_l^{(i)}$ and $P_{l-1}^{(i)}$ are intensely positively correlated. Consequently:

2. Giles' Complexity Theorem

Theorem (Giles 2008): Let $C_l = O(h_l^{-1})$ be the computational cost per path at level $l$, and suppose the variance decays as $V_l = O(h_l^\beta)$ with $\beta > 0$. The optimal allocation of sample sizes minimizing total variance subject to fixed cost is:

Under the Milstein discretization scheme, strong convergence is of order 1, yielding $\beta = 2$. Because $\beta > 1$, the variance decays faster than the cost increases. The sum of computational work across all levels is dominated by the coarsest level $l = 0$.

The total computational cost to achieve root mean squared error $\epsilon$ is bounded by:

MLMC eliminates the penalty of time-discretization, achieving the theoretical complexity bound of an exact continuous sampler.

In reliability engineering, telecommunications network buffering, and nuclear safety verification, simulation models must quantify the probability of catastrophic failure events occurring with probability $p \le 10^{-9}$. Crude Monte Carlo requires $n \approx 10^{11}$ runs to record 100 failure events, demanding years of supercomputer time.

While Importance Sampling provides one route to variance reduction, finding an effective proposal density $g(\mathbf{x})$ in high-dimensional dynamic network systems is often impossible. Multilevel Splitting, originally conceived by Herman Kahn and Theodore E. Harris (1951) and refined into the RESTART (Repetitive Simulation Trials After Reaching Thresholds) method by Manuel Villén-Altamirano et al. (1991), provides a powerful alternative that requires zero analytical knowledge of likelihood functions.

1. Nested Threshold Decomposition

Let $S(\mathbf{X}_t)$ be a continuous importance function or "score" measuring the proximity of system state $\mathbf{X}_t$ to the failure domain $\mathcal{F} = \{ \mathbf{x} : S(\mathbf{x}) \ge \gamma_M \}$.

The continuous state space is partitioned by a sequence of nested intermediate thresholds:

Define intermediate events $E_k = \{ \max_{0 \le t \le T} S(\mathbf{X}_t) \ge \gamma_k \}$, such that $E_M = \mathcal{F}$ and $E_M \subset E_{M-1} \subset \dots \subset E_1$.

By the definition of conditional probability, the target rare-event probability factors into a product of conditional transition probabilities:

The thresholds $\gamma_k$ are chosen such that each conditional probability is comfortably large:

Estimating ten conditional probabilities of magnitude $0.10$ requires only thousands of simulation draws per level, transforming an intractable $10^{-10}$ estimation problem into a series of easily solvable sub-problems.

2. The Splitting Mechanism

The algorithm executes iteratively across the threshold levels:

1. Level 1: Generate $N$ independent trajectories from the system initial state. Record the trajectories that successfully cross threshold $\gamma_1$. Let the empirical success fraction be $\hat{p}_1 = N_1 / N$.

2. Branching (Splitting): For each trajectory that reached $\gamma_1$, split the process into $R$ cloned sub-trajectories starting from the exact state $\mathbf{X}_{\tau_1}$ at the threshold-crossing time $\tau_1$.

3. Independent Continuation: Propagate each cloned sub-trajectory forward independently using fresh random numbers. Record the fraction of clones that successfully cross the next threshold $\gamma_2$:

1. Repeat the branching and propagation steps at each intermediate boundary until the failure domain $\gamma_M$ is reached.

The grand estimator of the failure probability is:

P. L'Ecuyer, F. Le Gland, and B. Tuffin (2007) proved that with optimal threshold spacing and clone branching ratios $R_k \approx 1 / p_k$, the computational effort scales as:

rather than the exponential scaling $O(1 / P(\mathcal{F}))$ of Crude Monte Carlo, enabling routine estimation of events with probabilities as small as $10^{-50}$.

In spatial network simulation, percolation theory, and image segmentation, statisticians simulate random graph structures, notably Uniform Spanning Trees (UST): drawing a spanning tree uniformly at random from the set of all spanning trees of a connected graph $G = (V, E)$.

While Kirchhoff's Matrix Tree Theorem allows computing the total count of spanning trees via the determinant of the graph Laplacian, direct combinatorial sampling from all trees is non-trivial.

David B. Wilson (1996) formulated an exact, cycle-erased simulation algorithm:

1. Pick an arbitrary root vertex $r \in V$ and add it to the initialized tree $T = \{r\}$.
2. While there exist unvisited vertices in $V \setminus T$:

• Pick an unvisited vertex $u \in V \setminus T$.
• Start a simple random walk from $u$.
• Whenever the random walk visits a vertex already in its current path, erase the created loop immediately (Loop-Erased Random Walk).
• Continue the walk until it hits any vertex already belonging to the existing tree $T$.
• Add the loop-erased path to tree $T$.

1. Output tree $T$.

Theorem (Wilson): The tree generated by Wilson's algorithm is distributed strictly uniformly over the set of all spanning trees of graph $G$. The expected running time of the algorithm is equal to the mean hitting time of the graph, outperforming earlier Aldous-Broder random walk algorithms by orders of magnitude.

To quantify the practical advantage of Quasi-Monte Carlo over pseudo-random Monte Carlo, we consider pricing a 40-dimensional arithmetic Asian basket option under multi-asset geometric Brownian motion.

Tracking convergence across sample sizes from $n = 10^3$ to $n = 10^6$:

• Crude Monte Carlo (PCG64): The standard error decays at the classic theoretical rate $O(n^{-0.50})$, requiring $n = 10^6$ paths to achieve a pricing precision of $\pm \$0.015$.
• Sobol' Low-Discrepancy Sequence with Brownian Bridge: By utilizing Sobol' points paired with Brownian bridge path construction (which allocates the most uniform low-discrepancy dimensions to the dominant early time steps), the pricing error decays at an observed empirical rate of $O(n^{-0.88})$.

The Sobol' simulation achieves the target accuracy of $\pm \$0.015$ with only $n = 18,500$ paths—reducing the required computational sample size by a factor of 54.

Unit 4 established the theoretical framework of Monte Carlo integration and its asymptotic convergence properties. We proved that Crude Monte Carlo integration achieves an error rate of $\sigma / \sqrt{n}$, rendering its convergence rate entirely independent of spatial dimension $d$ and overcoming the curse of dimensionality that cripples deterministic quadrature grids. We formulated Quasi-Monte Carlo (QMC) integration using low-discrepancy point sets (Halton, Sobol', and Faure sequences), proving via the Koksma-Hlawka inequality that deterministic integration error scales as $O((\log n)^d / n)$. We resolved the apparent paradox of high-dimensional QMC performance through ANOVA functional decompositions, demonstrating that models dominated by low truncation or superposition dimensions integrate with near-linear speed. We analyzed Randomized Quasi-Monte Carlo (RQMC) via Cranley-Patterson rotations and Owen's scrambled nets, which combine QMC speed with valid statistical standard error computation. Finally, we surveyed Giles' Multilevel Monte Carlo framework, which decouples time-discretization errors in stochastic differential equations.

When planning a production Monte Carlo study, the required sample size $n$ must be chosen rigorously to satisfy predefined relative precision bounds.

Let the relative error tolerance be $\delta > 0$ with nominal confidence level $1 - \alpha$:

Using the asymptotic normality of the Monte Carlo estimator, this condition translates to:

where $\mathrm{CV} = \sigma / |\theta|$ is the coefficient of variation of the integrand.

In preliminary exploratory pilot runs, an initial sample of size $n_0 \approx 1,000$ draws is evaluated to obtain a consistent plug-in estimate of the coefficient of variation $\widehat{\mathrm{CV}} = s_0 / |\bar{X}_0|$. If the estimated coefficient of variation is large (as occurs in rare-event models or highly volatile derivatives where $\mathrm{CV} \ge 10$), achieving a modest relative precision of $\delta = 1\%$ ($0.01$) with $95\%$ confidence requires:

Recognizing when raw sample size expansion becomes computationally prohibitive is the exact signal to transition to the variance reduction methods of Unit 5.