В панельных данных по каждому объекту у нас есть данные по нескольким переменным, меняющимся во времени.
Загружаем нужные пакеты:
library("knitr") # используем Rmarkdown
library("texreg") # вывод результатов регрессии в тех и html
library("tidyverse") # подключаем ggplot2 (графики), dplyr, tidyr, etc (манипуляции с данными)
library("plm") # работа с панельками
# install.packages("plm") # может быть нужно установить пакеты...Берём встроенный в пакет plm набор данных Grunfeld
data(Grunfeld)
head(Grunfeld)## firm year inv value capital
## 1 1 1935 317.6 3078.5 2.8
## 2 1 1936 391.8 4661.7 52.6
## 3 1 1937 410.6 5387.1 156.9
## 4 1 1938 257.7 2792.2 209.2
## 5 1 1939 330.8 4313.2 203.4
## 6 1 1940 461.2 4643.9 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 4900.9 402.2
## 13 568.9 3526.5 761.5
## 14 529.2 3254.7 922.4
## 15 555.1 3700.2 1020.1
## 16 642.9 3755.6 1099.0h[12:16, 3:5]## inv value capital
## 1-1946 688.1 4900.9 402.2
## 1-1947 568.9 3526.5 761.5
## 1-1948 529.2 3254.7 922.4
## 1-1949 555.1 3700.2 1020.1
## 1-1950 642.9 3755.6 1099.0Добавим лагированные инвестиции
h$linv <- stats::lag(h$inv)
h[c(1:3, 21:23), ]## firm year inv value capital linv
## 1-1935 1 1935 317.6 3078.5 2.8 NA
## 1-1936 1 1936 391.8 4661.7 52.6 317.6
## 1-1937 1 1937 410.6 5387.1 156.9 391.8
## 2-1935 2 1935 209.9 1362.4 53.8 NA
## 2-1936 2 1936 355.3 1807.1 50.5 209.9
## 2-1937 2 1937 469.9 2676.3 118.1 355.3Тут есть опасный момент! В R куча пакетов определяют функцию lag. Нам нужна версия команды из пакета stats, которая берёт лаг отдельно для каждого объекта, а не просто сдвигает весь вектор значений на одно наблюдение. Например, lag из dplyr не сработает!
Иллюстрируем опасный момент:
h_bad <- mutate(h, linv = dplyr::lag(inv))
h_bad[18:22, ] ## firm year inv value capital linv
## 18 1 1952 891.2 4924.9 1430.5 755.9
## 19 1 1953 1304.4 6241.7 1777.3 891.2
## 20 1 1954 1486.7 5593.6 2226.3 1304.4
## 21 2 1935 209.9 1362.4 53.8 1486.7
## 22 2 1936 355.3 1807.1 50.5 209.9С какого-то перепугу на панелях плохо сработает и классический dplyr:
h_bad_2 <- mutate(h, linv = stats::lag(inv))
head(h_bad_2)## firm year inv value capital linv
## 1 1 1935 317.6 3078.5 2.8 317.6
## 2 1 1936 391.8 4661.7 52.6 391.8
## 3 1 1937 410.6 5387.1 156.9 410.6
## 4 1 1938 257.7 2792.2 209.2 257.7
## 5 1 1939 330.8 4313.2 203.4 330.8
## 6 1 1940 461.2 4643.9 207.2 461.2Оценим три модели:
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")Все три модели можно записать в виде: \[ y_{it}=\alpha + x'_{it}\beta + z'_i \gamma + c_i + u_{it} \]
Здесь \(z_i\) — вектор характеристик, не меняющихся во времени, а \(c_i\) и \(u_{it}\) — случайные составляющие, \(E(c_i)=0\), \(E(u_{it})=0\).
В модели со случайными эффектами (Random Effects, RE) предполагается, что \(E(c_i | z_i, X_i) = 0\).
В модели с фиксированными эффектами (Fixed Effects, FE) допускается, что \(E(c_i| X_i)\) зависит от \(X_i\). Модель с фиксированными эффектами не позволяет оценить \(\alpha\) и \(\gamma\).
В сквозной регресси (pooling) предполагается, что \(c_i=0\).
Посмотрим отчёт по модели со случайным эффектами:
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.46 52.77 0.282
## individual 7089.80 84.20 0.718
## theta: 0.8612
##
## 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.834415 28.898935 -2.0013 0.04674 *
## capital 0.308113 0.017180 17.9339 < 2e-16 ***
## value 0.109781 0.010493 10.4627 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2381400
## Residual Sum of Squares: 548900
## R-Squared: 0.7695
## Adj. R-Squared: 0.75796
## F-statistic: 328.837 on 2 and 197 DF, p-value: < 2.22e-16Все три модели рядом:
htmlreg(list(m.pooled, m.re, m.fe),
custom.model.names = c("Pooling", "RE", "FE"))| Pooling | RE | FE | ||
|---|---|---|---|---|
| (Intercept) | -42.71*** | -57.83* | ||
| (9.51) | (28.90) | |||
| capital | 0.23*** | 0.31*** | 0.31*** | |
| (0.03) | (0.02) | (0.02) | ||
| value | 0.12*** | 0.11*** | 0.11*** | |
| (0.01) | (0.01) | (0.01) | ||
| R2 | 0.81 | 0.77 | 0.77 | |
| Adj. R2 | 0.80 | 0.76 | 0.72 | |
| Num. obs. | 200 | 200 | 200 | |
| p < 0.001, p < 0.01, p < 0.05 | ||||
На самом деле модель сквозной регрессии и модель с фиксированными эффектами оцениваются обычным МНК. Эти результаты можно воспроизвести без пакета 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.7143694 0.2306785 0.1155622coefficients(m.pooled.lm)## (Intercept) capital value
## -42.7143694 0.2306785 0.1155622Аналогично совпадут и коэффициенты в моделях с фиксированными эффектами.
Проведём три теста на сравнение моделей.
pFtest(m.fe, m.pooled)##
## F test for individual effects
##
## data: inv ~ capital + value
## F = 49.177, 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.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistentЗдесь нулевая гипотеза о состоятельности коэффициентов в обеих моделях (FE и RE) отвергается в пользу гипотезы о том, что в RE модели коэффициенты несостоятельны.
Кстати, этот тест может дать отрицательное значение \(\chi^2\) статистики, т.к. хи-квадрат распределение она имеет только асимптотическое и только при верной \(H_0\). А если наблюдений мало, или если \(H_a\) верна, то можно получить отрицательные числа. Можно по приколу заметить, что модуль статистики Хаусмана также асимптотически распределён по хи-квадрату.
plmtest(m.re, type = "bp")##
## Lagrange Multiplier Test - (Breusch-Pagan)
##
## data: inv ~ capital + value
## chisq = 798.16, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effectsЗдесь нулевая гипотеза о верности сквозной регресси, то есть о том, что \(c_i\) тождественно равны \(0\), отвергается.
Что называют остатками RE модели в R и в Stata
Оценки по методу FE невозможно получить оценки коэффициентов при регрессорах, не изменяющихся во времени. Здесь для глупого примера включим номер фирмы :) А R его благополучно выкинет как неидентифицируемый.
m_fe_bad <- plm(inv ~ capital + firm, data = h, model = "within")
summary(m_fe_bad)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = inv ~ capital + firm, data = h, model = "within")
##
## Balanced Panel: n=10, T=20, N=200
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -191.000 -20.800 -0.459 21.400 294.000
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## capital 0.370750 0.019368 19.143 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2244400
## Residual Sum of Squares: 763680
## R-Squared: 0.65973
## Adj. R-Squared: 0.62345
## F-statistic: 366.446 on 1 and 189 DF, p-value: < 2.22e-16m_re_negative <- plm(inv ~ capital + value + linv, data = h, model = "random")
summary(m_re_negative)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = inv ~ capital + value + linv, data = h, model = "random")
##
## Balanced Panel: n=10, T=19, N=190
##
## Effects:## Warning in sqrt(sigma2): NaNs produced## var std.dev share
## idiosyncratic 1587.91 39.85 1.013
## individual -20.70 NA -0.013
## theta: -0.1529
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -200.00 -13.30 3.52 15.60 285.00
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -15.1075451 4.6776513 -3.2297 0.001465 **
## capital 0.0463882 0.0144357 3.2134 0.001546 **
## value 0.0237762 0.0045178 5.2628 3.884e-07 ***
## linv 0.9007120 0.0344028 26.1813 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 11543000
## Residual Sum of Squares: 451970
## R-Squared: 0.96085
## Adj. R-Squared: 0.94062
## F-statistic: 1521.48 on 3 and 186 DF, p-value: < 2.22e-16Скорее всего это говорит о том, что имеющиеся данные плохо описываются моделью RE. Вероятнее всего в данных другая структура корреляции ошибок. Например, там присутствует автокорреляция первого порядка…
Почиташки:
plm