置换检验：为事件研究创建多日期7天窗口虚拟变量生成函数

数据老陈

2026-8-14

问题描述

本人正在开展本科毕业论文的事件研究，需检验交互效应系数的显著性。当前核心需求与难点：

运行固定效应回归，其中交互项treatment*event对应的事件窗口日期需动态调整
需创建函数生成以不同日期为起始点的连续7天事件虚拟变量
用这些虚拟变量运行固定效应回归，提取交互项系数估计值，对比真实事件日与“伪”事件日的系数差异

此前仅找到用于检验两组均值的置换检验函数，自行编写了一段R代码，现寻求技术指导与优化方案：

set.seed(1000)

N <- 10^3-1 # Number of permutations
resultX <- numeric(N)
for(i in 1:N){ 
index1 <- sample(length(merged_df), size = 7, replace = FALSE)
df_c$eventdummy <- ifelse(merged_df$day == index1, 1, 0)

resultX[i] <- plm(ret ~ sse*eventdummy + 
                    sse + 
                    eventdummy + 
                    roa + 
                    leverage +  
                    mtb + 
                    tangibility, 
                    data = merged_df,
                    index = c("entity"), 
                    model = "within")
}

代码问题分析与优化方案

原代码核心问题

采样逻辑错误：sample(length(merged_df), size=7)抽取的是数据集行号，不是连续7天的日期；且merged_df$day == index1的判断逻辑不成立，行号和日期值不匹配。
变量赋值混乱：给df_c$eventdummy赋值却用了merged_df的日期列，数据集不统一。
结果提取错误：直接把plm回归对象赋值给resultX[i]，没有提取所需的交互项系数。
冗余代码：回归公式中sse*eventdummy已经自动包含主效应sse和eventdummy，无需单独列出。

优化后的完整代码

假设merged_df包含entity（个体标识）、day（统一格式的日期列）、ret（被解释变量）、sse（处理组标识）及控制变量roa/leverage/mtb/tangibility：

1. 封装生成连续7天事件虚拟变量的函数

# 生成以指定起始日期为起点的连续7天事件虚拟变量
create_event_dummy <- function(data, start_day) {
  # 生成连续7天的日期序列
  event_days <- seq(start_day, by = "day", length.out = 7)
  # 创建虚拟变量：日期在事件窗口内为1，否则为0
  data$eventdummy <- as.integer(data$day %in% event_days)
  return(data)
}

2. 优化置换检验循环逻辑

set.seed(1000)
N <- 999  # 置换次数（10^3-1）

# 提取所有不重复的日期作为伪事件日候选池
unique_days <- unique(merged_df$day)
# 替换为你的真实事件起始日期，排除真实事件窗口避免干扰
real_start <- as.Date("2023-01-01")
real_event_days <- seq(real_start, by = "day", length.out = 7)
candidate_days <- unique_days[!unique_days %in% real_event_days]

# 初始化结果向量，存储每次置换的交互项系数
result_coef <- numeric(N)

for(i in 1:N) {
  # 随机抽取一个伪起始日期（确保能生成完整7天窗口）
  pseudo_start <- sample(candidate_days[candidate_days <= max(unique_days) - 6], size = 1)
  # 生成事件虚拟变量
  df_temp <- create_event_dummy(merged_df, pseudo_start)
  # 运行固定效应回归（自动包含主效应，无需单独写）
  model <- plm(ret ~ sse*eventdummy + roa + leverage + mtb + tangibility,
               data = df_temp,
               index = "entity",
               model = "within")
  # 提取交互项sse:eventdummy的系数
  result_coef[i] <- coef(model)["sse:eventdummy"]
}

# 运行真实事件窗口的回归，获取真实系数
df_real <- create_event_dummy(merged_df, real_start)
model_real <- plm(ret ~ sse*eventdummy + roa + leverage + mtb + tangibility,
                  data = df_real,
                  index = "entity",
                  model = "within")
real_coef <- coef(model_real)["sse:eventdummy"]

# 计算显著性p值（单侧检验：伪系数大于等于真实系数的比例）
p_value_one_sided <- mean(result_coef >= real_coef)
# 双侧检验：伪系数绝对值大于等于真实系数绝对值的比例
p_value_two_sided <- mean(abs(result_coef) >= abs(real_coef))

# 可视化伪系数分布与真实系数对比
hist(result_coef, main = "伪事件日交互项系数分布", xlab = "系数值")
abline(v = real_coef, col = "red", lwd = 2)

3. 额外效率与稳健性建议

并行加速：若数据集较大，可用furrr包替换for循环实现并行计算，大幅提升速度。
数据预处理：提前将merged_df转换为pdata.frame格式（pdata.frame(merged_df, index = c("entity", "day"))），避免plm重复处理数据。
边界校验：抽取伪起始日期时，过滤掉无法生成完整7天窗口的日期（如数据最后6天），避免缺失值。

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