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

Я пытаюсь запустить регрессию в пакете 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

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

2 ответа

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.

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


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 до трех цифр:

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

licensed under cc by-sa 3.0 with attribution.