第 10 章 Monte Carlo 方法
上一章介绍了如何生成随机数。本章讨论如何用随机数解决计算问题。Monte Carlo 方法的核心思想很简单：如果一个量可以写成某个随机变量函数的期望，就可以用大量随机样本的平均值来近似它。
在经济统计中，很多问题都可以用这种思路处理。例如，复杂收入分布下的贫困率、税收规则的平均负担、预测模型的未来误差、抽样调查估计量的抽样分布，都可以通过模拟来近似。Monte Carlo 方法不一定给出解析公式，但它给出了可执行、可诊断、可复现的计算路径。
本章的学习重点有四个。第一，把积分、概率和统计量分布写成期望；第二，理解 Monte Carlo 误差及其标准误；第三，会设计重复模拟实验来研究估计量性质；第四，初步掌握重要性抽样和方差缩减思想。
10.1 Monte Carlo 积分
设 \(X\) 的密度为 \(f(x)\)，我们希望计算
\[ \mu=E\{g(X)\}=\int g(x)f(x)\,dx. \]
其中 \(X\) 是随机变量，\(f(x)\) 是其概率密度，\(g(\cdot)\) 是我们关心的函数，\(\mu\) 是目标期望或积分。若可以从 \(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 \mathcal{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)\) 的样本标准差。这个标准误衡量的是“模拟平均值因为随机数而产生的波动”，不是原始经济数据的抽样误差。
很多看起来不同的问题都可以写成这个形式。例如，若要估计事件 \(A\) 的概率，可以令
\[ g(X)=I(X\in A), \]
则
\[ P(X\in A)=E\{I(X\in A)\}. \]
若要计算普通积分
\[ I=\int_a^b h(x)\,dx, \]
可以令 \(U\sim U(a,b)\)，则
\[ I=(b-a)E\{h(U)\}. \]
这说明 Monte Carlo 不只是“做很多次模拟”，而是把难算的量改写成可以由样本平均近似的期望。
Monte Carlo 估计
- 明确目标量 \(\mu=E\{g(X)\}\)。
- 设计随机变量生成方法，从目标分布或合适的辅助分布中生成 \(X_1,\ldots,X_N\)。
- 计算模拟值 \(z_i=g(X_i)\)，\(i=1,\ldots,N\)。
- 用 \(\hat\mu_N=\bar z\) 估计 \(\mu\)。
- 用 \(s_z/\sqrt{N}\) 报告 Monte Carlo 标准误，并说明模拟次数 \(N\) 和随机种子。
下面先写一个小函数，用于汇总 Monte Carlo 估计值和模拟误差。
10.2 例：补贴规则的平均支出
延续第 7 章的模拟收入分布。设家庭收入 \(Y\) 服从对数正态分布，低收入补贴函数为
\[ g(y)= \begin{cases} 2(1-y/4), & 0<y<4,\\ 0, & y\geq 4. \end{cases} \]
其中 \(Y\) 可以理解为家庭年收入，单位为万元；\(g(y)\) 是该收入水平下的补贴金额，单位同样为万元。我们希望估计每户平均补贴支出
\[ \mu=E\{g(Y)\}. \]
如果直接积分不方便，可以用 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_summary(subsidy)
#> estimate mc_se lower upper
#> 0.0586503 0.0006198 0.0574354 0.0598651这个区间不是传统意义上由抽样调查误差形成的置信区间，而是 Monte Carlo 计算误差的近似区间。它回答的是：在给定收入模型和模拟次数下，随机模拟本身带来的误差有多大。若模型设定改变，例如收入分布的均值、离散程度或补贴规则改变，估计目标 \(\mu\) 本身也会改变。
类似地，贫困率也可以写成指示函数的期望。若贫困线设为 4 万元，则贫困率为
\[ p=P(Y<4)=E\{I(Y<4)\}. \]
poverty_indicator <- as.numeric(income < 4)
mc_summary(poverty_indicator)
#> estimate mc_se lower upper
#> 0.124300 0.001043 0.122255 0.126345对于二元指示变量，Monte Carlo 标准误也可以写成
\[ \sqrt{\frac{\hat p(1-\hat p)}{N}}, \]
其中 \(\hat p\) 是模拟得到的比例。这个公式与样本比例的标准误形式相同，但这里的随机性来自模拟次数。
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_real_, mc_se = NA_real_)
for (j in seq_along(Ns)) {
y <- rlnorm(Ns[j], meanlog = log(8), sdlog = 0.6)
z <- ifelse(y < 4, 2 * (1 - y / 4), 0)
out$estimate[j] <- mean(z)
out$mc_se[j] <- sd(z) / sqrt(Ns[j])
}
out_show <- out
out_show$N <- format(out_show$N, scientific = FALSE, trim = TRUE)
knitr::kable(out_show, row.names = FALSE, digits = c(0, 4, 4))| N | estimate | mc_se |
|---|---|---|
| 100 | 0.0370 | 0.0148 |
| 1000 | 0.0707 | 0.0070 |
| 10000 | 0.0605 | 0.0020 |
| 100000 | 0.0589 | 0.0006 |
实际工作中，不能只报告模拟均值，也要报告 Monte Carlo 标准误或至少说明模拟次数。否则读者无法判断结果的小数位是否有意义。例如，如果 Monte Carlo 标准误为 0.01，那么报告 0.376241 并没有太多意义，通常报告到 0.38 或 0.376 已经足够。
有时可以用一次较小规模的预模拟来估计需要的模拟次数。若希望 95% Monte Carlo 误差区间的半宽度不超过 \(\epsilon\)，可近似要求
\[ 1.96\frac{s_g}{\sqrt{N}}\leq \epsilon, \]
即
\[ N\geq \left(\frac{1.96s_g}{\epsilon}\right)^2. \]
这里 \(s_g\) 可由预模拟得到。下面以平均补贴支出为例，估计若希望半宽度不超过 0.002 万元，大约需要多少次模拟。
set.seed(7)
N_pilot <- 3000
y_pilot <- rlnorm(N_pilot, meanlog = log(8), sdlog = 0.6)
z_pilot <- ifelse(y_pilot < 4, 2 * (1 - y_pilot / 4), 0)
epsilon <- 0.002
required_N <- ceiling((1.96 * sd(z_pilot) / epsilon)^2)
sample_size_plan <- data.frame(
pilot_sd = sd(z_pilot),
epsilon = epsilon,
required_N = required_N
)
knitr::kable(sample_size_plan, row.names = FALSE,
digits = c(4, 4, 0))| pilot_sd | epsilon | required_N |
|---|---|---|
| 0.1887 | 0.002 | 34215 |
这个公式只是近似规划工具。实际报告时，仍应基于最终模拟结果重新计算 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\) 的抽样分布。
设第 \(r\) 次模拟得到估计量 \(\hat\theta^{(r)}\)，真实值为 \(\theta\)，重复次数为 \(R\)。常见评价指标为
\[ \widehat{\operatorname{Bias}} =\frac{1}{R}\sum_{r=1}^R\left(\hat\theta^{(r)}-\theta\right), \]
\[ \widehat{\operatorname{Var}} =\frac{1}{R-1}\sum_{r=1}^R \left(\hat\theta^{(r)}-\bar{\hat\theta}\right)^2, \]
以及
\[ \widehat{\operatorname{MSE}} =\frac{1}{R}\sum_{r=1}^R \left(\hat\theta^{(r)}-\theta\right)^2. \]
这里 \(\bar{\hat\theta}=R^{-1}\sum_{r=1}^R\hat\theta^{(r)}\)。偏差衡量估计量平均是否偏离真实值，方差衡量重复抽样下是否稳定，均方误差同时考虑偏差和方差。
set.seed(2)
simulate_beta1 <- function(n, R = 1000) {
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(n = n,
mean = mean(beta1_hat),
bias = mean(beta1_hat) - 0.6,
sd = sd(beta1_hat),
mse = mean((beta1_hat - 0.6)^2))
}
ols_out <- as.data.frame(
do.call(rbind, lapply(c(50, 200, 1000), simulate_beta1))
)
knitr::kable(ols_out, row.names = FALSE, digits = 4)| n | mean | bias | sd | mse |
|---|---|---|---|---|
| 50 | 0.6001 | 1e-04 | 0.0624 | 0.0039 |
| 200 | 0.6001 | 1e-04 | 0.0303 | 0.0009 |
| 1000 | 0.5999 | -1e-04 | 0.0134 | 0.0002 |
这个实验帮助我们把“估计量是随机变量”这句话具体化。每次样本不同，估计结果也不同；重复模拟后，我们可以近似观察估计量的偏差、方差和均方误差。注意这里有两层误差：每个样本中的回归估计有抽样误差；我们用有限的 \(R\) 次模拟总结这些性质，又会产生 Monte Carlo 误差。
10.5 例：信用组合损失模拟
Monte Carlo 在经济金融风险分析中也很常见。考虑一个由 \(m\) 个小微企业贷款构成的组合。第 \(j\) 个企业的贷款暴露为 \(E_j\)，违约指示变量为 \(D_j\)，违约损失率为 \(\ell\)。组合损失为
\[ L=\sum_{j=1}^m E_jD_j\ell. \]
其中 \(D_j=1\) 表示违约，\(D_j=0\) 表示未违约。若已知或已估计每个企业的违约概率 \(p_j\)，可以通过模拟 \(D_j\sim \operatorname{Bernoulli}(p_j)\) 得到组合损失分布。下面数据是模拟的，不对应任何真实金融机构。
set.seed(2028)
m <- 300
R <- 5000
risk_score <- rnorm(m)
exposure <- runif(m, min = 0.5, max = 3.0)
default_prob <- plogis(-3 + 0.8 * risk_score)
lgd <- 0.45
loss_threshold <- 22
loss <- numeric(R)
default_count <- numeric(R)
for (r in seq_len(R)) {
default <- rbinom(m, size = 1, prob = default_prob)
loss[r] <- sum(exposure * default * lgd)
default_count[r] <- sum(default)
}
c(expected_loss = mean(loss),
mc_se_expected_loss = sd(loss) / sqrt(R),
sd_loss = sd(loss),
VaR_95 = unname(quantile(loss, 0.95)),
prob_loss_gt_22 = mean(loss > loss_threshold),
mean_default_count = mean(default_count))
#> expected_loss mc_se_expected_loss sd_loss VaR_95
#> 15.06973 0.04976 3.51850 21.05371
#> prob_loss_gt_22 mean_default_count
#> 0.02980 19.36920这里 expected_loss 是平均损失，VaR_95 是模拟损失分布的 95% 分位数，prob_loss_gt_22 是损失超过 22 个单位的概率。与平均补贴支出不同，风险管理常常更关心尾部，而尾部概率和分位数通常比均值更难稳定估计。因此在报告这类结果时，模拟次数、随机种子和尾部估计的稳定性尤其重要。
10.6 重要性抽样
有时我们关心的是罕见事件，例如高收入尾部概率、极端亏损概率或违约概率。如果直接从原分布抽样，罕见事件出现次数很少，估计会很不稳定。重要性抽样通过改变抽样分布来提高效率。
设目标积分为
\[ \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)\) 贡献大的区域给出足够概率。如果 \(q(x)\) 在某些重要区域太小，权重就会非常大，估计可能由少数样本支配。
下面估计收入超过 20 万元的尾部概率。直接模拟可能需要很大样本，重要性抽样则从均值更高的对数正态分布中抽样，再用权重修正。
set.seed(3)
N <- 20000
threshold <- 20
# 直接 Monte Carlo
y_direct <- rlnorm(N, meanlog = log(8), sdlog = 0.6)
z_direct <- as.numeric(y_direct > threshold)
p_direct <- mean(z_direct)
se_direct <- sd(z_direct) / sqrt(N)
# 重要性抽样：从更偏向高收入区域的分布中抽样
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)
z_imp <- as.numeric(y_imp > threshold) * w
p_imp <- mean(z_imp)
se_imp <- sd(z_imp) / sqrt(N)
ess <- sum(w)^2 / sum(w^2)
c(direct = p_direct,
direct_se = se_direct,
importance = p_imp,
importance_se = se_imp,
max_weight = max(w),
effective_sample_size = ess)
#> direct direct_se importance
#> 6.455e-02 1.738e-03 6.396e-02
#> importance_se max_weight effective_sample_size
#> 9.438e-04 1.078e+01 1.274e+04重要性抽样不是简单地“换一个分布抽样”。如果 \(q(x)\) 选得不好，权重可能极端不稳定，估计反而变差。因此，实际使用时要检查权重分布，例如最大权重、权重分位数和有效样本量。上面代码中的
\[ \operatorname{ESS}=\frac{\left(\sum_{i=1}^N w_i\right)^2}{\sum_{i=1}^N w_i^2} \]
称为有效样本量的一个常用近似。它反映加权样本大约相当于多少个均匀权重样本。若 ESS 远小于 \(N\)，说明少数样本权重过大，应重新考虑建议分布 \(q(x)\)。
10.7 方差缩减的基本思想
Monte Carlo 方法的成本主要来自随机误差。除了增加 \(N\)，还可以通过更聪明的抽样降低方差。方差缩减方法的目标不是改变估计对象，而是在估计同一个 \(\mu\) 时，让估计量更稳定。
10.7.1 对偶变量法
最简单的例子是对偶变量法。若要估计
\[ E\{h(U)\},\quad U\sim U(0,1), \]
可以同时使用 \(U\) 和 \(1-U\)：
\[ \hat\mu_{\text{anti}}=\frac{1}{N}\sum_{i=1}^N \frac{h(U_i)+h(1-U_i)}{2}. \]
当 \(h\) 是单调函数时，\(h(U_i)\) 和 \(h(1-U_i)\) 往往负相关，从而降低方差。
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),
ratio = sd(anti) / sd(plain))
#> plain_sd antithetic_sd ratio
#> 0.034165 0.004358 0.127542这里估计的是 \(E\{\exp(U)\}\)。ratio 小于 1 表示对偶变量法降低了重复模拟中的波动。
10.7.2 控制变量法
控制变量法利用一个期望已知、且与目标变量相关的辅助变量。设 \(Z=g(X)\) 是目标模拟值，\(H\) 是辅助变量，且 \(E(H)=\eta\) 已知。对任意常数 \(a\)，
\[ \tilde\mu=\bar Z-a(\bar H-\eta) \]
仍然估计 \(E(Z)\)。如果 \(H\) 与 \(Z\) 高度相关，选择合适的 \(a\) 可以显著降低方差。理论上的最优系数为
\[ a^\ast=\frac{\operatorname{Cov}(Z,H)}{\operatorname{Var}(H)}. \]
在平均补贴支出的例子中，补贴金额 \(Z=g(Y)\) 与收入 \(Y\) 负相关，而对数正态分布的均值
\[ E(Y)=\exp\left(m+\frac{s^2}{2}\right) \]
是已知的，其中 \(m\) 是 meanlog，\(s\) 是 sdlog。因此可以把收入本身作为控制变量。
set.seed(5)
R <- 1000
N <- 500
meanlog <- log(8)
sdlog <- 0.6
income_mean <- exp(meanlog + sdlog^2 / 2)
plain <- control <- numeric(R)
for (r in seq_len(R)) {
y <- rlnorm(N, meanlog = meanlog, sdlog = sdlog)
z <- ifelse(y < 4, 2 * (1 - y / 4), 0)
a_hat <- cov(z, y) / var(y)
plain[r] <- mean(z)
control[r] <- mean(z) - a_hat * (mean(y) - income_mean)
}
c(plain_sd = sd(plain),
control_sd = sd(control),
ratio = sd(control) / sd(plain))
#> plain_sd control_sd ratio
#> 0.008697 0.008366 0.961954控制变量法提醒我们：模拟设计可以利用统计模型中已知的结构。若某个辅助量期望已知，并且与目标量相关，就可能用它修正随机波动。
10.8 Monte Carlo 结果如何报告
一个可复现的 Monte Carlo 分析至少应说明以下内容：
- 目标量是什么，例如均值、概率、分位数、回归系数偏差或均方误差；
- 随机变量如何生成，包括分布、参数和样本量；
- 模拟次数 \(N\) 或重复次数 \(R\)；
- 随机种子和主要代码；
- Monte Carlo 标准误或其他稳定性检查；
- 若使用重要性抽样或方差缩减方法，应说明建议分布、权重诊断或辅助变量。
尤其要区分“模型不确定性”“真实抽样误差”和“Monte Carlo 误差”。本章主要处理第三类误差：在模型和目标量已经明确时，计算近似本身还有多不稳定。下一章将转向 Bootstrap、Jackknife 和置换检验，它们更多从已经观测到的数据出发，用重采样近似统计量的抽样分布。
10.9 本章小结
Monte Carlo 方法用随机样本平均值近似期望和积分。它的理论基础是大数定律和中心极限定理，实际报告中应给出模拟次数和 Monte Carlo 标准误。Monte Carlo 不仅可以计算复杂积分，也可以通过重复抽样研究估计量的偏差、方差和均方误差。重要性抽样、对偶变量法和控制变量法进一步说明，模拟设计本身也是统计计算的重要部分。
10.10 练习
- 用 Monte Carlo 方法估计标准正态分布下 \(P(X>2)\)，并与
pnorm()的结果比较。 - 在平均补贴支出例子中，把收入分布的
sdlog从 0.6 改为 0.3 和 0.9，比较平均补贴支出和贫困率。 - 在消费函数模拟中，把样本量从 50 改为 1000，观察 \(\hat\beta_1\) 的标准差如何变化。
- 用重要性抽样估计 \(P(Y>25)\)，尝试不同的建议分布均值，并比较权重的最大值和有效样本量。
- 在信用组合损失例子中，把违约损失率
lgd从 0.45 改为 0.6，观察期望损失和 95% 分位数如何变化。 - 用控制变量法估计 \(E(Y^2)\)，其中 \(Y\) 服从对数正态分布，并尝试选择 \(Y\) 作为控制变量。
- 解释为什么 Monte Carlo 误差下降速度是 \(1/\sqrt{N}\)，这对模拟次数选择意味着什么？