3.3 Composition and Table Look-Up Methods

Composition (Mixture) Method

where each component density \(f_j\) is easy to sample from, then:

1. Generate an index \(J\) with \(P(J = j) = p_j\).
2. Generate \(X\) from density \(f_J\).

The resulting \(X\) has density \(f\).

Table Look-Up (Alias Method)

For discrete distributions with finitely many outcomes, the alias method (due to A. J. Walker, 1977) allows sampling in \(O(1)\) time after an \(O(k)\) preprocessing step, where \(k\) is the number of possible values.

The idea is to construct a table with \(k\) entries, each containing a probability \(q_j\) and an “alias” index \(a_j\). To generate a variate:

1. Generate an integer \(J\) uniformly from \(\{1, 2, \ldots, k\}\).
2. Generate \(U \sim \text{Uniform}(0,1)\).
3. If \(U \leq q_J\), return \(J\); otherwise return \(a_J\).

The preprocessing step constructs the table so that the correct PMF is reproduced. This method is particularly valuable when many variates from the same distribution are needed.

3.4 Generating from Specific Distributions

Normal Distribution

Several specialised methods exist for generating standard normal variates:

are independent \(N(0,1)\) random variables.

The acceptance probability is \(\pi/4 \approx 0.785\).

Gamma Distribution

For \(X \sim \text{Gamma}(\alpha, 1)\) with integer \(\alpha\), the sum-of-exponentials method works: \(X = -\sum_{i=1}^{\alpha} \ln U_i\). For non-integer \(\alpha > 1\), acceptance-rejection methods are used. A widely-used algorithm is due to Ahrens and Dieter (1974).

For \(\alpha < 1\), one can use the identity: if \(Y \sim \text{Gamma}(\alpha + 1, 1)\) and \(U \sim \text{Uniform}(0,1)\) are independent, then \(X = Y \cdot U^{1/\alpha} \sim \text{Gamma}(\alpha, 1)\).

Poisson Distribution

For \(X \sim \text{Poisson}(\lambda)\), the inverse transform via the CDF is natural: starting from \(p = e^{-\lambda}\), \(F = p\), and \(x = 0\), generate \(U\); while \(U > F\), update \(p \leftarrow p \cdot \lambda / (x+1)\), \(F \leftarrow F + p\), and \(x \leftarrow x + 1\). Return \(x\). This is efficient for small \(\lambda\) but slow for large \(\lambda\), where one switches to acceptance-rejection or normal approximation methods.

3.5 Generating Multivariate Distributions

Multivariate Normal

To generate \(\mathbf{X} \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) in \(\mathbb{R}^d\):

1. Compute the Cholesky decomposition \(\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}^T\), where \(\mathbf{L}\) is lower triangular.
2. Generate \(\mathbf{Z} = (Z_1, \ldots, Z_d)^T\) with each \(Z_i \sim N(0,1)\) independently.
3. Set \(\mathbf{X} = \boldsymbol{\mu} + \mathbf{L}\mathbf{Z}\).

Then \(E[\mathbf{X}] = \boldsymbol{\mu}\) and \(\text{Cov}(\mathbf{X}) = \mathbf{L}\, E[\mathbf{Z}\mathbf{Z}^T]\, \mathbf{L}^T = \mathbf{L}\mathbf{I}\mathbf{L}^T = \boldsymbol{\Sigma}\).

Conditional Distribution Method

Sample \(X_1\) from its marginal, then \(X_2\) from its conditional given \(X_1\), and so on.

3.6 Copulas

A copula provides a way to model the dependence structure between random variables separately from their marginal distributions. This is invaluable in finance, insurance, and risk management.

Gaussian Copula

where \(\Phi_R\) is the CDF of the standard multivariate normal with correlation \(\mathbf{R}\) and \(\Phi^{-1}\) is the standard normal quantile function.

To sample from a Gaussian copula with given marginals \(F_1, \ldots, F_d\):

1. Generate \(\mathbf{Z} \sim N(\mathbf{0}, \mathbf{R})\) using the Cholesky method.
2. Set \(U_i = \Phi(Z_i)\) for \(i = 1, \ldots, d\). The vector \((U_1, \ldots, U_d)\) has the Gaussian copula.
3. Set \(X_i = F_i^{-1}(U_i)\) to obtain the desired marginals.

Student-\(t\) Copula

The \(t\) copula allows for heavier tails and stronger tail dependence than the Gaussian copula. If \(\mathbf{R}\) is a correlation matrix and \(\nu > 0\) is the degrees of freedom:

