R中如何实现带AR(1)扰动的固定效应分析？panelAR包已下架

问题背景 #

需要执行带AR(1)序列相关扰动的面板固定效应分析，原依赖的panelAR包已从CRAN下架，无法直接调用，且偏好使用R而非STATA的xtregar命令。原尝试代码如下：

pdata20032020_2 <- panelAR(RDItplusone ~
                         PBA
                       +PAA
                       #+CASH
                       #+COMMON_SHARES_OUTSTANDING
                       +SIC_CODE
                       +MS
                       +HHI
                       +AbsorbedSlack
                       +UnabsorbedSlack
                       +PotentialSlack
                       #+SchumacherSlack_Bromiley
                       #+CURRENT_LIABILITIESTOTAL
                       #+CURRENT_RATIO
                       #+MARKET_VALUE
                       #+NET_INCOME_DOLLAR
                       +NET_SALES_OR_REVENUES
                       #+QUICK_RATIO
                       #+SHORTTERM_INVESTMENTS
                       +TOTAL_ASSETS
                       #+CEO_TOTAL_WRDS
                       #+Analyst_TOTAL_WRDS
                       #+Analyst_CCC_standardized_avg
                       , data=pdata20032020
                       ,timeVar = 'year'
                       ,panelVar = 'company'
                       ,autoCorr = c("ar1","none", "psar1")
                       ,panelCorrMethod = c("none","phet","pcse","pwls", "parks")
                       , rhotype ="breg"
                       , bound.rho = FALSE
                       , rho.na.rm = FALSE
                       , panel.weight = c("t-1", "t")
                       , dof.correction = FALSE
                       , complete.case = FALSE
                       , seq.times = FALSE
                       , singular.ok=TRUE) 

替代方案 #

方案1：从GitHub安装旧版panelAR包 #

虽然CRAN已下架，部分开发者会在GitHub留存包的源码，可通过devtools安装：

# 安装devtools（若未安装）
install.packages("devtools")
# 从GitHub安装存档版本
devtools::install_github("cran/panelAR")

注意：GitHub存档可能未维护，若出现编译错误，建议使用下方其他方案。

方案2：使用plm包实现固定效应+AR(1)扰动 #

plm是R中主流面板数据分析包，支持通过可行广义最小二乘（FGLS）处理AR(1)序列相关：

# 安装并加载plm包
install.packages("plm")
library(plm)

# 将数据转换为面板数据格式
pdata <- pdata.frame(pdata20032020, index = c("company", "year"))

# 估计固定效应模型并处理AR(1)序列相关
ar1_fe_model <- pggls(RDItplusone ~ PBA + PAA + SIC_CODE + MS + HHI + 
                        AbsorbedSlack + UnabsorbedSlack + PotentialSlack + 
                        NET_SALES_OR_REVENUES + TOTAL_ASSETS,
                      data = pdata,
                      model = "within",
                      correlation = corAR1(form = ~ year | company))

# 查看结果
summary(ar1_fe_model)

• corAR1(form = ~ year | company)指定每个个体（company）的残差遵循AR(1)过程，时间维度为year。
• model="within"对应固定效应模型，与原panelAR的固定效应设定一致。

方案3：使用lfe包+手动拟合AR(1) #

若需要更灵活的固定效应设定（如多维度固定效应），可先用lfe估计固定效应，再手动拟合AR(1)并执行FGLS：

# 安装并加载lfe包
install.packages("lfe")
library(lfe)

# 第一步：估计固定效应模型，提取残差
fe_model <- felm(RDItplusone ~ PBA + PAA + SIC_CODE + MS + HHI + 
                   AbsorbedSlack + UnabsorbedSlack + PotentialSlack + 
                   NET_SALES_OR_REVENUES + TOTAL_ASSETS | company,
                 data = pdata20032020)
residuals <- residuals(fe_model)

# 第二步：对每个个体的残差拟合AR(1)，估计rho值
install.packages("dplyr")
library(dplyr)
rho_estimates <- pdata20032020 %>%
  group_by(company) %>%
  mutate(lag_resid = lag(residuals)) %>%
  na.omit() %>%
  summarise(rho = coef(lm(residuals ~ lag_resid - 1))[1])

# 第三步：基于估计的rho执行FGLS
ar1_fe_model_fgls <- pggls(RDItplusone ~ PBA + PAA + SIC_CODE + MS + HHI + 
                             AbsorbedSlack + UnabsorbedSlack + PotentialSlack + 
                             NET_SALES_OR_REVENUES + TOTAL_ASSETS,
                           data = pdata,
                           model = "within",
                           correlation = corAR1(value = rho_estimates$rho, form = ~ year | company))

summary(ar1_fe_model_fgls)

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