Chapter 3 Monte Carlo Methods
3.1 Variance Reduction Techniques
3.1.1 Control Variates
Let \(X\) be our target estimator of \(\ell = E[X]\). Suppose we have another random variable \(Y\) with a known expectation \(\mu_Y = E[Y]\). The control variate estimator is defined as: \[X_c = X + c(Y - \mu_Y)\]
Suppose you are estimating the average temperature of a room using a cheap sensor \(X\). If you also have a high-precision thermometer \(Y\) that fluctuates in sync with \(X\) (correlated) and whose average reading \(\mu_Y\) is known, you can correct your reading: if \(Y\) is running higher than its known average (\(Y - \mu_Y > 0\)), you subtract a proportional correction from \(X\). This “control” adjusts the reading to significantly reduce measurement noise (variance). where \(c\) is a constant. The optimal coefficient \(c^*\) that minimizes \(\text{Var}(X_c)\) is: \[c^* = -\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\] Substituting \(c^*\) back into the variance formula: \[\text{Var}(X_c^*) = \text{Var}(X) \left( 1 - \rho_{XY}^2 \right)\] where \(\rho_{XY}\) is the correlation coefficient.
- Step-by-Step Proof:
Expand the Variance: Using properties of variance for a sum of random variables: \[\text{Var}(X_c) = \text{Var}\big(X + c(Y - \mu_Y)\big) = \text{Var}(X) + 2c\text{Cov}(X, Y - \mu_Y) + c^2\text{Var}(Y - \mu_Y)\] Since \(\mu_Y\) is a constant, shifting by it does not affect variance or covariance: \[\text{Var}(X_c) = \text{Var}(X) + 2c\text{Cov}(X, Y) + c^2\text{Var}(Y)\]
Minimize with Respect to \(c\): To find the value of \(c\) that minimizes the variance, we differentiate \(\text{Var}(X_c)\) with respect to \(c\) and set the derivative to zero: \[\frac{d}{d c}\text{Var}(X_c) = 2\text{Cov}(X, Y) + 2c\text{Var}(Y) = 0\] Solving for \(c\) yields the optimal coefficient: \[c^* = -\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\]
Substitute \(c^*\) to Find Minimum Variance: Now we plug \(c^*\) back into the variance expression: \[\text{Var}(X_c^*) = \text{Var}(X) + 2\left(-\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\right)\text{Cov}(X, Y) + \left(-\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\right)^2 \text{Var}(Y)\] \[\text{Var}(X_c^*) = \text{Var}(X) - 2\frac{\text{Cov}(X, Y)^2}{\text{Var}(Y)} + \frac{\text{Cov}(X, Y)^2}{\text{Var}(Y)} = \text{Var}(X) - \frac{\text{Cov}(X, Y)^2}{\text{Var}(Y)}\]
Express in Terms of Correlation Coefficient (\(\rho_{XY}\)): By definition, the correlation coefficient is \(\rho_{XY} = \frac{\text{Cov}(X, Y)}{\sqrt{\text{Var}(X)\text{Var}(Y)}}\), which means \(\text{Cov}(X, Y)^2 = \rho_{XY}^2 \text{Var}(X)\text{Var}(Y)\). Substituting this: \[\text{Var}(X_c^*) = \text{Var}(X) - \frac{\rho_{XY}^2 \text{Var}(X)\text{Var}(Y)}{\text{Var}(Y)} = \text{Var}(X)\left(1 - \rho_{XY}^2\right)\] Since \(\rho_{XY}^2 \ge 0\), we have \(\text{Var}(X_c^*) \le \text{Var}(X)\), guaranteeing that the variance is reduced or remains unchanged. \(\blacksquare\)
3.1.1.1 Code Verification
We estimate \(\ell = E[e^U]\) where \(U \sim \text{Unif}(0, 1)\) using Control Variates with \(Y = U\). The optimal control coefficient is \(c^* = -\frac{(3 - e)/2}{1/12} \approx -1.6903\). The exact analytical variance reduction ratio is: \[\text{Ratio}_{\text{theory}} = \frac{1}{1 - \rho_{XY}^2} \approx 61.4269\]
u_samples <- runif(10000)
x_est <- exp(u_samples)
y_est <- u_samples
c_opt <- -(3 - exp(1)) / (2 * (1/12)) # Analytical optimal c
x_c <- x_est + c_opt * (y_est - 0.5)
emp_reduction <- var(x_est) / var(x_c)The empirical variance reduction ratio is 61.2169 (theoretical: 61.4269).
3.1.2 Antithetic Variates
If \(U_i \sim \text{Unif}(0, 1)\) are used to generate random variables \(X_1 = h(U_1)\), we can generate \(X_2 = h(1 - U_1)\). If \(\text{Cov}(X_1, X_2) < 0\), then \(\text{Var}(\bar{X}) < \frac{\text{Var}(X_1)}{2}\).
Suppose you are estimating the average height of a hill using random steps. If you take a step far to the right, you are likely on high ground; if your next step is deliberately taken at the mirror-opposite direction far to the left, you are likely on low ground. Averaging these two “opposite” (negatively correlated) steps cancels out the extreme values, giving a much more stable estimate of the average height than two independent steps.
- Step-by-Step Proof:
Variance of the Mean of Two Variables: Let \(X_1, X_2\) be identical in distribution, so \(\text{Var}(X_1) = \text{Var}(X_2) = \sigma^2\). The variance of their sample mean is: \[\text{Var}\left(\frac{X_1 + X_2}{2}\right) = \frac{1}{4}\text{Var}(X_1 + X_2) = \frac{1}{4}\left(\text{Var}(X_1) + \text{Var}(X_2) + 2\text{Cov}(X_1, X_2)\right)\] \[\text{Var}\left(\frac{X_1 + X_2}{2}\right) = \frac{1}{4}\left(2\sigma^2 + 2\text{Cov}(X_1, X_2)\right) = \frac{\sigma^2}{2} + \frac{\text{Cov}(X_1, X_2)}{2}\]
Comparison with Independent Sampling: If \(X_1, X_2\) were sampled independently, their covariance would be 0, yielding a variance of \(\frac{\sigma^2}{2}\). If we can design the coupling such that \(\text{Cov}(X_1, X_2) < 0\): \[\text{Var}\left(\frac{X_1 + X_2}{2}\right) = \frac{\sigma^2}{2} + \frac{\text{Cov}(X_1, X_2)}{2} < \frac{\sigma^2}{2}\] This proves that introducing a negative covariance between pairs yields a lower estimator variance than drawing independent samples. \(\blacksquare\)
3.1.2.1 Code Verification
We estimate \(E[e^U]\) using 5000 antithetic pairs (10000 function evaluations total) and compare it to 10000 independent samples:
u1 <- runif(5000)
x_anti <- (exp(u1) + exp(1 - u1)) / 2
var_anti <- var(x_anti) / 5000 # Variance of mean of 5000 pairs
u2 <- runif(10000)
var_ind <- var(exp(u2)) / 10000
anti_ratio <- var_ind / var_antiThe variance reduction ratio achieved by Antithetic Variates is 30.5032, showing a substantial variance reduction compared to independent sampling.
3.1.3 Conditional Monte Carlo (Rao-Blackwellization)
By the law of total variance, conditioning on a helper variable \(Y\) always reduces variance: \[\text{Var}(E[X \mid Y]) = \text{Var}(X) - E[\text{Var}(X \mid Y)] \le \text{Var}(X)\]
- Step-by-Step Proof of the Law of Total Variance:
Decompose the Variance: By definition, the variance of \(X\) is: \[\text{Var}(X) = E[X^2] - (E[X])^2\]
Condition the Expectation inside the Squared Term: Using the Law of Total Expectation, \(E[X] = E[E[X \mid Y]]\). We can expand \(E[X^2]\) by conditioning on \(Y\): \[E[X^2] = E[E[X^2 \mid Y]]\]
Rewrite the Conditional Second Moment: For any conditional random variable, the conditional variance is: \[\text{Var}(X \mid Y) = E[X^2 \mid Y] - (E[X \mid Y])^2 \implies E[X^2 \mid Y] = \text{Var}(X \mid Y) + (E[X \mid Y])^2\]
Substitute Back and Take Expectations: Taking expectations with respect to \(Y\): \[E[X^2] = E[E[X^2 \mid Y]] = E\big[\text{Var}(X \mid Y)\big] + E\big[(E[X \mid Y])^2\big]\]
Assemble the Variance of \(X\): \[\text{Var}(X) = E\big[\text{Var}(X \mid Y)\big] + E\big[(E[X \mid Y])^2\big] - (E[E[X \mid Y]])^2\] Notice that the last two terms are exactly the variance of the random variable \(E[X \mid Y]\) (since \(\text{Var}(Z) = E[Z^2] - (E[Z])^2\) where \(Z = E[X \mid Y]\)): \[\text{Var}(X) = E\big[\text{Var}(X \mid Y)\big] + \text{Var}\big(E[X \mid Y]\big)\] Rearranging terms: \[\text{Var}\big(E[X \mid Y]\big) = \text{Var}(X) - E\big[\text{Var}(X \mid Y)\big]\] Since variance is always non-negative, the expectation of conditional variance is non-negative (\(E[\text{Var}(X \mid Y)] \ge 0\)). Thus: \[\text{Var}\big(E[X \mid Y]\big) \le \text{Var}(X)\] This completes the proof. \(\blacksquare\)
3.1.3.1 Code Verification
We simulate \(X \sim N(0, 1)\) and \(Y \mid X \sim N(\rho X, 1 - \rho^2)\) with \(\rho = 0.8\). We compare the empirical variance of \(X\) with the conditional expectation \(E[X \mid Y] = \rho Y\).
x_sim <- rnorm(10000)
y_sim <- 0.8 * x_sim + rnorm(10000, sd = sqrt(1 - 0.8^2))
cond_exp <- 0.8 * y_sim
var_x <- var(x_sim)
var_cond <- var(cond_exp)
cond_ratio <- var_x / var_condThe variance ratio of the raw variable to the conditional expectation is 1.5638 (theoretical target: \(1 / \rho^2 = 1.5625\)).
3.1.4 Stratified Sampling (Neyman Optimal Allocation)
For a fixed total sample size \(N\), the allocation \(n_h\) that minimizes Stratified variance is the Neyman Allocation: \[n_h = N \frac{W_h \sigma_h}{\sum_{k=1}^L W_k \sigma_k}\]
3.1.4.1 Code Verification
We estimate the mean of a population divided into two strata (weights \(W_1 = 0.3, W_2 = 0.7\), and standard deviations \(\sigma_1 = 1, \sigma_2 = 4\)). We draw a total of \(N=1000\) samples.
w1 <- 0.3; w2 <- 0.7
s1 <- 1; s2 <- 4
# Neyman allocation
n1 <- round(1000 * (w1 * s1) / (w1 * s1 + w2 * s2))
n2 <- 1000 - n1
# Draw stratified samples
y1 <- rnorm(n1, mean = 0, sd = s1)
y2 <- rnorm(n2, mean = 0, sd = s2)
var_strat <- (w1^2 * var(y1)/n1) + (w2^2 * var(y2)/n2)
# SRS samples
srs_samples <- rnorm(1000, mean = 0, sd = sqrt(w1*s1^2 + w2*s2^2))
var_srs <- var(srs_samples) / 1000
strat_deff <- var_srs / var_stratThe design effect (variance ratio of SRS to Stratified) is 1.3456, confirming Stratified sampling achieves lower variance.
3.2 Multilevel Monte Carlo (MLMC)
By discretization on levels \(l=0,\dots,L\), MLMC estimates the finest level \(P_L\) as a telescoping sum: \[E[P_L] = E[P_0] + \sum_{l=1}^L E[P_l - P_{l-1}]\]
3.2.1 Code Verification
We estimate \(E[X_T]\) for geometric Brownian motion \(dX_t = 0.05 X_t dt + 0.2 X_t dW_t\) with \(x_0 = 1, T = 1\) using the 2-level Euler-Maruyama discretization from helpers.R.
mlmc_data <- run_mlmc_2level(5000)
var_coarse <- var(mlmc_data$coarse)
var_difference <- var(mlmc_data$fine - mlmc_data$coarse)The coarse path variance is 0.0402, whereas the variance of the difference between fine and coarse paths is 0.0004, demonstrating that the variance of the correction term is significantly smaller.
3.3 Importance Sampling & Sequential Methods
3.3.1 Importance Sampling
To estimate \(\ell = \int h(x) f(x) \, dx\), we sample from proposal \(g(x)\) and weight by \(w(x) = f(x)/g(x)\).
Suppose you want to estimate the average depth of shipwrecks in the ocean, but they are concentrated in a few deep, rare trenches. If you drop random search lines across the entire ocean (standard Monte Carlo), you will almost never hit a trench, leading to high variance. Instead, you bias your search by targeting the known trench areas (sampling from a proposal \(g(x)\)) and correct for this deliberate bias by weighting down the discoveries in proportion to how much more likely you were to search there (\(w(x) = f(x)/g(x)\)). This focuses your sampling effort where it actually matters.
- Step-by-Step Proof of the Optimal Proposal Density:
Expression for Variance of Importance Sampling Estimator: Let \(\hat{\ell}_{\text{IS}} = \frac{1}{N} \sum_{i=1}^N \frac{h(X_i)f(X_i)}{g(X_i)}\) where \(X_i \sim g\). The variance of the estimator is: \[\text{Var}(\hat{\ell}_{\text{IS}}) = \frac{1}{N} \text{Var}_g\left( \frac{h(X)f(X)}{g(X)} \right) = \frac{1}{N} \left[ E_g\left[ \left(\frac{h(X)f(X)}{g(X)}\right)^2 \right] - \ell^2 \right]\]
Minimize the Second Moment: Since \(\ell^2\) is constant (independent of \(g\)), minimizing the variance is equivalent to minimizing the first term: \[E_g\left[ \left(\frac{h(X)f(X)}{g(X)}\right)^2 \right] = \int \left(\frac{h(x)f(x)}{g(x)}\right)^2 g(x) \, dx = \int \frac{h(x)^2 f(x)^2}{g(x)} \, dx\]
Apply Cauchy-Schwarz Inequality: Recall that Cauchy-Schwarz states \(\left(\int u(x)v(x) \, dx\right)^2 \le \left(\int u(x)^2 \, dx\right)\left(\int v(x)^2 \, dx\right)\). Let \(u(x) = \frac{|h(x)|f(x)}{\sqrt{g(x)}}\) and \(v(x) = \sqrt{g(x)}\). Then: \[\left( \int \frac{|h(x)|f(x)}{\sqrt{g(x)}} \sqrt{g(x)} \, dx \right)^2 \le \left( \int \frac{h(x)^2 f(x)^2}{g(x)} \, dx \right) \left( \int g(x) \, dx \right)\] Since \(g(x)\) is a probability density, \(\int g(x) \, dx = 1\). The left side simplifies: \[\left( \int |h(x)|f(x) \, dx \right)^2 \le \int \frac{h(x)^2 f(x)^2}{g(x)} \, dx\]
Condition for Equality: The lower bound is achieved (equality holds in Cauchy-Schwarz) if and only if \(u(x) / v(x)\) is constant: \[\frac{\frac{|h(x)|f(x)}{\sqrt{g(x)}}}{\sqrt{g(x)}} = c_{\text{const}} \implies \frac{|h(x)|f(x)}{g(x)} = c_{\text{const}}\] \[g^*(x) = \frac{|h(x)|f(x)}{\int |h(y)|f(y) \, dy}\] This proves that the optimal proposal density \(g^*(x)\) is directly proportional to \(|h(x)|f(x)\). \(\blacksquare\)
3.3.1.1 Code Verification
We estimate \(\int_0^1 \cos(x) \, dx\) (analytical value: \(\sin(1) \approx 0.8415\)) using a linear proposal density \(g(x) = \frac{4 - 2x}{3}\) on \([0, 1]\), which is a close match for the cosine function:
u <- runif(10000)
x_proposal <- 2 - sqrt(4 - 3 * u) # Inverse transform of g(x)
weight <- cos(x_proposal) / ((4 - 2 * x_proposal) / 3)
imp_mean <- mean(weight)
imp_var <- var(weight) / 10000
srs_var <- var(cos(runif(10000))) / 10000
imp_ratio <- srs_var / imp_varThe Importance Sampling estimate is 0.8414 (theoretical: 0.8415), and the variance reduction factor (ratio of SRS variance to IS variance) is 12.3293, demonstrating a substantial variance reduction (approx. 12.5x) over standard random sampling.
3.3.2 Sequential Importance Sampling (SIS) Degeneracy
At step \(t\), the weight update is: \[w_t = w_{t-1} \frac{f(x_t \mid x_{1:t-1})}{q(x_t \mid x_{1:t-1})}\]
3.3.2.1 Code Verification
We run a 20-step SIS algorithm on a state space without resampling and print the weight Effective Sample Size (ESS) to illustrate weight collapse:
N <- 1000
w <- rep(1/N, N)
for (t in 1:20) {
# Stochastically update weights representing degeneracy
w <- w * runif(N, 0, 2)
w <- w / sum(w)
}
sis_ess <- 1 / sum(w^2)The final weight Effective Sample Size (ESS) after 20 steps is 11.43 (starting at 1000.00), demonstrating severe weight degeneracy.
3.3.3 Nonlinear Filtering (Particle Filtering)
We filter a state-space model using the Bootstrap particle filter from helpers.R.
3.3.3.1 Code Verification
We run the particle filter on a simulated 50-step trajectory:
# Generate true state and observations
x_true <- numeric(50)
y_obs <- numeric(50)
x_true[1] <- rnorm(1)
y_obs[1] <- x_true[1] + rnorm(1)
for (t in 2:50) {
x_true[t] <- 0.5 * x_true[t-1] + rnorm(1)
y_obs[t] <- x_true[t] + rnorm(1)
}
x_est <- run_particle_filter(y_obs, N = 500)
rmse_filter <- sqrt(mean((x_est - x_true)^2))
rmse_obs <- sqrt(mean((y_obs - x_true)^2))The RMSE of the filtered state estimates is 0.7638, which is lower than the raw observation noise RMSE of 1.0173.
3.4 Transform Likelihood Ratio Methods
The sensitivity gradient is: \[\nabla_\theta E_\theta[h(X)] = E_\theta[h(X) \nabla_\theta \log f(X; \theta)]\]
3.4.1 Code Verification
We estimate the derivative of \(E_\theta[X^2]\) where \(X \sim N(\theta, 1)\) at \(\theta = 2\) (analytical derivative: \(2\theta = 4.0000\)) using the Likelihood Ratio estimator:
x_samples <- rnorm(10000, mean = 2, sd = 1)
# Score function for mean mu: (x - mu) / sigma^2
score_fun <- x_samples - 2
lr_grad_est <- mean(x_samples^2 * score_fun)The Likelihood Ratio gradient estimate is 4.0828 (analytical value: 4.0000).
3.5 Resampling & Empirical Likelihood
3.5.1 Bootstrap
We estimate the standard error of the correlation coefficient for a small dataset.
3.5.1.1 Code Verification
We draw 100 samples with correlation 0.5 and run 1000 bootstrap replicates:
x <- rnorm(100)
y <- 0.5 * x + rnorm(100, sd = sqrt(1 - 0.5^2))
orig_cor <- cor(x, y)
boot_cor <- bootstrap_resample(1:100, function(idx) cor(x[idx], y[idx]), R = 1000)
boot_se <- sd(boot_cor)The estimated bootstrap standard error of the correlation coefficient is 0.0871.
3.5.2 Jackknife
We estimate the bias of the sample variance \(S^2\) for a small sample size.
3.5.2.1 Code Verification
We draw 15 samples from a Normal distribution:
x_data <- rnorm(15)
n <- length(x_data)
jack_est <- numeric(n)
for (i in 1:n) {
jack_est[i] <- var(x_data[-i])
}
jack_bias <- (n - 1) * (mean(jack_est) - var(x_data))The estimated Jackknife bias of the sample variance is 0.0000 (theoretical expectation is near 0).
3.5.3 Cross-Validation
Cross-Validation is a resampling method used to assess how well a statistical model generalizes to independent, unseen data.
Suppose you are a student preparing for a final exam. If you study using only the mock exam questions, you might memorize the specific questions without truly understanding the concepts, leading to a poor grade on the actual exam (overfitting). To test your true preparation, you hide one chapter’s mock questions, study the rest, and then test yourself on the hidden chapter. In statistical modeling, \(k\)-fold Cross-Validation partitions the dataset into \(k\) equal parts: the model is trained on \(k-1\) folds and tested on the remaining fold. By cycling this process, we estimate how well the model generalizes to unseen data.
Practical Applications: Used to select optimal model hyper-parameters (such as regularization parameters in Lasso/Ridge regression or number of trees in random forests) and prevent model overfitting in clinical predictive models.
Mathematical Formulation: Let the dataset be partitioned into \(k\) non-overlapping folds \(D_1, \dots, D_k\). For each fold \(j \in \{1, \dots, k\}\), the model is trained on the dataset excluding fold \(j\) (\(D \setminus D_j\)) to obtain prediction function \(\hat{f}^{(-j)}(\mathbf{x})\). The cross-validation estimate of the prediction error (e.g. Mean Squared Error) is: \[\text{CV}_k = \frac{1}{n} \sum_{i=1}^n \left(y_i - \hat{f}^{(-j(i))}(\mathbf{x}_i)\right)^2\] where \(j(i)\) is the index of the fold containing sample \(i\), and \(n\) is the total sample size.
Manual Walkthrough (By-Hand Simulation): We trace \(k\)-fold cross-validation with \(k=3\) folds for a dataset of \(n=3\) points: \(d_1 = (1, 2)\), \(d_2 = (2, 4)\), and \(d_3 = (3, 5)\). The model is a simple average predictor \(\hat{f}(x) = \bar{y}\).
Fold 1 (\(D_1 = \{d_1\}\)): Train on \(D \setminus D_1 = \{d_2, d_3\} = \{(2, 4), (3, 5)\} \implies \hat{f}^{(-1)}(x) = \frac{4+5}{2} = 4.5\). Test on \(d_1\): Squared error is \((y_1 - 4.5)^2 = (2 - 4.5)^2 = 6.25\).
Fold 2 (\(D_2 = \{d_2\}\)): Train on \(D \setminus D_2 = \{d_1, d_3\} = \{(1, 2), (3, 5)\} \implies \hat{f}^{(-2)}(x) = \frac{2+5}{2} = 3.5\). Test on \(d_2\): Squared error is \((y_2 - 3.5)^2 = (4 - 3.5)^2 = 0.25\).
Fold 3 (\(D_3 = \{d_3\}\)): Train on \(D \setminus D_3 = \{d_1, d_2\} = \{(1, 2), (2, 4)\} \implies \hat{f}^{(-3)}(x) = \frac{2+4}{2} = 3.0\). Test on \(d_3\): Squared error is \((y_3 - 3.0)^2 = (5 - 3.0)^2 = 4.0\).
Estimated CV Mean Squared Error: \[\text{CV}_3 = \frac{6.25 + 0.25 + 4.0}{3} \approx 3.50\]
3.5.3.1 Cross-Validation Code Verification
We simulate \(n=50\) samples from a linear model \(Y = 2X + \epsilon\) with \(\epsilon \sim N(0, 1)\) (noise variance is exactly 1.0). We run a custom 5-fold cross-validation in R to estimate the prediction MSE:
# Generate simulated data
set.seed(12345)
x_cv <- runif(50, 0, 10)
y_cv <- 2 * x_cv + rnorm(50, mean = 0, sd = 1)
data_cv <- data.frame(x = x_cv, y = y_cv)
# 5-Fold Cross-Validation using helper function
cv_mse <- run_kfold_cv(data_cv, k = 5)The estimated 5-fold Cross-Validation Mean Squared Error (MSE) is 1.5390 (theoretical expected prediction error is near 1.0, representing the irreducible noise variance \(\sigma^2 = 1.0\)).
3.5.4 Empirical Likelihood
We evaluate the Empirical Profile Likelihood ratio for the mean.
3.5.4.1 Code Verification
We evaluate the profile likelihood ratio statistic \(-2 \log R(\theta)\) for the true mean of a normal sample:
x_data <- rnorm(50, mean = 0, sd = 1)
el_stat <- calculate_empirical_likelihood_mean(x_data, mu_test = 0)The empirical profile likelihood ratio statistic \(-2 \log R(\theta_0)\) is 0.0034, which is well below the critical value of \(\chi^2_1(0.95) = 3.8415\) (confirming the hypothesis \(H_0: \mu = 0\) is accepted).
3.6 How to Report This in a Paper
When publishing Monte Carlo variance reduction studies:
- Report the variance reduction factor (ratio of standard Monte Carlo variance to the reduced variance).
- Report the theoretical correlation \(\rho_{XY}\) that justified the choice of the control variate.
- For bootstrap confidence intervals, specify whether percentile, studentized, or BCa intervals were used.
Key Formulas Summary
| Technique / Estimator | Formula / Definition | Key Constraints / Details |
|---|---|---|
| Control Variates | \(X_c = X + c(Y - \mu_Y)\) | \(\mu_Y = E[Y]\) must be known; optimal \(c^* = -\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\). |
| CV Variance Reduction | \(\text{Var}(X_c^*) = \text{Var}(X)(1 - \rho_{XY}^2)\) | Variance reduces by factor of correlation squared. |
| Antithetic Variates | \(\bar{X}_{\text{anti}} = \frac{h(U) + h(1 - U)}{2}\) | Requires \(\text{Cov}(h(U), h(1 - U)) < 0\). |
| Importance Sampling | \(\hat{\ell}_{\text{IS}} = \frac{1}{N} \sum_{i=1}^N h(X_i)\frac{f(X_i)}{g(X_i)}\) | Samples drawn from \(g(x)\); weight function \(w(x) = \frac{f(x)}{g(x)}\). |
| Optimal IS Proposal | \(g^*(x) \propto |h(x)|f(x)\) | Achieves zero variance if \(h(x) \ge 0\). |
| Bootstrap Filter Update | \(w_t^{(i)} \propto w_{t-1}^{(i)} \cdot p(y_t \mid x_t^{(i)})\) | Sequential MC update weight in state-space models. |
| Empirical Likelihood | \(R(\theta) = \max \left\{ \prod_{i=1}^N N p_i \right\}\) | Constraints: \(\sum p_i = 1\), \(\sum p_i (x_i - \theta) = 0\). |
| Cross-Validation | \(\text{CV}_k = \frac{1}{n} \sum_{i=1}^n \left(y_i - \hat{f}^{(-j(i))}(\mathbf{x}_i)\right)^2\) | Partitions dataset into \(k\) folds to estimate generalization error. |
Common Mistakes
Control Variates Sign Confusion: The formula for the control variates estimator can be written as \(X_c = X + c(Y - \mu_Y)\) or \(X_c = X - c(Y - \mu_Y)\).
- If you use the addition form (\(X + c(\dots)\)), then the optimal coefficient is \(c^* = -\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\).
- If you use the subtraction form (\(X - c(\dots)\)), then the optimal coefficient is \(c^* = \frac{\text{Cov}(X, Y)}{\text{Var}(Y)}\). Mixing these signs will flip the correction direction, adding noise and potentially doubling the variance instead of reducing it!
Antithetic Monotonicity: Antithetic sampling is only guaranteed to reduce variance if the transformation function \(h(u)\) is monotonic in \(u\). If \(h(u)\) is non-monotonic (e.g. a symmetric quadratic function), \(h(u)\) and \(h(1-u)\) can be positively correlated, which actually increases the estimator’s variance.
Fat-Tailed Importance Proposal: In Importance Sampling, the proposal distribution \(g(x)\) must have “fatter tails” than the target \(f(x)\). If \(g(x)\) goes to zero faster than \(f(x)\) in the tails, the weights \(w(x) = f(x)/g(x)\) will blow up to infinity, rendering the estimator’s variance infinite.
3.7 Exercises
3.7.1 Subsection 3.1: Variance Reduction Techniques
- Derive the optimal control coefficient \(c^*\) when using two control variates \(Y_1, Y_2\).
- Implement a Stratified sampling program in R for a 3-strata population and compare it with Neyman allocation vs proportional allocation.
3.7.2 Subsection 3.2: Importance Sampling & HMMs
- Prove that the optimal proposal density in Importance Sampling is proportional to \(|h(x)| f(x)\).
- Implement a bootstrap filter in R for a state-space model with non-linear transitions.
3.7.3 Subsection 3.3: Resampling & Likelihoods
- Code a double-bootstrap algorithm in R to estimate the bias of a bootstrap standard error.
- Compare the empirical likelihood confidence interval for a population median with the standard Student’s t-interval.
Chapter 3 covered:
Monte Carlo integration: Approximating \(\int h(x)f(x)\,dx\) by the sample mean \(\hat{\ell}_N = \frac{1}{N}\sum h(X_i)\), with standard error \(\sigma_h / \sqrt{N}\)
Control variates: Using a correlated variable \(Y\) with known mean \(\mu_Y\) to correct the estimate; optimal coefficient \(c^* = - ext{Cov}(X,Y)/ ext{Var}(Y)\); variance reduction factor \((1 - ho^2)\)
Antithetic variates: Pairing \(U\) and \(1-U\) to introduce negative correlation between estimates, halving the effective number of samples needed
Importance sampling: Sampling from an alternative density \(g(x)\) concentrated in high-impact regions; self-normalised estimator for unknown normalizing constants
Stratified sampling: Partitioning the sample space and sampling proportionally to reduce between-stratum variance
| Technique | Key Idea | When to Use |
|---|---|---|
| Control variates | Exploit a correlated known-mean variable | When a good correlated control is available |
| Antithetic variates | Pair \(U\) with \(1-U\) | Monotone integrands |
| Importance sampling | Sample from a better proposal | Rare events, heavy tails |
| Stratified sampling | Partition the domain | When variance varies across regions |
Coming up: Chapter 4 goes beyond Monte Carlo integration to full Bayesian inference: Markov Chain Monte Carlo allows you to draw samples from arbitrary posterior distributions that cannot be directly simulated.
Exercises
Tier 1 — Conceptual
- Explain in words why the standard error of a Monte Carlo estimate decreases at the rate \(1/\sqrt{N}\), not \(1/N\). What does this mean for doubling accuracy?
- The importance sampling estimator uses weights \(w(x) = f(x)/g(x)\). What happens to the estimator if you choose \(g(x)\) to be very different from \(f(x)\) in the tails?
- Antithetic variates exploit the fact that \(U\) and \(1-U\) are negatively correlated when passed through a monotone function \(h\). Give an example of a non-monotone \(h\) for which antithetic variates would fail to reduce variance.
Tier 2 — Applied
- Estimate \(\int_0^2 \sqrt{1 + x^4}\,dx\) using plain Monte Carlo with \(N = 50{,}000\) samples. Report your estimate, standard error, and 95% confidence interval.
- Repeat Question 4 using antithetic variates with \(N/2 = 25{,}000\) antithetic pairs. Compare the standard errors.
- Implement importance sampling to estimate \(P(Z > 4)\) where \(Z \sim N(0,1)\), using an exponential shifted proposal. Compare with the crude Monte Carlo estimate using \(N = 10^5\) samples.
Tier 3 — Challenge
- Option Pricing: The Black-Scholes price of a European call option can be computed as \(C = e^{-rT} E[\max(S_T - K, 0)]\) where \(S_T = S_0 \exp((r - \sigma^2/2)T + \sigma\sqrt{T}Z)\) and \(Z \sim N(0,1)\). With \(S_0 = 100, K = 105, r = 0.05, \sigma = 0.20, T = 1\), estimate \(C\) using: (a) plain Monte Carlo, (b) control variates with \(S_T\) as the control, (c) antithetic variates. Use \(N = 100{,}000\) and compare standard errors. The exact Black-Scholes price is 8.02.
Solutions to Tier 1 and 2: see Appendix B.