1. Generate \(\mathbf{Z} \sim N(\mathbf{0}, \mathbf{R})\).
2. Generate \(W \sim \chi^2_\nu / \nu\) independently.
3. Set \(\mathbf{T} = \mathbf{Z}/\sqrt{W}\). Each \(T_i\) has a univariate \(t_\nu\) distribution.
4. Set \(U_i = t_\nu(T_i)\) where \(t_\nu\) is the CDF of the \(t_\nu\) distribution.
5. Set \(X_i = F_i^{-1}(U_i)\).

Archimedean Copulas

where \(\psi : [0, \infty) \to [0, 1]\) is a generator function satisfying \(\psi(0) = 1\), \(\psi(\infty) = 0\), and appropriate monotonicity conditions.

Common examples include:

Sampling from Archimedean copulas in two dimensions can be done using the conditional method: generate \(U_1 \sim \text{Uniform}(0,1)\), then generate \(U_2\) from the conditional distribution \(C(u_2 \mid u_1) = \partial C / \partial u_1\) evaluated at \(u_1\).

3.7 Simulating Stochastic Processes

Brownian Motion

A standard Brownian motion (Wiener process) \(\{W(t) : t \geq 0\}\) satisfies: (i) \(W(0) = 0\), (ii) the increments \(W(t) - W(s) \sim N(0, t-s)\) for \(0 \leq s < t\), and (iii) non-overlapping increments are independent.

To simulate \(W\) at times \(0 = t_0 < t_1 < \cdots < t_n\):

1. For \(i = 1, \ldots, n\), generate \(Z_i \sim N(0,1)\) independently.
2. Set \(W(t_i) = W(t_{i-1}) + \sqrt{t_i - t_{i-1}}\, Z_i\).

For equal spacing \(\Delta t = t_i - t_{i-1}\), this reduces to \(W(t_i) = W(t_{i-1}) + \sqrt{\Delta t}\, Z_i\).

where \(\mu\) is the drift and \(\sigma\) is the volatility.

Poisson Process

A Poisson process \(\{N(t) : t \geq 0\}\) with rate \(\lambda\) satisfies: (i) \(N(0) = 0\), (ii) independent increments, and (iii) \(N(t) - N(s) \sim \text{Poisson}(\lambda(t-s))\) for \(s < t\).

The key simulation insight is that the inter-arrival times \(T_1, T_2, \ldots\) are iid \(\text{Exp}(\lambda)\). To simulate the arrival times on \([0, T]\):

1. Set \(t = 0\), \(n = 0\).
2. Generate \(U \sim \text{Uniform}(0,1)\) and set \(t \leftarrow t - \ln(U)/\lambda\).
3. If \(t > T\), stop. Otherwise, \(n \leftarrow n + 1\), record \(S_n = t\), and return to step 2.

For a non-homogeneous Poisson process with time-varying rate \(\lambda(t) \leq \lambda^*\), one can use thinning (Lewis and Shedler, 1979): generate arrivals from a homogeneous Poisson process at rate \(\lambda^*\), then independently accept each arrival at time \(t\) with probability \(\lambda(t)/\lambda^*\).

Chapter 4: Monte Carlo Estimation

4.1 Monte Carlo Estimators for Integrals

Multidimensional Integration

achieves standard error \(\sigma/\sqrt{n}\) regardless of dimension, where \(\sigma^2 = \text{Var}(g(\boldsymbol{X}))\). For \(s = 10\), even 10 points per dimension requires \(10^{10}\) function evaluations, while MC can achieve a 1% relative error with roughly \(10^4/\sigma^2\) samples.

This advantage makes Monte Carlo the method of choice for integrals over more than 3–4 dimensions, which arise naturally in Bayesian posterior computations, high-dimensional finance models, and particle transport simulations.

Hit-or-Miss Estimator

An alternative approach for estimating \(\theta = \int_0^1 h(x)\, dx\) (assuming \(0 \leq h(x) \leq 1\)) uses the geometric interpretation:

1. Generate \((U_1, U_2)\) uniformly on the unit square \([0,1]^2\).
2. If \(U_2 \leq h(U_1)\), count a “hit.”

Then \(\hat{\theta}_n = (\text{number of hits})/n\) is an unbiased estimator because \(P(U_2 \leq h(U_1)) = \int_0^1 h(x)\, dx = \theta\).

4.2 Properties of Monte Carlo Estimators

Unbiasedness

Consistency

By the strong law of large numbers, \(\hat{\theta}_n \to \theta\) almost surely as \(n \to \infty\). This means that with enough samples, the estimator converges to the true value.

Mean Squared Error

For the unbiased sample-mean estimator, \(\text{MSE}(\hat{\theta}_n) = \text{Var}(\hat{\theta}_n) = \sigma^2/n\), where \(\sigma^2 = \text{Var}(g(X))\).

Central Limit Theorem for MC Estimators

Confidence Intervals

