R语言用for循环迭代pi1抽样计算OLS与IV回归偏误问题排查

存在的错误 #

• 步长设置不符合要求：需求是pi1以0.0001为步长从-1递增到0，代码中seq(from = -1, to = 0, by = 0.001)的步长参数为0.001，仅为需求的1/10，遍历的pi1值数量大幅不足
• 回归数据集未定义：调用lm和ivreg时传入的data = data中的data对象没有提前生成，你模拟得到的d、y、z等变量均为独立向量，未打包为数据框传入回归函数，运行会直接报错
• 结果存储逻辑错误：当前beta_OLS和beta_IV每次外层循环都会覆盖上一轮pi1对应的1000次模拟结果，最终仅能得到pi1=0时的偏误，无法得到所有pi1取值下的平均偏误
• 依赖包未加载：ivreg函数来自AER扩展包，代码未提前加载该包会报函数未找到错误
• 真实值定义矛盾：DGP中y方程的d的系数b1=0，但你设置的real_value为-1/1000，两者不一致会导致偏误计算完全错误，这里按DGP设置修正真实值为0，如果你需求的真实值为其他数值可自行调整

修正后可运行代码 #

# 加载依赖包
library(AER)

# 初始化存储容器，每个pi1对应一组偏误结果
n_pi1 <- length(seq(from = -1, to = 0, by = 0.0001))
biased_OLS_all <- numeric(n_pi1)
biased_IV_all <- numeric(n_pi1)
pi1_seq <- seq(from = -1, to = 0, by = 0.0001)
index <- 1:1000
real_value <- 0 # 与DGP中b1保持一致，若有其他需求可自行修改

for(i in 1:n_pi1) {
  beta_OLS <- numeric(length(index))
  beta_IV <- numeric(length(index))
  pi1 <- pi1_seq[i]
  # 固定参数
  b0 <- 0
  b1 <- 0
  b2 <- -1/1000
  b3 <- 1/5
  n <- 2000
  
  for(k in index) {
    z <- runif(n, 0, 25)
    ov <- rnorm(n, 0, 1)
    
    d <- -1/2 + pi1 * z + 1/2 * ov + rnorm(n, 0, 1) > 0
    y <- b0 + b1 * d + b2 * z + b3 * ov + rnorm(n, 0, 1/10)
    # 生成回归用数据框
    sim_data <- data.frame(y = y, d = d, z = z)
    
    # OLS Regression
    model12 = lm(y ~ d, data = sim_data)
    beta_OLS[k] = model12$coefficients[2]
      
    # IV Regression
    model12_1 = ivreg(y ~ d | z, data = sim_data)
    beta_IV[k] = model12_1$coefficients[2]
  }
  # 计算当前pi1下的偏误
  average_OLS <- mean(beta_OLS)
  average_IV <- mean(beta_IV)
  biased_OLS_all[i] <- average_OLS - real_value 
  biased_IV_all[i] <- average_IV - real_value 
}

# 输出所有pi1下的平均偏误
mean(biased_OLS_all)
mean(biased_IV_all)

内容的提问来源于stack exchange，提问作者freddywit