Анализ панельных данных. Модели с фиксированными эффектами и модель со случайными эффектами
Будем знать панельные данные “в лицо”
Для начала посмотрим на то, как панельные данные могут быть представлены в массиве. Возьмем “кусочек” данных (включает только часть стран исходного массива), с которыми будем работать далее на практическом занятии. Источник данных – проект Всемирного банка Worldwide Governance Indicators. В текущем массиве представлены данные по политической стабильности и контролю коррупции за 2006 – 2015 гг. для Австралии, Великобритании, Канады, Новой Зеландии и США.
pol_stab – Политическая стабильность и отсутствие насилия/терроризма. Отражает склонность к политической нестабильности и/или политически мотивированному насилию, включая терроризм. Значения показателя варьируются от −2.5 до 2.5, причем более высокие значения означают более высокий уровень политической стабильности.
con_cor – Показатель контроля коррупции. Отражает восприятие степени использования государственной власти для личной выгоды, включая как мелкие, так и крупные формы коррупции. Показатель изменяется от −2.5 до 2.5, причем более высокие значения соответствуют более низкому уровню коррупции (т.е. более высокому контролю).
Посмотрим, как представлены данные для первой страны в массиве – Австралии. В качестве наблюдения теперь можно рассматривать не пространственную единицу, а комбинацию: “страна – год”:
## # A tibble: 10 × 4
## year country pol_stab con_cor
## <dbl> <chr> <dbl> <dbl>
## 1 2006 Australia 0.935 1.96
## 2 2007 Australia 0.929 2.01
## 3 2008 Australia 0.956 2.04
## 4 2009 Australia 0.856 2.05
## 5 2010 Australia 0.889 2.03
## 6 2011 Australia 0.936 2.04
## 7 2012 Australia 0.998 1.99
## 8 2013 Australia 1.03 1.79
## 9 2014 Australia 1.03 1.85
## 10 2015 Australia 0.885 1.88Точно таким же образом организованы данные и для других стран. Так для США – последней страны в массиве – видим аналогичную структуру:
## # A tibble: 10 × 4
## year country pol_stab con_cor
## <dbl> <chr> <dbl> <dbl>
## 1 2006 USA 0.494 1.35
## 2 2007 USA 0.376 1.39
## 3 2008 USA 0.585 1.45
## 4 2009 USA 0.448 1.29
## 5 2010 USA 0.438 1.27
## 6 2011 USA 0.591 1.27
## 7 2012 USA 0.632 1.41
## 8 2013 USA 0.643 1.31
## 9 2014 USA 0.582 1.38
## 10 2015 USA 0.678 1.40В общем виде можно использовать следующее обозначение: \(x_{it}\), где \(x\) – сам показатель, \(i\) – субиндекс, отвечающий за пространственную единицу, \(t\) – субиндекс, отвечающий за временной период. Далее будем также обозначать как N - общее количество пространственных единиц, а T – общее количество временных периодов.
А что если, как уже умеем, применить pooled model?
Такие данные отличаются от тех данных, которые мы рассматривали в рамках предыдущих курсов по регрессионному анализу, тем, что наблюдения в данном случае не являются независимыми. Условно массив можно разделить на N подгрупп – N стран. Внутри каждой такой подгруппы данные связаны: наблюдается зависимость между значениями показателей в разные временные периоды. Кстати, аналогично мы могли бы разделить массив на T подгрупп, выделенных на основе количества временных периодов, и посмотреть, а как же связаны данные по странам, собранные за один и тот же год.
Посмотрим, можно ли применить классическую линейную регрессионную модель без всяких поправок (pooled model – то есть, будем рассматривать массив как единый, закрывая глаза на то, что есть и пространственная, и временная перспектива) – ту модель, которую мы рассматривали ранее – для анализа таких данных. Применить-то можем, но то, что получится, нас не порадует. В результате такого подхода “в лоб” получим вот такие результаты, картинку ниже можно было бы считать, как то, что с ростом контроля коррупции увеличивается и показатель политической стабильности.
Однако если учесть разные подгруппы наблюдений, выделенные по пространственным единицам (разным цветом на графике обозначим разные страны), то увидим, что не все так однозначно. Картинка же выше “съела” имеющиеся различия.
Таким образом, можно отметить, что необдуманное применение pooled model может привести к смещению в результатах. Кроме того, возникнут проблемы со значимостью результатов. Так как массив воспринимается как единый, без разграничения на подгруппы, N велико, а значит стандартные ошибки малы, и вследствие этого мы будем с большей вероятностью интерпретировать результаты как значимые, когда это не так.
Что же делать?
В качестве одной из возможных альтернатив для анализа панельных данных выступает модель с фиксированными эффектами (FE-model). Мы исходим из простой предпосылки о том, что стартовые условия разные, а взаимосвязь рассматриваемых показателей одинакова во всех странах.
1. Далекая от элегантной LSDV-модель
Для начала запишем спецификацию FE-модели в форме LSDV (least-squares dummy-variable model, то есть, первая часть названия в явном виде указывает на то, что модель будет оцениваться с помощью хорошо Вам знакомого МНК, а вторая часть – на то, что в модели стоит ожидать набор дамми-переменных)
\[y_{it} = b_{0} + \gamma_{1}*D_{1i} + ... \gamma_{n-1}*D_{(n-1)i} + b_{1}*x_{it} + e_{it}\],
где \(D_{i}\) – дамми для i-ой пространственной единицы (принимает значение 1 для i-ой пространственной единицы, 0 – для всех остальных). Обратите внимание на то, что таких дамми в модели \(N-1\), так как базовая категория (та страна, с которой мы будем сравнивать все остальные) вынесена в константу. При такой спецификации все \(N\) дамми в модель включить не получится по причине строгой мультиколлинеарности (dummy-variable trap).
\(\hat{b_{0}}\) – чему в среднем равно значение зависимой переменной в базовой категории при равенстве предикторов 0;
\(\hat{\gamma_{i}}\) – на сколько в среднем отклоняется значение зависимой переменной в i-ой пространственной единице в отличие от базовой категории при прочих равных
Для понимания изобразим схематично на картинке (пока без привязки к предыдущему массиву, оставим его до практического занятия). Пусть первая страна, выделенная красным цветом, выступает базовой категорией. Константа для нее равна \(\hat{b_{0}}\), а вот для второй “зеленой” страны значение зависимой переменной в среднем выше на \(\hat{\gamma_{1}}\) при прочих равных; для третьей страны, обозначенной синим цветом – выше на \(\hat{\gamma_{2}}\) по сравнению с первой “красной” страной – базовой категорией – при прочих равных.
2. Схлопываем модель с помощью внутригруппового преобразования
Согласитесь, что такая спецификация модели довольно громоздкая. К примеру, если бы у нас было 30 стран, то 29 параметров (коэффициентов) нам пришлось бы оценивать только для дамми-переменных. Если Вы оцените такую модель, мера \(R^2\) “вздуется” из-за такого обилия параметров, но это отнюдь не означает, что модель прекрасна и на ней стоит остановиться. Не будем тащить за собой “хвост” параметров и представим модель в более экономном виде, для этого воспользуемся внутригрупповым преобразованием. Для этого для каждой подгруппы (страна) рассчитаем среднее (как по зависимой переменной, так и по предикторам), а затем заменим исходные переменные на ее отклонение от среднего группового. Посредством такого центрирования избавляемся от константы и шлейфа параметров для дамми.
\[y_{it} - \overline{y_i.} = b_{1}*(x_{it} - \overline{x_i.}) + (e_{it}-\overline{e_i.})\]
3. Откуда берется оценка коэффициента при предикторе в FE-модели на пространственные единицы?
Начнем с несколько неожиданного вопроса, а можно ли получить оценку коэффициента при предикторе в FE-модели на основе результатов оценивания моделей по подгруппам? На первый взгляд, вопрос кажется праздным: мы же можем получить оценку коэффициентов автоматически с помощью Python, вряд ли кто-то будет в стиле “Очумелых ручек” мастерить оценку на основе результатов по подгруппам. Однако понимание того, откуда берется оценка коэффициента, позволит нам продвинуться в интерпретации результатов.
Мы по-прежнему работаем с FE-моделью, запишем ее, к примеру, в виде LSDV. Для простоты представим, что у нас один предиктор и сбалансированные данные (соответствующая процедура для случая множественной регрессии с предварительным “очищением” вариации приведена в Python-ноутбуке с практического занятия):
\[\hat{y}_{it}=\hat{b}_0+\hat\gamma_1D_{1i}+...+\hat{\gamma}_{N−1}D_{(N−1)i}+\hat{b}_1x_{it}\]
Представьте, что мы поделили общий массив данных на N кусочков, каждый из которых представляет i-ую пространственную единицу (к примеру, страну или регион). На данных каждой такой подвыборки оценили отдельную регрессионную модель, где \(y\) – по-прежнему зависимая переменная, а \(x\) – объясняющая переменная. Итого таких моделей у нас будет N штук:
\[\hat{y}_{1t}=\hat{\alpha}_{01}+\hat{\alpha}_{11}x_{1t}\]
\[⋯\]
\[\hat{y}_{Nt}=\hat{\alpha}_{0N}+\hat{\alpha}_{1N}x_{Nt}\]
Выгрузим все оценки коэффициентов при предикторе (с \(\hat{\alpha}_{11}\) до \(\hat{\alpha}_{1N}\) ). Далее посчитаем их взвешенную сумму. При этом в качестве веса будет выступать доля внутристрановой изменчивости (variation как сумма квадратов отклонений исходного значения случайной величины от своего среднего значения) предиктора i-ой страны от суммы всех внутристрановых вариаций по всем странам.
\[\hat{b}_1=\sum_{i=1}^{N}\alpha_{1i}×\dfrac{\hat{Var}(x|unit=i)}{\sum_{i=1}^{N}\hat{Var}(x|unit=i)}\] Таким образом, оценка коэффициента при предикторе в модели c FE на пространственные единицы – это взвешенная сумма оценок соответствующих коэффициентов, полученных на предварительном шаге в результате оценивания отдельных регрессионных моделей y на x на каждой из подвыборок, заданных принадлежностью к i-ой пространственной единице. Тем странам, в которых наблюдается наибольшая изменчивость по x , будет присваиваться больший вес при формировании итоговой оценки коэффициента в FE-модели. И наоборот, те страны, в которых x остается постоянным, неизменяется во времени, не будут учитываться при формировании данной оценки. Знакомство с такой процедурой взвешивания позволит понять границы интерпретации результатов оценивания FE-модели, выявить, на какие пространственные единицы интерпретацию распространять поспешно.
4. Модель со случайными эффектами: спецификация
Для моделирования разных стартовых различий при одинаковой взаимосвязи существует альтернатива -– модель со случайными эффектами (RE-модель). Спецификация выглядит следующим образом:
\[y_{it} = b_{0} + b_{1}*x_{it} + \alpha_{i} + e_{it}\]
В этой модели две составляющие ошибки: те значимые характеристики на уровне стран (неизменяющиеся во времени), которые мы пропускаем, уходят в ошибку \(\alpha_{i}\), в то время как пропущенные изменяющиеся во времени и по странам характеристики оказываются во второй составляющей ошибки – \(e_{it}\) . Так, теперь различия в стартовых условиях представлены как случайная величина, в FE-модели неоднородность задается как набор параметров. Такая спецификация RE-модели предполагает следующие допущения:
1.\(Cov(\alpha_{i}, e_{it})=0\)
2.\(Cov(x_{it}, e_{it})=0\)
3.\(Cov(x_{it}, \alpha_{i})=0\)
4.\(Cov(e_{it}, e_{is}) = 0\)
5.Также предполагаем, что \(\alpha_{i} \sim N(0; Var = \sigma^2_{\alpha})\), \(e_{it} \sim N(0; Var = \sigma^2_{e}\))
Уже на основе содержательных соображений, исходя из того, с какими данными мы работаем, можно сообразить, какая модель c FE или RE нам в большей степени подходит, какую альтернативу предпочтем. Если мы работаем с некоторой специфической выборкой (к примеру, если мы работаем с посткоммунистическими странами, такую выборку нельзя назвать случайной), то FE-модель нам подойдет больше.
Метод оценивания: (F)GLS
Гетероскедастичность, зашитая в модель
Здесь появляются две составляющие ошибки: в \(\alpha_i\) закладываются межстрановые различия (отклонения групповых средних от общего среднего), а \(e_{it}\) отражает внутригрупповую вариацию (отклонения от среднего по стране). В случайный эффект попадают значимые пропущенные страновые характеристики, при такой спецификации было бы наивно надеяться на постоянство условной вариации ошибок. Гетероскедастичность уже “зашита” в эту модель.
Что применять вместо OLS в условиях гетероскедастичности? Альтернатива
Для получения эффективных оценок вместо привычного нам OLS в качестве альтернативы будет использоваться метод оценивания GLS (generalized least squares) – обобщенный метод наименьших квадратов. Если ранее мы исходили их допущения гомоскедастичности и предполагали постоянной условную вариацию ошибку = \(\sigma^2\), то теперь предварительно для того, чтобы учесть разную условную вариацию ошибки для разных пространственных единиц, мы преобразуем формулу оценивания коэффициентов, таким образом, что каждое наблюдение предварительно делится на стандартное отклонение этого наблюдения. Такая стандартизация приводит к постоянной дисперсии ошибки, что нам и нужно было.
Пусть матрица \(\Sigma\) – ковариационная матрица ошибок. Мы находим матрицу P такую, что \(\Sigma^{-1}=P^{T}P\)
Вместо X рассматриваем PX, а вместо y используем Py. Подставляем в классическую формулу оценки \(\beta_{OLS}\): \(((PX)^{T}PX)^{−1}(PX)^{T}Py=(X^{T}P^{T}PX)^{-1}X^{T}P^{T}Py\) То есть, в результате минимизации суммы квадратов остатков преобразованной модели (то есть, с учетом предварительного шага – стандартизации) получается следующая формула для расчета оценок коэффициентов: \[\hat{\beta}_{GLS} = (X^T\Sigma^{-1}X)^{-1}X^T\Sigma^{-1}y\].
Сравните с соответствующей оценкой, полученной без применения предварительного преобразования модели: \[\hat{\beta}_{OLS} = (X^TX)^{-1}X^Ty\]
FGLS: “предательское” появление F
И здесь нужно было радоваться: в идеальном мире при применении GLS все было бы хорошо, и мы бы получили в результате оценки с желательными для нас свойствами, в том числе, и эффективные, к чему стремились изначально. Но это в идеальном мире…, если у нас есть истинная ковариационная матрица ошибок. На практике же у нас есть только ее оценка, то есть, остатки вместо ошибок \(\hat{\Sigma}\). И здесь GLS превращается в FGLS – feasible generalized least squares, то есть, GLS, реализуемый на практике, доступный для нас.
\[\hat{\beta}_{FGLS} = (X^T\hat{\Sigma}^{-1}X)^{-1}X^T\hat{\Sigma}^{-1}y\].
В таком случае при ограниченной выборке мы ничего точно про свойства оценок \(\hat{\beta}\) сказать не можем. Только асимптотически на оценки FGLS можно опираться: \(\hat{\beta}_{FGLS}\) становится эквивалентной при большом количестве наблюдений \(\hat{\beta}_{GLS}\), и по эффективности выигрывает у \(\hat{\beta}_{OLS}\).
5. FE- VS RE-: формальный тест
Мы уже говорили о содержательных основаниях для выбора подходящей модели. Кроме того, для определения согласованности оценок FE- и RE-моделей существует формальный тест Хаусмана.
Нулевая гипотеза в этом тесте соответствует допущению RE-модели: а именно, \(cor(x_{it}, \alpha_{i})=0\)
Статистика критерия устроена таким образом, чтобы учесть и разницу между оценками коэффициентов в FE- и RE-модели, и разницу в вариации этих оценок. Выглядит она следующим образом:
\[S = (\hat{\beta}_{FE} - \hat{\beta}_{RE})^T(\hat{C}_{FE}-\hat{C}_{RE})^{-1}(\hat{\beta}_{FE} - \hat{\beta}_{RE})\] Такая статистика при верной нулевой гипотезе имеет распределение \(\chi^2\) с количеством степеней свободы \(k\), где \(k\) – это количество изменяющихся во времени предикторов. Если “перевести” запись: берем квадрат разницы оценок коэффициентов FE- и RE-моделей и делим его на разницу вариаций соответствующих оценок. Получается, что статистика будет иметь высокое значение (нетипичное для нулевой гипотезы), если
- оценки коэффициентов FE- и RE-моделей несогласованны между собой, наблюдается большая разница между ними, а также
- при маленькой разнице дисперсий оценок коэффициентов.
Значимое различие между FE- и RE-моделями надо трактовать отнюдь не в пользу RE-модели, так как в условиях альтернативы (нетипичное значение статистики критерия) нарушается предпосылка RE-модели, и случайный эффект и предикторы становятся скоррелированными (привет эндогенности!), что говорит о том, что мы пропустили какие-то значимые факторы на уровне стран. Поэтому при отвержении нулевой гипотезы теста Хаусмана мы отдаем предпочтение FE-модели.
И наоборот, при типичном для \(H_0\) значении статистики критерия оценки коэффициентов FE- и RE-моделей схожи, однако наблюдаются значимые различия в оценках их вариаций. Следует помнить, что в этом случае оценки FE-модели неэффективны
6. Модель с FE на временные периоды
Запишем модель с фиксированными эффектами теперь уже не на пространственные единицы, а на временные периоды. \[\hat{y}_{it} = \hat{b}_{0} + \hat{\gamma}_{1}*D_{1t} + ... \hat{\gamma}_{T-1}*D_{(T-1)i} + \hat{b}_{1}*x_{it}\],
В такой модели \(\hat{b}_0\) показывает, чему в среднем равно значение зависимой переменной во временном периоде –- базовой категории -– при равенстве предикторов 0.
\(\hat{\gamma}\) –- на сколько в среднем отклоняется значение зависимой переменной в t-ом временном периоде в отличие от временного периода –- базовой категории –- при прочих равных.
Как мы уже говорили, дамми для временных периодов “поглощают” все неизменяющиеся в пространстве, но при этом изменяющиеся во времени характеристики. Такая модель позволит отследить, как в среднем изменяется зависимая переменная в разные временные периоды. Допустим, нас интересует, как объявление COVID-19 пандемией повлияло на потребительские расходы на услуги турагентств. Мы можем сравнить соответствующий показатель в период до объявления COVID-19 пандемией и после, оценка коэффициента при дамми на временной период покажет нам характер эффекта COVID-19 на потребительские расходы на услуги турагентств.
Если такие коэффициенты при дамми на временные периоды не представляют особого интереса, и в фокусе Вашего внимания – оценка коэффициента при предикторе, тогда модель можно также “схлопнуть”, то есть, сделать более экономной посредством внутригруппового преобразования:
\[y_{it} - \overline{y_{.t}} = b_{1}*(x_{it} - \overline{x_{.t}}) + (e_{it}-\overline{e_{.t}})\]
Для этого мы оцениваем модель, предварительно вычитая из исходных значений переменных среднее по каждому временному периоду. То есть, к примеру, зафиксировали первый временной период, рассчитали среднее по выборке стран за этот временной период, также зафиксировали второй период и рассчитали среднее только по этому временному периоду, и так далее. В итоге получаем массив, разрезанный на T подвыборок. После внутригруппового преобразования получим T подмассивов со значениями переменных в терминах отклонения от среднего по временному периоду. Внутри каждого из этих T подмассивов выявляем взаимосвязь отклика и предикторов. То есть, коэффициент при предикторе теперь показывает, как в среднем при увеличении предиктора на 1 в межстрановой перспективе(!) изменяется значение отклика при прочих равных.
Сравнение с FE на пространственные единицы
А что же показывает коэффициент при предикторе в случае, если мы работаем с моделью с FE на пространственные единицы? В такую модель мы в явном виде включаем дамми, охватывающие характеристики, изменяющиеся только между странами, но при этом неизменяющиеся во времени. Таким образом, та изменчивость, которая остается для оценивания коэффициента при предикторе(-ах) – это различия между временными периодами.
Если ту же идею передать через призму внутригруппового преобразования, то получается, что мы вычитаем из исходных значений переменных среднее, рассчитанное по каждой пространственной единице. Допустим, у нас в массиве N стран. После внутригруппового преобразования получим N подмассивов со значениями переменных в терминах отклонения от среднего по стране. В каждом подмассиве (а таких подмассивов у нас N) выявляем взаимосвязь отклика и предикторов. На выходе получаем, как в среднем при увеличении предиктора на 1 во временной перспективе(!) изменяется значение отклика при прочих равных.
Таким образом, если мы хотим, чтобы коэффициент при предикторе был результатом сравнения стран в фиксированный временной период, тогда мы обращаемся к FE-модели на временные периоды. Если же, наоборот, хотим получить коэффициент, отражающий сравнение страны в разные временные периоды, тогда наш выбор – модель с FE на пространственные единицы.
7. А если учесть все вместе?
На данном этапе вполне закономерно возникает вопрос, а почему бы в модели не учесть сразу и временные периоды, и страны? В целом такое тоже возможно. Такую twoways-модель нередко называют tricky model, и дело отнюдь не в том, что выглядит она внушительно (ведь можно точно так же, как и раньше, прибегнуть внутригрупповому преобразованию). Проблема в том, что могут возникнуть проблемы с идентификацией оценок, кроме того, несколько смазывается интерпретация коэффициента при предикторе. Оценивание этой модели равносильно оцениванию при предварительном центрировании как по стране, так и по временному периоду. Берем исходный массив, совершаем внутригрупповое преобразование как в классической FE-модели, а потом берем срез по каждому временному периоду. Внутри каждой такой \(t\)-ой подвыборки содержатся уже значения переменных в терминах отклонений от среднего по стране. Внутри каждого из этих T кусочков выявляем взаимосвязь отклика и предикторов. То есть, теперь коэффициент при предикторе показывает, как в среднем при увеличении предиктора на 1 как в межстрановой, так и временной перспективе(!) изменяется значение отклика при прочих равных. Согласитесь, интерпретация самая что ни на есть tricky.
8. А если не только стартовые условия, но и взаимосвязь различается?
До этого момента мы исходили из простого допущения, что у нас разные стартовые условия, но при этом взаимосвязь предиктора и отклика одинакова во всех рассматриваемых единицах. Мы рассмотрели, как можно задать спецификацию FE- и RE-модели для этого случая. Но нередко, такое допущение не соответствует действительности, не всегда правдоподобно, что характер взаимосвязи одинаковый во всех странах или во все временные периоды.
Адаптируем модель для того, чтобы она позволила выявить и различия во взаимосвязи, к примеру, в разных группах стран. Введем дамми на наклоны:
\[\hat{y}_{it}=\hat{b}_0+\hat{\gamma}_1D_{1i}+...+\hat{\gamma}_{N-1}D_{(N-1)i}+\hat{b}_1x_{it}+\hat{c}_1D_{1i}x_{it}+...\hat{c}_{N-1}D_{(N-1)i}x_{it}\] , где
\(\hat{b}_0\) – то, чему в среднем равен \(y_{it}\) в группе стран – базовой категории – при всех предикторах равных 0;
\(\hat{\gamma}\) – отклонение \(y_{it}\) в среднем в i-ой группе стран от \(y_{it}\) в базовой категории при \(x_{it}\) равном 0 (теперь наклоны разные, поэтому фиксируем определенную точку);
\(\hat{b}_1\) показывает в среднем взаимосвязь \(x_{it}\) и \(y_{it}\) в погруппе стран, выбранной базовой категорией
\(\hat{c}\) – насколько в среднем отличается взаимосвязь \(x_{it}\) и \(y_{it}\) в i-ой стране в отличие от базовой категории при прочих равных