Загружаем нужные пакеты:
library(knitr)
opts_chunk$set(cache = TRUE)
library(ggplot2)
library(plm)
# install.packages('plm') # может быть нужно установить пакеты...
Берём встроенный в пакет plm набор данных Grunfeld
data(Grunfeld)
head(Grunfeld)
## firm year inv value capital
## 1 1 1935 317.6 3078 2.8
## 2 1 1936 391.8 4662 52.6
## 3 1 1937 410.6 5387 156.9
## 4 1 1938 257.7 2792 209.2
## 5 1 1939 330.8 4313 203.4
## 6 1 1940 461.2 4644 207.2
Указываем R, какая переменная отвечает за \( i \), какая — за \( t \).
h <- pdata.frame(Grunfeld, index = c("firm", "year"), row.names = TRUE)
Посмотрим на кусок таблички:
Grunfeld[12:16, 3:5]
## inv value capital
## 12 688.1 4901 402.2
## 13 568.9 3526 761.5
## 14 529.2 3255 922.4
## 15 555.1 3700 1020.1
## 16 642.9 3756 1099.0
h[12:16, 3:5]
## inv value capital
## 1-1946 688.1 4901 402.2
## 1-1947 568.9 3526 761.5
## 1-1948 529.2 3255 922.4
## 1-1949 555.1 3700 1020.1
## 1-1950 642.9 3756 1099.0
Добавим лагированные инвестиции
h$linv <- lag(h$inv)
head(h)
## firm year inv value capital linv
## 1-1935 1 1935 317.6 3078 2.8 NA
## 1-1936 1 1936 391.8 4662 52.6 317.6
## 1-1937 1 1937 410.6 5387 156.9 391.8
## 1-1938 1 1938 257.7 2792 209.2 410.6
## 1-1939 1 1939 330.8 4313 203.4 257.7
## 1-1940 1 1940 461.2 4644 207.2 330.8
Оценим три модели:
m.pooled <- plm(inv ~ capital + value, data = h, model = "pooling")
m.re <- plm(inv ~ capital + value, data = h, model = "random")
m.fe <- plm(inv ~ capital + value, data = h, model = "within")
Посмотрим отчёт по модели со случайным эффектами:
summary(m.re)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = inv ~ capital + value, data = h, model = "random")
##
## Balanced Panel: n=10, T=20, N=200
##
## Effects:
## var std.dev share
## idiosyncratic 2784.5 52.8 0.28
## individual 7089.8 84.2 0.72
## theta: 0.861
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -178.00 -19.70 4.69 19.50 253.00
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -57.8344 28.8989 -2.0 0.047 *
## capital 0.3081 0.0172 17.9 <2e-16 ***
## value 0.1098 0.0105 10.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2380000
## Residual Sum of Squares: 549000
## R-Squared : 0.77
## Adj. R-Squared : 0.758
## F-statistic: 328.837 on 2 and 197 DF, p-value: <2e-16
На самом деле модель сквозной регрессии и модель с фиксированными эффектами оцениваются обычным МНК. Эти результаты можно воспроизвести без пакета plm своими руками:
m.pooled.lm <- lm(inv ~ capital + value, data = h)
m.fe.lm <- lm(inv ~ capital + value + factor(firm), data = h)
Сравним коэффициенты:
coefficients(m.pooled)
## (Intercept) capital value
## -42.7144 0.2307 0.1156
coefficients(m.pooled.lm)
## (Intercept) capital value
## -42.7144 0.2307 0.1156
Аналогично совпадут и коэффициенты в моделях с фиксированными эффектами.
Проведём три теста на сравнение моделей.
pFtest(m.fe, m.pooled)
##
## F test for individual effects
##
## data: inv ~ capital + value
## F = 49.18, df1 = 9, df2 = 188, p-value < 2.2e-16
## alternative hypothesis: significant effects
phtest(m.fe, m.re)
##
## Hausman Test
##
## data: inv ~ capital + value
## chisq = 2.33, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistent
plmtest(m.re, type = "bp")
##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: inv ~ capital + value
## chisq = 798.2, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects
Типичные проблемы при работе с панельными данными.
Например,
…
…
m.re <- plm(inv ~ capital + value + linv, data = h, model = "random")
## Error: the estimated variance of the individual effect is negative
Скорее всего это говорит о том, что имеющиеся данные плохо описываются моделью RE. Вероятнее всего в данных другая структура корреляции ошибок. Например, там присутствует автокорреляция первого порядка…