Курица или яйцо: причинность по Грэнджеру / forpes.ru

Главная
Курица или яйцо: причинность по Грэнджеру

Курица или яйцо: причинность по Грэнджеру

27.04.2017 08:24

kxx 0 330 Источник

caption

Хорошо известно, что корреляция далеко не всегда подразумевает причинность. И примеров тому великое множество. В этой статье пойдет речь о статистическом тесте, который является необходимым условием для определения причинно-следственной связи между величинами (как правило, стационарными временными рядами). Содержание:

Что такое причинность и как ее установить

aristotel

Рассуждения о причинности восходят своими корнями к древним философам. Аристотель, изучая проблему бытия и всего такого, пришел к четырем вопросам, ответив на которые, можно получить наиболее полное представление об исследуемом предмете. Мы не будем погружаться в пучины метафизики, а ограничимся более простым утверждением, что причина — это "явление, вызывающее или обусловливающее возникновение другого явления". К сожалению, даже при таком определении выявить причинность не так-то и просто — для этого может понадобиться строго контролируемый эксперимент. Но и он не всегда гарантирует установление причинно-следственных связей — чаще же лишь придает вес той или иной гипотезе.

emission

Пример из физики: нагрев катода в вакуумной лампе увеличивает эмиссию свободных электронов. Тут ясно вырисовывается зависимость: повышение температуры катода является причиной (возможно, одной из) роста количества излученных электронов. Доказать утверждение "эмиссия электронов является причиной нагрева катода" вряд ли будет просто.

Но часто проведение эксперимента по многим причинам невозможно, особенно, когда исследуются макроэкономические показатели. Поэтому определение причинности мы редуцируем до уровня простого эконометрического теста на причинность, который был предложен Грэнджером в 1969 г. Идея такова: причина X предшествует следствию Y и влияет на его будущие значения; в то же время Y не оказывает существенное влияние на будущие значения X. Для формализации сказанного запишем векторную авторегрессию (VAR):

$y_t = a_0 + a_1 \cdot y_{t-1} + \ldots + a_p \cdot y_{t-p} + b_1 \cdot x_{t-1} + \ldots + b_p \cdot x_{t-p} + \epsilon_t \\ x_t = c_0 + c_1 \cdot x_{t-1} + \ldots + c_p \cdot x_{t-p} + d_1 \cdot y_{t-1} + \ldots + d_p \cdot y_{t-p} + \tau_t$

Отвергнув гипотезу H₀¹: b₁ =… = b_p = 0, мы сможем утверждать, что X является Грэнджер-причиной Y.
Аналогично, если отвергнуть гипотезу H₀²: d₁ =… = d_p = 0, можно будет сказать, что Y является Грэнджер-причиной X.

Однако этот подход не применим к нестационарным рядам; кроме того, необходимо как-то определить количество лагов p. Тода и Ямамото в 1995 г. предложили такой алгоритм:

Определить максимальный порядок разностей рядов m (тест KPSS).
Построить VAR-модель на исходных данных и определить параметр p на основании критерия AIC/BIC.
Проверить модель VAR(p) — особенно на автокорреляцию ошибок.
Построить неограниченную модель VAR(p+m).
Выполнить тест Вальда для модели VAR(p+m) для коэффициентов первых p переменных, которые мы считаем причиной (т.е. проверить H₀).

Методология Тода-Ямамото в R

В этой статье рассматривалась зависимость между курсом рубля, стоимостью нефти и международными резервами. Можно задать закономерный вопрос: что является Грэнджер-причиной курса рубля и объемов международных резервов?

Коды

EIA/PET_RBRTE_D — спотовая цена на нефть марки Brent в $.
BOE/XUDLBK69 — стоимость 1 доллара в российских рублях .
BANKRUSSIA/RESRV — международные резервы.

Данные от Quandl

library(Quandl)
library(forecast)
library(vars)
library(aod)
oil.ts <-Quandl("EIA/PET_RBRTE_D", trim_start="2005-04-03", trim_end="2017-04-15", type="zoo", collapse="daily")
usd.ts <-Quandl("BOE/XUDLBK69", trim_start="2005-04-03", trim_end="2017-04-15", type="zoo", collapse="daily")
res.ts <- Quandl("BANKRUSSIA/RESRV", trim_start="2005-04-03", trim_end="2017-04-15", type="zoo", collapse="daily")[, 1]
rub.ts <- 1 / usd.ts

