第 10 章 Monte Carlo 方法
上一章介绍了如何生成随机数。本章讨论如何用随机数解决计算问题。Monte Carlo 方法的核心思想很简单：如果一个量可以写成某个随机变量函数的期望，就可以用大量随机样本的平均值来近似它。
在经济统计中，很多问题都可以用这种思路处理。例如，复杂收入分布下的贫困率、税收规则的平均负担、预测模型的未来误差、抽样调查估计量的抽样分布，都可以通过模拟来近似。Monte Carlo 方法不一定给出解析公式，但它给出了可执行、可诊断、可复现的计算路径。
10.1 Monte Carlo 积分
设 \(X\) 的密度为 \(f(x)\)，我们希望计算
\[ \mu=E\{g(X)\}=\int g(x)f(x)\,dx. \]
若可以从 \(f(x)\) 中独立抽取样本 \(X_1,\ldots,X_N\)，则 Monte Carlo 估计量为
\[ \hat\mu_N=\frac{1}{N}\sum_{i=1}^N g(X_i). \]
其中 \(N\) 是模拟次数，\(g(X_i)\) 是第 \(i\) 次模拟得到的函数值。根据大数定律，当 \(N\) 增大时，\(\hat\mu_N\) 会趋近于 \(\mu\)。根据中心极限定理，
\[ \hat\mu_N \approx N\left(\mu,\frac{\sigma_g^2}{N}\right), \]
其中 \(\sigma_g^2=\operatorname{Var}\{g(X)\}\)。因此，Monte Carlo 标准误可以估计为
\[ \widehat{\operatorname{se}}(\hat\mu_N) = \frac{s_g}{\sqrt{N}}, \]
这里 \(s_g\) 是 \(g(X_1),\ldots,g(X_N)\) 的样本标准差。
10.2 例：补贴规则的平均支出
延续第 7 章的模拟收入分布。设家庭收入 \(Y\) 服从对数正态分布，低收入补贴函数为
\[ g(y)= \begin{cases} 2(1-y/4), & 0<y<4,\\ 0, & y\geq 4. \end{cases} \]
我们用 Monte Carlo 方法估计每户平均补贴支出。
set.seed(2026)
N <- 100000
income <- rlnorm(N, meanlog = log(8), sdlog = 0.6)
subsidy <- ifelse(income < 4, 2 * (1 - income / 4), 0)
mc_est <- mean(subsidy)
mc_se <- sd(subsidy) / sqrt(N)
c(estimate = mc_est,
mc_se = mc_se,
lower = mc_est - 1.96 * mc_se,
upper = mc_est + 1.96 * mc_se)
#> estimate mc_se lower upper
#> 0.0586503 0.0006198 0.0574354 0.0598652这个区间不是传统意义上由抽样调查误差形成的置信区间，而是 Monte Carlo 计算误差的近似区间。它回答的是：在给定模型和模拟次数下，随机模拟本身带来的误差有多大。
10.3 模拟次数与误差
Monte Carlo 误差通常按 \(1/\sqrt{N}\) 的速度下降。也就是说，如果希望误差大约减半，模拟次数需要增加到原来的 4 倍。这一点和许多数值积分方法不同：Monte Carlo 在低维问题中未必最快，但在高维问题中非常有价值。
set.seed(1)
Ns <- c(100, 1000, 10000, 100000)
out <- data.frame(N = Ns, estimate = NA, mc_se = NA)
for (j in seq_along(Ns)) {
y <- rlnorm(Ns[j], meanlog = log(8), sdlog = 0.6)
g <- ifelse(y < 4, 2 * (1 - y / 4), 0)
out$estimate[j] <- mean(g)
out$mc_se[j] <- sd(g) / sqrt(Ns[j])
}
out
#> N estimate mc_se
#> 1 1e+02 0.03701 0.0147815
#> 2 1e+03 0.07073 0.0069786
#> 3 1e+04 0.06045 0.0020103
#> 4 1e+05 0.05891 0.0006216实际工作中，不能只报告模拟均值，也要报告 Monte Carlo 标准误或至少说明模拟次数。否则读者无法判断结果的小数位是否有意义。
10.4 用模拟研究估计量性质
Monte Carlo 方法还可以用来研究估计量在重复抽样下的表现。假设我们关心一个简单消费函数
\[ C_i=\beta_0+\beta_1 Y_i+\varepsilon_i, \]
其中 \(C_i\) 是消费支出，\(Y_i\) 是收入，\(\varepsilon_i\) 是随机扰动。若真实边际消费倾向为\(\beta_1=0.6\)，我们可以反复生成样本并估计线性回归，观察 \(\hat\beta_1\) 的抽样分布。
set.seed(2)
R <- 1000
n <- 200
beta1_hat <- numeric(R)
for (r in seq_len(R)) {
income <- rlnorm(n, meanlog = log(8), sdlog = 0.5)
consumption <- 1 + 0.6 * income + rnorm(n, sd = 2)
beta1_hat[r] <- coef(lm(consumption ~ income))[2]
}
c(mean = mean(beta1_hat),
sd = sd(beta1_hat),
true_value = 0.6)
#> mean sd true_value
#> 0.60085 0.03018 0.60000这个实验帮助我们把“估计量是随机变量”这句话具体化。每次样本不同，估计结果也不同；重复模拟后，我们可以近似观察估计量的偏差、方差和均方误差。
10.5 重要性抽样
有时我们关心的是罕见事件，例如高收入尾部概率、极端亏损概率或违约概率。如果直接从原分布抽样，罕见事件出现次数很少，估计会很不稳定。重要性抽样通过改变抽样分布来提高效率。
设目标积分为
\[ \mu=\int g(x)f(x)\,dx. \]
如果从另一个更容易覆盖重要区域的密度 \(q(x)\) 抽样，则
\[ \mu=\int g(x)\frac{f(x)}{q(x)}q(x)\,dx =E_q\left\{g(X)\frac{f(X)}{q(X)}\right\}. \]
其中 \(w(x)=f(x)/q(x)\) 称为重要性权重。关键是 \(q(x)\) 要在 \(g(x)f(x)\) 贡献大的区域给出足够概率。
下面估计收入超过 20 万元的尾部概率。直接模拟可能需要很大样本，重要性抽样则从均值更高的对数正态分布中抽样，再用权重修正。
set.seed(3)
N <- 20000
threshold <- 20
# 直接 Monte Carlo
y_direct <- rlnorm(N, meanlog = log(8), sdlog = 0.6)
p_direct <- mean(y_direct > threshold)
# 重要性抽样：从更偏向高收入区域的分布中抽样
y_imp <- rlnorm(N, meanlog = log(12), sdlog = 0.6)
log_w <- dlnorm(y_imp, log(8), 0.6, log = TRUE) -
dlnorm(y_imp, log(12), 0.6, log = TRUE)
w <- exp(log_w)
p_imp <- mean((y_imp > threshold) * w)
se_imp <- sd((y_imp > threshold) * w) / sqrt(N)
c(direct = p_direct,
importance = p_imp,
importance_se = se_imp)
#> direct importance importance_se
#> 0.0645500 0.0639562 0.0009438重要性抽样不是简单地“换一个分布抽样”。如果 \(q(x)\) 选得不好，权重可能极端不稳定，估计反而变差。因此，实际使用时要检查权重分布，例如是否有少数样本权重特别大。
10.6 方差缩减的基本思想
Monte Carlo 方法的成本主要来自随机误差。除了增加 \(N\)，还可以通过更聪明的抽样降低方差。最简单的例子是对偶变量法。若要估计
\[ E\{h(U)\},\quad U\sim U(0,1), \]
可以同时使用 \(U\) 和 \(1-U\)：
\[ \hat\mu=\frac{1}{N}\sum_{i=1}^N \frac{h(U_i)+h(1-U_i)}{2}. \]
当 \(h\) 是单调函数时，两者往往负相关，从而降低方差。
set.seed(4)
R <- 1000
N <- 200
plain <- anti <- numeric(R)
for (r in seq_len(R)) {
u <- runif(N)
plain[r] <- mean(exp(u))
anti[r] <- mean((exp(u) + exp(1 - u)) / 2)
}
c(plain_sd = sd(plain), antithetic_sd = sd(anti))
#> plain_sd antithetic_sd
#> 0.034165 0.004358本书不系统展开所有方差缩减技术，但重要性抽样和对偶变量法说明了一个共同原则：模拟不是越蛮力越好，合理利用问题结构可以显著提高效率。