面板数据循环回归与F统计量最优起始年份求解

1. 模拟示例面板数据 #

先构建符合需求的简化示例数据，包含核心字段：

library(dplyr)
library(plm)

set.seed(123)
df <- expand.grid(country = c("USA", "CHN", "DEU"), year = 2010:2020) %>%
  mutate(potential = ifelse(year %in% c(2012,2015,2018), 1, 0),
         y = rnorm(n(), 50, 10),
         x1 = rnorm(n(), 10, 2),
         x2 = rnorm(n(), 20, 3))

2. 回归与F统计量计算函数 #

针对每个国家的potential=1年份执行指定回归，提取F统计量：

calc_fstat <- function(sub_df) {
  sub_df_filtered <- sub_df %>% filter(potential == 1)
  if(nrow(sub_df_filtered) == 0) return(tibble(year = NA, fstat = NA))
  
  # 按实际需求调整回归模型（此处为固定效应面板回归）
  model <- plm(y ~ x1 + x2, data = sub_df_filtered, index = c("country", "year"), model = "within")
  f_stat <- summary(model)$fstatistic[1]
  
  sub_df_filtered %>%
    mutate(fstat = f_stat) %>%
    select(country, year, fstat)
}

3. 分组计算并标记最高F值年份 #

按国家分组计算F值，筛选每组F值最高的年份并标记：

# 分组计算F统计量
fstat_results <- df %>%
  group_by(country) %>%
  do(calc_fstat(.)) %>%
  ungroup()

# 标记每组F值最高的年份
fstat_results <- fstat_results %>%
  group_by(country) %>%
  mutate(max = ifelse(fstat == max(fstat, na.rm = TRUE), 1, 0)) %>%
  ungroup()

4. 合并回原数据集 #

将计算结果整合到原始数据，非potential=1的年份填充NA：

final_df <- df %>%
  left_join(fstat_results, by = c("country", "year")) %>%
  mutate(
    fstat = ifelse(potential == 1, fstat, NA),
    max = ifelse(potential == 1, max, NA)
  )

# 查看输出结果
head(final_df, 15)

Python版本替代方案 #

若使用Python处理，核心逻辑一致，用linearmodels做面板回归：

import pandas as pd
import numpy as np
from linearmodels import PanelOLS

# 模拟数据
np.random.seed(123)
countries = ["USA", "CHN", "DEU"]
years = range(2010, 2021)
df = pd.MultiIndex.from_product([countries, years], names=["country", "year"]).to_frame(index=False)
df["potential"] = df["year"].isin([2012,2015,2018]).astype(int)
df["y"] = np.random.normal(50, 10, len(df))
df["x1"] = np.random.normal(10, 2, len(df))
df["x2"] = np.random.normal(20, 3, len(df))

# 定义计算函数
def calc_fstat(group):
    filtered = group[group["potential"] == 1]
    if len(filtered) == 0:
        return pd.DataFrame({"year": [], "fstat": []})
    filtered = filtered.set_index(["country", "year"])
    model = PanelOLS.from_formula("y ~ x1 + x2 + EntityEffects", data=filtered)
    results = model.fit(cov_type="clustered", cluster_entity=True)
    f_stat = results.f_statistic.stat
    return filtered.reset_index()[["country", "year"]].assign(fstat=f_stat)

# 分组计算+标记最高F值
fstat_results = df.groupby("country").apply(calc_fstat).reset_index(drop=True)
fstat_results["max"] = fstat_results.groupby("country")["fstat"].transform(lambda x: x == x.max()).astype(int)

# 合并回原数据
final_df = df.merge(fstat_results, on=["country", "year"], how="left")
final_df[["fstat", "max"]] = final_df[["fstat", "max"]].where(df["potential"] == 1, np.nan)

关键提示 #

• 根据实际需求调整回归公式（如加入时间效应、工具变量等），F统计量提取方式需对应修改
• 若存在多组竞争起始年份（而非仅按国家分组），只需在分组时加入周期变量即可

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