dates <- as.Date(Reduce(intersect, Map(as.character, list(index(oil.ts), index(rub.ts), index(res.ts)))))
rub.h <- window(rub.ts, dates)
oil.h <- window(oil.ts, dates)
res.h <- window(res.ts, dates)

За последнее десятилетие сложилась такая картина:

Выполним процедуру Тода-Ямамото. Сначала определим максимальный порядок разностей:

sapply(list(oil.h, rub.h, res.h), ndiffs, alpha=0.05, test=c("kpss"))
## [1] 1 1 2

Найдем количество лагов p для каждой из пар переменных, причем выбирать будем на основании BIC (критерий Шварца — SC), т.к. он дает линейную модель, которая наилучшим образом объясняет данные и не так склонна к переобучению (AIC в свою очередь дает модель, которая лучше подходит для предсказаний). Заодно проверим остатки моделей на автокорреляцию и стационарность процессов:

Выбор VAR(p): курс рубля, стоимость нефти

rub.oil <- data.frame(rub=rub.h, oil=oil.h)
VARselect(rub.oil, lag.max=20, type="both")$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2
m.rub.oil <- VAR(rub.oil , p=1, type="both")
serial.test(m.rub.oil)
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object m.rub.oil
## Chi-squared = 50.85, df = 60, p-value = 0.794
roots(m.rub.oil)
## [1] 0.9506320 0.7437097

Выбор VAR(p): курс рубля, резервы

rub.res <- data.frame(rub=rub.h, res=res.h)
VARselect(rub.res, lag.max=20, type="both")$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2
m.rub.res <- VAR(rub.res, p=2, type="both")
serial.test(m.rub.res)
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object m.rub.res
## Chi-squared = 48.86, df = 56, p-value = 0.7395
roots(m.rub.res)
## [1] 0.9335624 0.8248357 0.3126882 0.2906614

Выбор VAR(p): резервы, стоимость нефти

res.oil <- data.frame(res=res.h, oil=oil.h) 
VARselect(res.oil, lag.max=20, type="both")$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      2      2      3
m.res.oil <- VAR(res.oil, p=2, type="both")
serial.test(m.res.oil)
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object m.res.oil
## Chi-squared = 49.41, df = 56, p-value = 0.7208
roots(m.res.oil)
## [1] 0.9395958 0.8012263 0.3500023 0.3500023

Нет статистически значимых свидетельств, что остатки моделей коррелированы. Собственные значения сопровождающей матрицы для всех уравнений меньше 1, т.е. процессы стабильны. Для наших целей этого анализа вполне достаточно. Построим неограниченные модели VAR(p+m) для каждой пары переменных и выполним тест Вальда на равенство нулю первых p коэффициентов при “причинных” переменных (для этого вручную придется задавать параметр Terms в тесте Вальда):

VAR(1+1): курс рубля, стоимость нефти

(m.rub.oil <- VAR(rub.oil, p=1+1, type="both"))
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation rub: 
## ======================================== 
## Call:
## rub = rub.l1 + oil.l1 + rub.l2 + oil.l2 + const + trend 
## 
##        rub.l1        oil.l1        rub.l2        oil.l2         const 
##  8.902715e-01  5.411910e-05 -1.221865e-02 -4.423024e-05  4.570473e-03 
##         trend 
## -3.550585e-05 
## 
## 
## Estimated coefficients for equation oil: 
## ======================================== 
## Call:
## oil = rub.l1 + oil.l1 + rub.l2 + oil.l2 + const + trend 
## 
##       rub.l1       oil.l1       rub.l2       oil.l2        const 
## 1270.6554905    1.0934022 -367.7720688   -0.3358575  -18.8349282 
##        trend 
##    0.2176249

Проверка H0: rub ~ oil

wald.test(Sigma=vcov(m.rub.oil$varresult[[1]]), 
          b=coef(m.rub.oil$varresult[[1]]), 
          Terms=c(2))
##Wald test:
##----------
##
##Chi-squared test:
##X2 = 6.5, df = 1, P(> X2) = 0.011

Проверка H0: oil ~ rub

wald.test(Sigma=vcov(m.rub.oil$varresult[[2]]), 
          b=coef(m.rub.oil$varresult[[2]]), 
          Terms=c(1)) 
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 1.8, df = 1, P(> X2) = 0.18

VAR(2+2): курс рубля, резервы