where \(z_{\alpha/2}\) is the \((1-\alpha/2)\)-quantile of the standard normal distribution. For a 95% confidence interval, \(z_{0.025} = 1.96\).

Determining Sample Size

In practice, one runs a pilot study to estimate \(S\), then determines the required total sample size.

4.3 Examples and Applications

Estimating \(\pi\)

This is a hit-or-miss estimator.

Option Pricing

Queueing Systems

Consider a single-server queue with Poisson arrivals (rate \(\lambda\)) and exponential service times (rate \(\mu\)). While the \(M/M/1\) queue has a known steady-state distribution, more complex systems (e.g., \(G/G/k\) queues with priorities) lack analytical solutions and must be studied via simulation. Typical output measures include average waiting time, server utilisation, and maximum queue length.

Reliability

For a system of \(d\) components with independent random lifetimes \(T_1, \ldots, T_d\), the system lifetime \(T_{\text{sys}} = h(T_1, \ldots, T_d)\) depends on the system structure function \(h\). For complex configurations (mixtures of series and parallel), Monte Carlo simulation estimates the system reliability \(P(T_{\text{sys}} > t)\) by generating many independent replications of component lifetimes.

4.4 Quasi-Monte Carlo Methods

Standard Monte Carlo uses pseudo-random sequences which, while statistically uniform, contain random clusters and gaps. Quasi-Monte Carlo (QMC) methods replace the random sequence with a low-discrepancy sequence that fills the unit hypercube more evenly.

For a pseudo-random sequence, \(D_n^* = O(\sqrt{\log\log n/n})\) almost surely (by the law of the iterated logarithm). Low-discrepancy sequences achieve the superior rate \(D_n^* = O((\log n)^s / n)\).

Halton Sequence

For \(s = 2\): use base 2 for dimension 1 and base 3 for dimension 2. The first 8 Halton points in 2D are far more evenly distributed than 8 pseudo-random points.

Sobol’ Sequences

Sobol’ sequences are the most widely used low-discrepancy sequences in high-dimensional problems. They are constructed over the base-2 field using direction numbers that optimise equidistribution. Key properties:

• Dimension \(s\) up to several thousand is supported by precomputed direction number tables.
• They are scrambled in modern implementations (using random digital shifts) to enable unbiased estimation and valid confidence intervals.
• Convergence rate of the QMC estimator variance: \(O((\log n)^{2s}/n^2)\), compared to \(O(1/n)\) for standard MC.

QMC vs. MC: When to Use Each

QMC is most beneficial when the integrand is smooth and of low effective dimension. For rough integrands, rare event problems, or Markov chain simulations, standard MC or specialised techniques (such as importance sampling) are preferred.

Chapter 5: Variance Reduction Techniques

5.1 Motivation

The standard error of the Monte Carlo estimator decreases as \(O(n^{-1/2})\). To reduce the error by a factor of 10, one needs 100 times as many samples – which can be prohibitively expensive. Variance reduction techniques (VRTs) aim to produce estimators with smaller variance without increasing the number of replications, or equivalently, to achieve the same precision with fewer replications.

5.2 Antithetic Variates

The idea of antithetic variates is to introduce negative correlation between pairs of replications to reduce the variance of the average.

Since \(1 - U_i \sim \text{Uniform}(0,1)\) whenever \(U_i\) does, both halves are valid estimators.

5.3 Stratified Sampling

Stratified sampling partitions the sample space into non-overlapping strata and samples from each stratum separately, guaranteeing that each region of the space is represented proportionally.

Method

where \(U_{ji} \sim \text{Uniform}(0,1)\) independently.

With proportional allocation (\(n_j = n/k\) for all \(j\)), the total number of samples is \(n\) and the estimator simplifies.

5.4 Importance Sampling

Importance sampling is arguably the most powerful and flexible variance reduction technique. The idea is to sample from a different distribution that places more probability mass in the “important” regions of the integrand.

Derivation

The ratio \(w(x) = f(x)/h(x)\) is called the likelihood ratio or importance weight.

5.5 Control Variates

for some constant \(\beta\). Since \(E[Y^{(c)}] = E[Y] = \theta\) regardless of \(\beta\), the estimator is unbiased.

In practice, \(\beta^*\) is unknown and must be estimated from the simulation output (e.g., using a pilot run or from the same data).

5.6 Combined Methods

Variance reduction techniques can often be combined for even greater effect.

Antithetic + Control Variates

One can apply antithetic variates to generate negatively correlated pairs, and then apply a control variate adjustment to each pair. Care must be taken to estimate the optimal control coefficient appropriately.

Stratification + Importance Sampling

Different importance sampling densities can be used within different strata, optimising the proposal locally. This is particularly useful when the integrand has different characteristics in different regions of the domain.

