Стандартные ошибки Дрисколла и Краая в регрессии FE: воспроизведение вывода Stata xtscc в R

Я пытаюсь воспроизвести результаты, предоставленные командой Stata xtscc в R с пакетом plm, но у меня возникают проблемы с тем, чтобы увидеть те же стандартные ошибки, которые я использую в Stata для репликации набора данных из пакета plm.

# code to obtain dataset
library(lmtest)
library(car)
library(tidyverse)
data("Produc", package="plm")
write.dta(Produc,"test.dta")

Моя цель состоит в том, чтобы выполнить двухстороннюю оценку модели панели с фиксированным эффектом со стандартными ошибками Дрисколла и Краая. Процедура в Stata следующая

use "test.dta", clear \\ to import data
** i declare the panel 
xtset state year
* create the dummies for the time fixed effects
quietly tab year, gen(yeardum)
* run a two way fixed effect regression model with Driscoll and Kraay standard errors
xi: xtscc gsp pcap emp unemp yeardum*,fe 
* results are the following
                    Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
        pcap |  -.1769881    .265713    -0.67   0.515    -.7402745    .3862983
         emp |   40.61522   2.238392    18.14   0.000     35.87004     45.3604
       unemp |   23.59849   85.10647     0.28   0.785    -156.8192    204.0161

В R я использую следующую процедуру:

# I declare the panel
Produc <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)
# run a two way fixed effect model
femodel <- plm(gsp~pcap+emp+unemp, data=Produc,effect = "twoway", 
               index = c("iso3c","year"), model="within")
# compute Driscoll and Kraay standard errors using vcovSCC
coeftest(femodel, vcovSCC(femodel))

pcap  -0.17699    0.25476 -0.6947   0.4874    
emp   40.61522    2.14610 18.9252   <2e-16 ***
unemp 23.59849   81.59730  0.2892   0.7725

Хотя точечные оценки такие же, как и в Stata, стандартные ошибки отличаются.

Чтобы проверить, не использую ли я неправильную настройку небольшой выборки для стандартных ошибок, я также попытался запустить coeftest со всеми доступными настройками, но ни одна из них не дает тех же значений, что и xtscc.

library(purrr)
results <- map(c("HC0", "sss", "HC1", "HC2", "HC3", "HC4"),~coeftest(femodel, vcovSCC(femodel,type = .x)))
walk(results,print)
# none of the estimated standard errors is the same as xtscc

Кто-нибудь знает, как я могу воспроизвести результаты Stata в R?


r panel stata panel-data plm
person Alex    schedule 09.12.2020    source источник
comment
Вероятно, поскольку Stata включает перехват в свою внутреннюю модель, а plm — нет, поэтому корректировки небольших выборок немного отличаются (конечный случай; нерелевантный в больших выборках). Есть дискуссии о различиях Stata и plm в настройке небольшого примера здесь или на Stackoverflow, может быть, это поможет. Вы также можете просмотреть кодировку (или документацию) plm::within_intercept, чтобы увидеть, как рассчитывается модель FE с перехватом (документация была ссылками). Я подозреваю, что xtssc каким-то образом использует Stata xtreg. Я не знаю, использует ли xtssc ту же настройку небольшого образца, что и xtreg.   -  person Helix123    schedule 10.12.2020
comment
Спасибо за этот комментарий, проблема с реальными данными заключается в том, что, хотя коэффициенты R значимы, коэффициенты Stata - нет. поправка 0,03. Также есть ли шанс, что вы знаете, какие типы стандартных ошибок используются в Stata? Из документации xtscc непонятно..   -  person Alex    schedule 10.12.2020
comment
xtscc на самом деле является пользовательской командой. Если его документация не раскрывает использованную небольшую примерную настройку, вы можете прочитать исходный код функции или связаться с ее автором. В plm введите = sss имитирует небольшую выборку настройки Stata (отсюда и ее название).   -  person Helix123    schedule 10.12.2020
comment
Обратите внимание, что аргумент index в вашем вызове plm игнорируется, поскольку ваши данные уже являются pdata.frame. Кроме того, аргумент индекса содержит переменную, которой нет в данных (iso3c).   -  person Helix123    schedule 13.12.2020

Ответы (1)


arrow_upward
2
arrow_downward

Начиная с plm версии 2.4, его функция within_intercept(., return.model = TRUE) может возвращать полную модель внутренней модели с перехватом, как в Stata. Благодаря этому можно точно воспроизвести результат пользовательской команды Stata xtscc.

Способ xtscc, по-видимому, работает, заключается в оценке двухсторонней модели КЭ как односторонней модели КЭ + манекены для временного измерения. Итак, давайте воспроизведем это с plm:

data("Produc", package="plm")
Produc <- pdata.frame(Produc, index = c("state","year"), drop.index = FALSE)

femodel <- plm(gsp ~ pcap + emp + unemp + factor(year), data = Produc, model="within")
femodelint  <- within_intercept(femodel,  return.model = TRUE)

lmtest::coeftest(femodelint, vcov. = function(x) vcovSCC(x, type = "sss"))
#                     Estimate  Std. Error t value              Pr(>|t|)    
# (Intercept)      -6547.68816  3427.47163 -1.9104             0.0564466 .  
# pcap                -0.17699     0.26571 -0.6661             0.5055481    
# emp                 40.61522     2.23839 18.1448 < 0.00000000000000022 ***
# unemp               23.59849    85.10647  0.2773             0.7816356    
# [...]
person Helix123    schedule 12.12.2020