(m.rub.res <- VAR(rub.res, p=2+2, type="both"))
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation rub: 
## ======================================== 
## Call:
## rub = rub.l1 + res.l1 + rub.l2 + res.l2 + rub.l3 + res.l3 + rub.l4 + res.l4 + const + trend 
## 
##        rub.l1        res.l1        rub.l2        res.l2        rub.l3 
##  1.102242e+00  1.871358e-09 -3.440431e-01  1.164672e-08  2.680034e-01 
##        res.l3        rub.l4        res.l4         const         trend 
## -1.931194e-08 -2.426553e-01  1.216213e-08  7.387591e-03 -6.701319e-05 
## 
## 
## Estimated coefficients for equation res: 
## ======================================== 
## Call:
## res = rub.l1 + res.l1 + rub.l2 + res.l2 + rub.l3 + res.l3 + rub.l4 + res.l4 + const + trend 
## 
##        rub.l1        res.l1        rub.l2        res.l2        rub.l3 
##  3.038377e+06  1.242807e+00 -3.368189e+06 -2.428137e-01  3.195951e+06 
##        res.l3        rub.l4        res.l4         const         trend 
## -2.034562e-01 -3.833291e+06  1.983447e-01  5.418858e+04 -3.381944e+02

Проверка H0: rub~res

wald.test(Sigma=vcov(m.rub.res$varresult[[1]]), 
          b=coef(m.rub.res$varresult[[1]]), 
          Terms=c(2, 4))
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 1.2, df = 2, P(> X2) = 0.55

Проверка H0: res ~ rub

wald.test(Sigma=vcov(m.rub.res$varresult[[2]]), 
          b=coef(m.rub.res$varresult[[2]]), 
          Terms=c(1, 3)) 
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 3.2, df = 2, P(> X2) = 0.2

VAR(2+2): резервы, стоимость нефти

(m.res.oil <- VAR(res.oil, p=2+2, type="both"))
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation res: 
## ======================================== 
## Call:
## res = res.l1 + oil.l1 + res.l2 + oil.l2 + res.l3 + oil.l3 + res.l4 + oil.l4 + const + trend 
## 
##        res.l1        oil.l1        res.l2        oil.l2        res.l3 
##  9.949240e-01  9.279481e+02 -9.462711e-02 -5.811121e+02  2.226545e-01 
##        oil.l3        res.l4        oil.l4         const         trend 
## -4.577351e+02 -1.496666e-01 -6.699770e+01  3.172796e+04 -4.402872e+01 
## 
## 
## Estimated coefficients for equation oil: 
## ======================================== 
## Call:
## oil = res.l1 + oil.l1 + res.l2 + oil.l2 + res.l3 + oil.l3 + res.l4 + oil.l4 + const + trend 
## 
##        res.l1        oil.l1        res.l2        oil.l2        res.l3 
##  3.080310e-05  1.180247e+00 -1.280450e-04 -1.755761e-01  1.364855e-04 
##        oil.l3        res.l4        oil.l4         const         trend 
## -1.817108e-01 -1.437033e-05  4.090393e-03  7.103714e+00 -5.810768e-02

Проверка H0: res ~ oil

wald.test(Sigma=vcov(m.res.oil$varresult[[1]]), 
          b=coef(m.res.oil$varresult[[1]]), 
          Terms=c(2, 4))
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 16.1, df = 2, P(> X2) = 0.00033

Проверка H0: oil ~ res

wald.test(Sigma=vcov(m.res.oil$varresult[[2]]), 
          b=coef(m.res.oil$varresult[[2]]), 
          Terms=c(1, 3)) 
## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 1.6, df = 2, P(> X2) = 0.45

Заключение

Алгоритм Тода-Ямамото довольно неочевидный; часто попадаются его некорректные трактовки ("проверить причинность на разностях временных рядов", "подбирать p наобум"). В тесте на причинность все модели строятся только на исходных данных (в том числе и неограниченные). При наличии автокорреляции в моделях надо увеличить параметр p.

Что касается интерпретации результатов тестов, то стоимость нефти является Грэнджер-причиной курса рубля (p-value = 0.011). Ни курс рубля, ни объемы международных резервов не являются "причинами" друг друга. Но в то же время цена на нефть является причиной по Грэнджеру объемов международных резервов (p-value = 0.00033).

Ссылки

Поделиться с друзьями

-->