Блог

Я пытаюсь понять, почему пакеты R «plm» и «fixest» дают мне разные стандартные ошибки, когда я оцениваю панельную модель с использованием стандартных ошибок, устойчивых к гетероскедастичности («HC1») и фиксированных эффектов состояния.

У кого-нибудь есть подсказка для меня?

Вот код:

 library(AER)       # For the Fatality Dataset
library(plm)       # PLM 
library(fixest)    # Fixest
library(tidyverse) # Data Management 


data("Fatalities")

# Create new variable : fatality rate
Fatalities <- Fatalities %>% 
  mutate(fatality_rate = (fatal/pop)*10000)

# Estimate Fixed Effects model using the plm package
plm_reg <- plm(fatality_rate ~ beertax,
               data = Fatalities,
               index = c("state", "year"),
               effect = "individual")

# Print Table with adjusted standard errors
coeftest(plm_reg, vcov. = vcovHC, type = "HC1")
# Output 
>t test of coefficients:

        Estimate Std. Error t value Pr(>|t|)  
beertax -0.65587    0.28880  -2.271  0.02388 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Estimate the very same model using the fixest package
# fixest is much faster and user friendly (in my opinion) 
fixest_reg <- feols(fatality_rate ~ beertax | state , 
                    data = Fatalities,
                    vcov = "HC1",
                    panel.id = ~ state   year)

# print table 
etable(fixest_reg)

#output
>                         fixest_reg
Dependent Var.:      fatality_rate
                                  
beertax         -0.6559** (0.2033)
Fixed-Effects:  ------------------
state                          Yes
_______________ __________________
S.E. type       Heteroskedas.-rob.
Observations                   336
R2                         0.90501
Within R2                  0.04075
 

В этом примере стандартная ошибка больше при использовании plm по сравнению с самыми фиксированными результатами (то же самое верно, если state year используются фиксированные эффекты). Кто — нибудь знает причину, по которой это произошло?

На самом деле ВКОВы другие.

По plm vcovHC умолчанию используется Arellano (1987), который также учитывает последовательную корреляцию. Смотрите документацию здесь.

Если вы добавите аргумент method = "white1" , вы получите тот же тип VCOV.

Наконец, вам также необходимо изменить способ учета исправленных эффектов fixest , чтобы получить те же стандартные ошибки (подробнее о небольшой коррекции выборки см. Здесь).

Вот результаты:

 # Requesting "White" VCOV
coeftest(plm_reg, vcov. = vcovHC, type = "HC1", method = "white1")
#> 
#> t test of coefficients:
#> 
#>         Estimate Std. Error t value  Pr(>|t|)    
#> beertax -0.65587    0.18815 -3.4858 0.0005673 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
# Changing the small sample correction in fixest (discarding the fixed-effects)
etable(fixest_reg, vcov = list("hc1", hc1 ~ ssc(fixef.K = "none")), fitstat = NA)
#>                         fixest_reg          fixest_reg
#> Dependent Var.:      fatality_rate       fatality_rate
#>                                                       
#> beertax         -0.6559** (0.2033) -0.6559*** (0.1882)
#> Fixed-Effects:  ------------------ -------------------
#> state                          Yes                 Yes
#> _______________ __________________ ___________________
#> S.E. type       Heteroskedas.-rob. Heteroskedast.-rob.

# Final comparison
rbind(se(vcovHC(plm_reg, type = "HC1", method = "white1")),
      se(fixest_reg, hc1 ~ ssc(fixef.K = "none")))
#>        beertax
#> [1,] 0.1881536
#> [2,] 0.1881536