工具变量双重差分法：如何用ivreg实现内生变量与外生变量交互项估计

问题背景 #

需要估计处理效应在两个分组的异质性，基准回归设定为：

$Y_i = \beta_0 + \beta_1 M_i + \beta_2 T_i + \beta_3 M_i\times T_i + u_i$
其中$T$是处理变量，$M$是分组虚拟变量，$T$因和遗漏变量相关存在内生性，已找到$T$的外生工具变量$Z$。
你此前手动实现的2SLS方法存在逻辑漏洞：交互项$M\times T$本身也是内生变量，不能直接用$T$的第一阶段拟合值乘$M$构造交互项的拟合值。正确的2SLS需要给两个内生变量$T$和$M\times T$分别配置工具变量，也就是外生的$Z$和$M\times Z$（$M$是外生虚拟变量，和外生工具$Z$的交互项依然满足外生性要求），这也是你之前手动结果和ivreg结果有偏差的核心原因。

问题解答 #

i) ivreg的正确参数写法 #

你最后尝试的写法就是符合计量理论的正确写法：

AER::ivreg(Y ~ M + Treat + M:Treat | M + Z + M:Z)

这个写法已经明确告知ivreg：

• 第二阶段的解释变量是M、Treat、M:Treat
• 所有外生变量（包括工具和外生控制）是M、Z、M:Z，刚好给两个内生变量Treat和M:Treat匹配了两个工具Z和M:Z，不存在识别不足的问题，估计结果是正确的，你此前认为的“偏差”实际上是手动实现错误导致的。
  你之前的错误写法工具数量少于内生变量数量，所以会报more regressors than instruments的警告。

ii) 支持手动指定一二阶段的2SLS估计包 #

推荐两个更灵活的常用包：

• fixest包的feols函数：语法简洁，支持多阶段IV、高维固定效应、聚类标准误等，写法示例：

  library(fixest)
  feols(Y ~ M | Treat + I(Treat*M) ~ Z + I(Z*M), data = dat)
  

  可以直接在竖线后分别指定内生变量和对应的工具变量，输出会自动调整2SLS标准误。
• lfe包的felm函数：同样支持多阶段IV和高维固定效应，适合大样本估计。

iii) 手动计算正确标准误的实现方法 #

如果不想依赖现成IV估计包，最简单的方法是使用自举法（Bootstrap），步骤如下：

1. 有放回抽取和原样本量一致的bootstrap样本
2. 在每个bootstrap样本上执行完整的2SLS估计（两个第一阶段+第二阶段）
3. 重复步骤1-2数百次，用所有bootstrap估计值的标准差作为系数的标准误
   示例代码：

library(boot)
boot_fun <- function(data, idx){
  d <- data[idx, ]
  # 两个第一阶段回归
  fs1 <- lm(Treat ~ M + Z + I(M*Z), data = d)
  fs2 <- lm(I(Treat*M) ~ M + Z + I(M*Z), data = d)
  # 第二阶段回归
  ss <- lm(Y ~ M + fitted(fs1) + fitted(fs2), data = d)
  return(coef(ss))
}
# 构造数据集
dat <- data.frame(Y, M, Treat, Z)
# 执行自举，重复500次
set.seed(123)
boot_res <- boot(dat, boot_fun, R = 500)
# 提取标准误
apply(boot_res$t, 2, sd)

如果需要直接推导计算标准误，可直接使用2SLS方差协方差矩阵公式：

# 构造X矩阵（包含常数、M、Treat、M*Treat）
X <- model.matrix(~ M + Treat + I(Treat*M), data = dat)
# 构造Z矩阵（包含常数、M、Z、M*Z）
Z_mat <- model.matrix(~ M + Z + I(M*Z), data = dat)
# 2SLS估计值
beta_2sls <- solve(t(X)%*%Z_mat%*%solve(t(Z_mat)%*%Z_mat)%*%t(Z_mat)%*%X) %*% t(X)%*%Z_mat%*%solve(t(Z_mat)%*%Z_mat)%*%t(Z_mat)%*%Y
# 计算残差
u_hat <- Y - X %*% beta_2sls
# 计算标准误
sigma2 <- sum(u_hat^2)/(nrow(X) - ncol(X))
vcov_2sls <- sigma2 * solve(t(X)%*%Z_mat%*%solve(t(Z_mat)%*%Z_mat)%*%t(Z_mat)%*%X)
sqrt(diag(vcov_2sls))

内容的提问来源于stack exchange，提问作者Igor Rigolon La-Cava Veiga