Кластерные стандартные ошибки в R с использованием plm (с фиксированными эффектами)

Question

Я пытаюсь запустить регрессию в R plm пакет с фиксированными эффектами и model = 'within', сгруппировав стандартные ошибки. С использованием Cigar набор данных из plm, Я бегу:

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales + factor(state), model = 'within', data = Cigar)
coeftest(model, vcovHC(model, type = 'HC0', cluster = 'group'))

  Estimate Std. Error t value Pr(>|t|)    
sales  -1.21956    0.21136 -5.7701 9.84e-09

Это (немного) отличается от того, что я получу, используя Stata (записав файл Cigar как.dta):

use cigar

xtset state year

xtreg price sales, fe vce(cluster state)


price   Coef.   Std. Err.   t   P>t [95% Conf.  Interval]

sales   -1.219563   .2137726    -5.70   0.000   -1.650124   -.7890033

А именно, стандартная ошибка и статистика Т отличаются. Я попытался перезапустить код R с разными "типами", но ни один из них не дал тот же результат, что и Stata. Я что-то пропустил?

Answer 1

Stata использует коррекцию конечной выборки, чтобы уменьшить смещение ошибок вниз из-за конечного числа кластеров. Это мультипликативный фактор на дисперсионно-ковариационной матрице, $c=\frac{G}{G-1} \cdot \frac{N-1}{NK}$, где G - количество групп, N - число число наблюдений, а K - количество параметров. Я думаю coeftest использует только $c'=\frac{N-1}{NK}$, так как, если я масштабирую стандартную ошибку R на квадрат первого члена в c, я получаю что-то очень похожее на стандартную ошибку Stata:

display 0.21136*(46/(46-1))^(.5)
.21369554

Вот как я бы повторил то, что Stata делает в R:

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales, model = 'within', data = Cigar)
G <- length(unique(Cigar$state))
c <- G/(G - 1)
coeftest(model,c * vcovHC(model, type = "HC1", cluster = "group"))

Это дает:

t test of coefficients:

       Estimate Std. Error  t value   Pr(>|t|)    
sales -1.219563   0.213773 -5.70496 1.4319e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

что согласуется с ошибкой Stata 0,2137726 и t-stat -5,70.

Этот код, вероятно, не идеален, поскольку число состояний в данных может отличаться от количества состояний в регрессии, но мне лень выяснять, как получить правильное количество панелей.

Answer 2

Stata использует специальную коррекцию малых выборок, которая была реализована в plm 1,5.

Попробуй это:

require(plm)
require(lmtest)
data(Cigar)
model <- plm(price ~ sales + factor(state), model = 'within', data = Cigar)
coeftest(model, function(x) vcovHC(x, type = 'sss'))

Который даст:

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
sales  -1.2196     0.2137  -5.707 1.415e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Это дает ту же оценку SE до 3 цифр:

x <- coeftest(model, function(x) vcovHC(x, type = 'sss'))
x[ , "Std. Error"]
## [1] 0.2136951