5.7 Applications to Option Pricing

European Options

Antithetic variates replace each \(Z_i\) with the pair \((Z_i, -Z_i)\). Since \(\max(S(T) - K, 0)\) is an increasing function of \(Z\), the antithetic pairs are negatively correlated.

Asian Options

An Asian call option has payoff \(\max(\bar{S} - K, 0)\) where \(\bar{S} = \frac{1}{m}\sum_{j=1}^{m} S(t_j)\) is the arithmetic average of the stock price at \(m\) monitoring dates. No closed-form price exists, so Monte Carlo is the standard tool. Variance reduction is essential: the geometric average Asian option (which has a closed-form price) serves as an excellent control variate because it is highly correlated with the arithmetic average version.

Barrier Options

A knock-out barrier option becomes worthless if the stock price crosses a barrier level during the option’s life. Simulating this requires tracking the entire price path, and importance sampling can shift the path distribution to make barrier-crossing events more frequent, improving efficiency for deep out-of-the-money barrier options.

5.8 Rare Event Simulation

which requires \(n = O(1/\theta)\) samples for a fixed relative error.

Importance sampling is the primary tool for rare event simulation. The idea is to change the probability measure so that the rare event becomes common, then correct via importance weights. For example, in estimating the probability of ruin in insurance, or the probability of buffer overflow in telecommunications, exponential tilting (also called exponential change of measure) shifts the distribution to make large losses more likely.

5.9 Conditional Monte Carlo

Conditional Monte Carlo reduces variance by analytically computing the conditional expectation over one or more random variables, thereby eliminating their contribution to the variance.

5.10 Common Random Numbers

With CRN, \(\text{Cov}(Y_1, Y_2) > 0\), reducing the variance of the estimated difference. This technique does not reduce the variance of estimating either \(E[Y_1]\) or \(E[Y_2]\) alone, but it dramatically improves the comparison.

5.11 Summary of Variance Reduction Techniques

Appendix: Key Algorithms and Results Summary

Random Variate Generation: Decision Tree

Choosing the right generation method depends on the form of the target density:

Acceptance-Rejection: Choosing the Envelope

The acceptance probability is \(1/c\) where \(c = \sup_x f(x)/f_e(x)\). A smaller \(c\) means fewer proposals are rejected. Optimal envelope density \(f_e^*(x) \propto f(x)\) gives \(c = 1\), but then \(f_e^* = f\) (which we cannot sample from). In practice, \(f_e\) is chosen from a family of easily sampled distributions, and \(c\) is minimized by matching the shape of \(f\) as closely as possible.

Monte Carlo Error Summary

Variance Reduction: Quantitative Benchmarks

Notation

Markov Chain Monte Carlo: A Brief Introduction

The methods covered in this course generate independent random variates. In many statistical problems — particularly Bayesian inference — the target distribution is a high-dimensional density \(\pi(x)\) that cannot be directly sampled. Markov chain Monte Carlo (MCMC) methods construct a Markov chain \(X_0, X_1, X_2, \ldots\) with \(\pi\) as its stationary distribution; after a burn-in period, the chain values \(X_t\) approximately follow \(\pi\).

Metropolis–Hastings Algorithm

The Metropolis–Hastings (MH) algorithm is the most general MCMC method. Given the current state \(X_t = x\):

1. Propose \(Y \sim q(y \mid x)\) from a proposal distribution \(q\).
2. Compute the acceptance ratio \[ \alpha(x, y) = \min\!\left(1, \frac{\pi(y)\,q(x \mid y)}{\pi(x)\,q(y \mid x)}\right). \]
3. With probability \(\alpha(x, y)\), set \(X_{t+1} = y\); otherwise set \(X_{t+1} = x\).

The chain satisfies detailed balance \(\pi(x) q(y \mid x) \alpha(x,y) = \pi(y)q(x \mid y)\alpha(y,x)\), which ensures \(\pi\) is stationary. A key advantage is that only the ratio \(\pi(y)/\pi(x)\) is needed, so any normalizing constant in \(\pi\) cancels.

Gibbs Sampling

When \(\pi\) is a joint distribution of \((X_1, \ldots, X_d)\) and all full conditionals \(\pi(X_j \mid X_{-j})\) are available in closed form, Gibbs sampling cycles through coordinates, sampling each from its full conditional. Gibbs sampling is a special case of MH with acceptance probability 1.

MCMC vs. Direct Simulation

The effective sample size (ESS) for MCMC is \(n / \tau\) where \(\tau = 1 + 2\sum_{k=1}^\infty \text{Corr}(X_t, X_{t+k})\). Highly correlated chains have large \(\tau\) and small ESS, requiring more iterations to achieve the same accuracy as direct MC.