R plm против Stata reghdfe

Question

Я нахожу немного другие результаты при оценке модели данных панели в Stata (используя команду , предоставленную сообществом *1002* reghdfe) и R.

Stata:

cls
webuse nlswork, clear
xtset idcode year
reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(year)  cluster(idcode)

R:

## import data
library(foreign)   
df = read_dta("http://www.stata-press.com/data/r14/nlswork.dta")

## estimate the model
model5 = plm( ln_wage ~   grade + age + ttl_exp + tenure+  not_smsa  + south + as.factor(year), data=df, index=c('idcode', 'year'), model="random")
summary(model5)[1:7,1:4]  #  <- this gives unclustered errors
coeftest(model5, vcov=vcovHC(model5,type="HC0",cluster="group"))[1:7,1:4] # <- this gives clustered errors

Я бы ожидал, что такие же коэффициенты (стандартные ошибки все еще нуждаются в коррекции степени свободы, я думаю).Чего мне не хватает?

Answer 1

После небольшой настройки я обнаружил, что пакет R's plm может использовать несколько фиксированных эффектов (по крайней мере, на обоих уровнях индекса)

## estimate the model
model5 = plm( ln_wage ~   grade + age + ttl_exp + tenure+    not_smsa  + south + as.factor(year), data=df, index=c('idcode',  'year'), model="with", effect="time")
summary(model5)[1:7,1:4]  #  <- this gives unclustered errors
coeftest(model5, vcov=vcovHC(model5,type="HC0",cluster="group"))   [1:7,1:4] # <- this gives clustered errors

Вышеуказанное равно фиксированным по времени эффектам и численно напоминает Statas reghdfe команда

 reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(year)  cluster(idcode)

Точно так же, если вы хотите оба фиксированных эффекта, где в Stata вы должны:

 reghdfe ln_w grade age ttl_exp tenure not_smsa south, abs(idcode year)  cluster(idcode) 

в R, вы можете использовать:

model5 = plm( ln_wage ~   grade + age + ttl_exp + tenure+    not_smsa  + south + as.factor(year), data=df, index=c('idcode',  'year'), model="with", effect="twoways")
summary(model5)  #  <- this gives unclustered errors
coeftest(model5, vcov=vcovHC(model5,type="HC0",cluster="group"))    # <- this gives clustered errors