Modele VAR — wektorowa autoregresja, IRF i dekompozycja wariancji

Zaawansowany
Streszczenie

Pełny wykład modeli VAR: forma towarzysząca i warunek stabilności (moduły wartości własnych < 1) z dowodem, reprezentacja MA(∞) i funkcje odpowiedzi na impuls, dekompozycja wariancji prognozy (FEVD), przyczynowość Grangera z dowodem, identyfikacja przez rozkład Choleskiego — z wykresami w skali szarości (2D i 3D), przykładem VAR(1) policzonym w całości, kodem R i Python, wkładem polskiej szkoły bayesowskiej i bibliografią z DOI.

Modele jednorównaniowe opisują pojedynczy szereg, lecz zmienne makroekonomiczne wzajemnie na siebie oddziałują — stopa procentowa wpływa na inflację, a inflacja na stopę. Model wektorowej autoregresji (VAR) traktuje wszystkie zmienne symetrycznie: każda jest funkcją własnych opóźnień i opóźnień pozostałych. Wykład obejmuje warunki stabilności z dowodem, reprezentację MA(∞), funkcje odpowiedzi na impuls, dekompozycję wariancji prognozy, przyczynowość Grangera oraz identyfikację strukturalną — uzupełnione przykładami liczbowymi, wykresami i kodem w językach R i Python.

Geneza i struktura

Model VAR wprowadził Christopher Sims (1980) jako odpowiedź na — jego zdaniem nierealistyczne — restrykcje wielkich modeli strukturalnych; praca Macroeconomics and Reality przyniosła mu Nagrodę Nobla w 2011 roku. Pojęcie przyczynowości predykcyjnej pochodzi od Grangera (1969).

Definicja
Model VAR(p)

Dla $k$ zmiennych i $p$ opóźnień, VAR($p$) to układ $k$ równań:

$$ \mathbf{y}_t=\mathbf{c}+\mathbf{A}_1\mathbf{y}_{t-1}+\dots+\mathbf{A}_p\mathbf{y}_{t-p}+\boldsymbol\varepsilon_t, $$

gdzie $\mathbf{y}_t$ jest wektorem $k\times1$, $\mathbf{A}_i$ — macierzami $k\times k$, a $\boldsymbol\varepsilon_t$ — wektorowym białym szumem o macierzy kowariancji $\boldsymbol\Sigma$. Dla $k=2$, $p=1$:

$$ \begin{aligned} y_{1,t}&=c_1+a_{11}y_{1,t-1}+a_{12}y_{2,t-1}+\varepsilon_{1,t},\\ y_{2,t}&=c_2+a_{21}y_{1,t-1}+a_{22}y_{2,t-1}+\varepsilon_{2,t}. \end{aligned} $$

Ponieważ każde równanie ma te same regresory (opóźnienia wszystkich zmiennych), estymacja równanie po równaniu zwykłą MNK jest zgodna i efektywna — uogólniony układ regresji pozornie niepowiązanych (SUR) pokrywa się tu z MNK.

Stabilność i reprezentacja MA(∞)

Twierdzenie
Warunek stabilności VAR(1)

Proces VAR(1) $\mathbf{y}_t=\mathbf{c}+\mathbf{A}\mathbf{y}_{t-1}+\boldsymbol\varepsilon_t$ jest stabilny (stacjonarny) wtedy i tylko wtedy, gdy wszystkie wartości własne macierzy $\mathbf{A}$ mają moduł mniejszy od $1$. Wówczas

$$ \mathbf{y}_t=\boldsymbol\mu+\sum_{j=0}^{\infty}\mathbf{A}^{j}\boldsymbol\varepsilon_{t-j},\qquad \boldsymbol\mu=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}. $$
Dowód
Stabilność, reprezentacja MA(∞) i średnia VAR(1)
  1. Podstawianie wsteczne. Iterując równanie: $\mathbf{y}_t=\mathbf{c}+\mathbf{A}(\mathbf{c}+\mathbf{A}\mathbf{y}_{t-2}+\boldsymbol\varepsilon_{t-1})+\boldsymbol\varepsilon_t$. Po $m$ krokach $\mathbf{y}_t=\sum_{j=0}^{m-1}\mathbf{A}^j\mathbf{c}+\mathbf{A}^m\mathbf{y}_{t-m}+\sum_{j=0}^{m-1}\mathbf{A}^j\boldsymbol\varepsilon_{t-j}$.
  2. Rola wartości własnych. Rozkładając $\mathbf{A}=\mathbf{Q}\boldsymbol\Lambda\mathbf{Q}^{-1}$, mamy $\mathbf{A}^m=\mathbf{Q}\boldsymbol\Lambda^m\mathbf{Q}^{-1}$, gdzie $\boldsymbol\Lambda^m=\operatorname{diag}(\lambda_i^m)$. Gdy każdy $|\lambda_i|<1$, to $\lambda_i^m\to0$, więc $\mathbf{A}^m\to\mathbf{0}$.
  3. Zanik członu początkowego. Zatem $\mathbf{A}^m\mathbf{y}_{t-m}\to\mathbf{0}$ (w sensie średniokwadratowym), a szereg $\sum_{j=0}^{\infty}\mathbf{A}^j\boldsymbol\varepsilon_{t-j}$ jest zbieżny, bo normy $\lVert\mathbf{A}^j\rVert$ maleją geometrycznie. Daje to reprezentację MA(∞).
  4. Suma deterministyczna. $\sum_{j=0}^{\infty}\mathbf{A}^j\mathbf{c}=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}$ — wektorowy szereg geometryczny zbieżny dokładnie wtedy, gdy $|\lambda_i|<1$ (wówczas $\mathbf{I}-\mathbf{A}$ jest odwracalna, bo nie ma wartości własnej $1$).
  5. Średnia. Z reprezentacji MA(∞): $\mathbb{E}[\mathbf{y}_t]=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}+\sum_j\mathbf{A}^j\mathbb{E}[\boldsymbol\varepsilon_{t-j}]=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}=\boldsymbol\mu$. Gdy któryś $|\lambda_i|\ge1$, szeregi rozbiegają się i proces jest niestacjonarny.
∎
Trójkąt stabilności VAR(1)
Trójkąt stabilności VAR(1) w przestrzeni śladu i wyznacznika macierzy $\mathbf{A}$. Obie wartości własne mają moduł poniżej $1$ tylko wewnątrz trójkąta wyznaczonego przez $\det<1$, $\det>\mathrm{tr}-1$ i $\det>-\mathrm{tr}-1$. Punkt przykładowego modelu $(\mathrm{tr},\det)=(0{,}9;\,0{,}18)$ leży wewnątrz — proces jest stabilny.
Przykład
Stabilność i średnia

Niech $\mathbf{A}=\begin{bmatrix}0{,}5&0{,}2\\0{,}1&0{,}4\end{bmatrix}$, $\mathbf{c}=(1,2)$. Wartości własne wyznacza równanie $\lambda^2-\mathrm{tr}\,\lambda+\det=0$, czyli $\lambda^2-0{,}9\lambda+0{,}18=0$:

$$ \lambda=\frac{0{,}9\pm\sqrt{0{,}81-0{,}72}}{2}=\frac{0{,}9\pm0{,}3}{2}=\{0{,}6;\ 0{,}3\}. $$

Oba moduły są poniżej $1$, więc proces jest stabilny. Średnia długookresowa:

$$ \boldsymbol\mu=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}=\frac{1}{0{,}28}\begin{bmatrix}0{,}6&0{,}2\\0{,}1&0{,}5\end{bmatrix}\begin{bmatrix}1\\2\end{bmatrix}=\frac{1}{0{,}28}\begin{bmatrix}1{,}0\\1{,}1\end{bmatrix}=\begin{bmatrix}3{,}571\\3{,}929\end{bmatrix}. $$

Funkcje odpowiedzi na impuls

Pojedyncze współczynniki VAR są trudne do interpretacji. Z reprezentacji MA(∞) wynika natomiast jasne narzędzie: macierze $\mathbf{A}^j$ to funkcje odpowiedzi na impuls (IRF), opisujące reakcję układu na szok.

Twierdzenie
IRF i identyfikacja Choleskiego
W reprezentacji $\mathbf{y}_t=\boldsymbol\mu+\sum_j\mathbf{A}^j\boldsymbol\varepsilon_{t-j}$ element $(i,l)$ macierzy $\mathbf{A}^h=\partial\mathbf{y}_{t+h}/\partial\boldsymbol\varepsilon_t$ to odpowiedź zmiennej $i$ po $h$ okresach na jednostkowy szok zmiennej $l$ (postać zredukowana). Ponieważ szoki zredukowane są skorelowane ($\boldsymbol\Sigma$ niediagonalna), w praktyce stosuje się szoki ortogonalne $\mathbf{u}_t=\mathbf{P}^{-1}\boldsymbol\varepsilon_t$, gdzie $\mathbf{P}$ pochodzi z rozkładu Choleskiego $\boldsymbol\Sigma=\mathbf{P}\mathbf{P}^\top$; strukturalna IRF to wtedy $\mathbf{A}^h\mathbf{P}$.

Kolejność zmiennych ma znaczenie.

Rozkład Choleskiego zależy od kolejności zmiennych: zakłada, że zmienna wcześniejsza oddziałuje na późniejszą natychmiast, lecz nie odwrotnie. Zmiana kolejności zmienia strukturalne IRF i FEVD. Zmienne porządkuje się od najbardziej egzogenicznej do najbardziej endogenicznej (polityka pieniężna zwykle na końcu) i zaleca się weryfikację wrażliwości wyników na permutacje.

Przykład
Przykład — funkcje odpowiedzi na impuls

Dla macierzy z przykładu 1 liczymy kolejne potęgi. $\mathbf{A}^0=\mathbf{I}$. $\mathbf{A}^1=\mathbf{A}$. Następnie

$$ \mathbf{A}^2=\mathbf{A}\mathbf{A}=\begin{bmatrix}0{,}5\cdot0{,}5+0{,}2\cdot0{,}1&0{,}5\cdot0{,}2+0{,}2\cdot0{,}4\\0{,}1\cdot0{,}5+0{,}4\cdot0{,}1&0{,}1\cdot0{,}2+0{,}4\cdot0{,}4\end{bmatrix}=\begin{bmatrix}0{,}27&0{,}18\\0{,}09&0{,}18\end{bmatrix}. $$

Odpowiedź $y_1$ na szok własny $\varepsilon_1$ (element $(1,1)$) w kolejnych okresach: $1{,}00$, $0{,}50$, $0{,}27$, $0{,}153$, $0{,}089$, … — monotoniczny zanik. Odpowiedź $y_1$ na szok $\varepsilon_2$ (element $(1,2)$): $0$, $0{,}20$, $0{,}18$, $0{,}126$, … — najpierw rośnie (szok dociera z opóźnieniem), potem wygasa. Wszystkie reakcje zbiegają do zera, bo proces jest stabilny.

Funkcje odpowiedzi na impuls VAR(1)
Cztery funkcje odpowiedzi na impuls modelu VAR(1). Reakcje na szoki własne ($y_1\!\leftarrow\!\varepsilon_1$, $y_2\!\leftarrow\!\varepsilon_2$) startują od jedności i maleją; reakcje krzyżowe ($y_1\!\leftarrow\!\varepsilon_2$, $y_2\!\leftarrow\!\varepsilon_1$) startują od zera, rosną z opóźnieniem i wygasają. Geometryczny zanik wszystkich krzywych jest konsekwencją stabilności (moduły wartości własnych $<1$).
Powierzchnia odpowiedzi na impuls a trwałość
Odpowiedź $y_1$ na szok własny w funkcji horyzontu $h$ oraz trwałości $a_{11}$. Przy małym $a_{11}$ reakcja gaśnie natychmiast; im bliżej $a_{11}=1$, tym wolniejszy zanik i dłuższa pamięć układu. Powierzchnia obrazuje, jak jeden parametr kształtuje cały profil dynamicznej odpowiedzi.

Dekompozycja wariancji prognozy

Twierdzenie
FEVD

Błąd prognozy $h$-okresowej to $\mathbf{y}_{t+h}-\hat{\mathbf{y}}_{t+h}=\sum_{j=0}^{h-1}\mathbf{A}^j\boldsymbol\varepsilon_{t+h-j}$. Dla szoków ortogonalnych o wariancji jednostkowej wariancja błędu prognozy zmiennej $i$ rozkłada się na wkłady poszczególnych szoków:

$$ \text{udział szoku } l=\frac{\sum_{j=0}^{h-1}\big([\mathbf{A}^j]_{il}\big)^2}{\sum_{l'}\sum_{j=0}^{h-1}\big([\mathbf{A}^j]_{il'}\big)^2}. $$

FEVD odpowiada na pytanie, jaka część niepewności prognozy danej zmiennej pochodzi od szoków każdej ze zmiennych.

Przykład
FEVD zmiennej y₁

Przyjmując ortogonalne szoki jednostkowe ($\boldsymbol\Sigma=\mathbf{I}$), z potęg $\mathbf{A}^j$ z przykładu 2 obliczamy wkład szoku własnego $\varepsilon_1$ i szoku $\varepsilon_2$ do wariancji prognozy $y_1$:

Horyzont $h$szok $\varepsilon_1$szok $\varepsilon_2$
1$100{,}0\%$$0{,}0\%$
2$96{,}9\%$$3{,}1\%$
4$93{,}8\%$$6{,}2\%$
12$93{,}2\%$$6{,}8\%$

W horyzoncie jednookresowym cała niepewność $y_1$ pochodzi z szoku własnego (szok $y_2$ jeszcze nie zdążył oddziałać). Z czasem udział szoku $y_2$ rośnie, stabilizując się na poziomie około $7\%$.

Dekompozycja wariancji prognozy FEVD
Dekompozycja wariancji prognozy $y_1$: ciemny obszar to udział szoku własnego $\varepsilon_1$, jasny — szoku $\varepsilon_2$. W krótkim horyzoncie dominuje szok własny; udział szoku $y_2$ narasta z horyzontem i stabilizuje się. Suma udziałów zawsze wynosi $100\%$.

Przyczynowość Grangera

Twierdzenie
Przyczynowość Grangera w VAR(1)
Zmienna $y_2$ nie jest przyczyną Grangera dla $y_1$ wtedy i tylko wtedy, gdy współczynnik $a_{12}=0$ — wówczas przeszłość $y_2$ nie poprawia prognozy $y_1$ ponad to, co daje przeszłość $y_1$.
Dowód
Brak przyczynowości Grangera ⟺ a₁₂ = 0
  1. Optymalna prognoza. Prognoza jednookresowa o minimalnym błędzie to warunkowa wartość oczekiwana $\mathbb{E}[y_{1,t}\mid\mathcal{F}_{t-1}]=c_1+a_{11}y_{1,t-1}+a_{12}y_{2,t-1}$.
  2. Rola $y_{2,t-1}$. Człon $a_{12}y_{2,t-1}$ jest jedynym miejscem, w którym przeszłość $y_2$ wchodzi do prognozy $y_1$. Jeżeli $a_{12}=0$, prognoza redukuje się do $c_1+a_{11}y_{1,t-1}$ — zależy wyłącznie od przeszłości $y_1$.
  3. Równoważność. Zatem dołączenie przeszłości $y_2$ zmniejsza błąd prognozy $y_1$ dokładnie wtedy, gdy $a_{12}\neq0$. Definicja Grangera (1969) — $y_2$ jest przyczyną $y_1$, jeśli poprawia jego prognozę — przekłada się więc na warunek $a_{12}\neq0$.
∎

Przyczynowość Grangera to nie przyczynowość.

Test Grangera mierzy predykcyjność, nie sprawczość strukturalną. Zmienna $X$ może być przyczyną Grangera dla $Y$ bez bycia jej rzeczywistą przyczyną — na przykład przez wspólną trzecią zmienną albo przez wyprzedzające oczekiwania (ceny aktywów „granger-powodują" przyszłe zyski, bo dyskontują oczekiwania). Wniosek przyczynowy wymaga argumentu ekonomicznego, nie samego testu.

Diagram przyczynowości Grangera w VAR
Struktura zależności przykładowego VAR(1). Strzałki krzyżowe odpowiadają współczynnikom $a_{12}$ i $a_{21}$, pętle — autoregresji $a_{11}$, $a_{22}$. Ponieważ $a_{12}=0{,}2\neq0$, przeszłość $y_2$ poprawia prognozę $y_1$, więc $y_2$ jest przyczyną Grangera dla $y_1$.
Przykład
Test przyczynowości
W przykładowym modelu $a_{12}=0{,}2\neq0$, więc $y_2$ jest przyczyną Grangera dla $y_1$ — co potwierdza dodatnia reakcja krzyżowa IRF $y_1\!\leftarrow\!\varepsilon_2$ z przykładu 2. Z kolei $a_{21}=0{,}1\neq0$, więc również $y_1$ jest przyczyną Grangera dla $y_2$ — mamy sprzężenie zwrotne. Gdyby któryś z tych współczynników był zerowy, kierunek przyczynowości byłby jednostronny.

Procedura modelowania i rozszerzenia

Praktyka VAR przebiega w pięciu krokach: (1) weryfikacja stacjonarności każdego szeregu testem ADF i różnicowanie lub model VECM przy kointegracji; (2) dobór rzędu $p$ kryteriami informacyjnymi (AIC skłania do większego $p$, BIC do mniejszego); (3) estymacja i diagnostyka reszt (autokorelacja, normalność, ARCH); (4) interpretacja przez IRF, FEVD i przyczynowość Grangera, nie surowe współczynniki; (5) prognozowanie z przedziałami.

RozszerzenieZastosowanie
VECMzmienne skointegrowane (relacja długookresowa)
SVARstrukturalny VAR — restrykcje ekonomiczne na szoki
VARXVAR ze zmiennymi egzogenicznymi
BVARwiele zmiennych, mało danych (priory Minnesota)
TVP-VARwspółczynniki zmienne w czasie

Wkład polskiej szkoły

Modele wektorowe są podstawowym narzędziem prognostycznym Narodowego Banku Polskiego, obok modeli klasy DSGE. Szczególny wkład ma polska szkoła ekonometrii bayesowskiej, rozwijana na Uniwersytecie Ekonomicznym w Krakowie przez Jacka Osiewalskiego i współpracowników; bayesowskie modele VAR (BVAR) z priorami typu Minnesota oraz bayesowskie modele zmienności (MSV, SV) są wykorzystywane do prognozowania makroekonomicznego i analizy rynków finansowych w warunkach krótkich szeregów i wielu zmiennych. Podejście bayesowskie łagodzi problem nadmiernej liczby parametrów, który ogranicza klasyczne VAR.

Podsumowanie

  • Struktura: $\mathbf{y}_t=\mathbf{c}+\sum_{i}\mathbf{A}_i\mathbf{y}_{t-i}+\boldsymbol\varepsilon_t$; estymacja równanie po równaniu MNK.
  • Stabilność: wszystkie $|\lambda_i(\mathbf{A})|<1$; wtedy reprezentacja MA(∞) i średnia $(\mathbf{I}-\mathbf{A})^{-1}\mathbf{c}$.
  • IRF: $\mathbf{A}^h$ (zredukowane) lub $\mathbf{A}^h\mathbf{P}$ (strukturalne, Cholesky — zależne od kolejności).
  • FEVD: podział wariancji prognozy na wkłady szoków.
  • Granger: $y_2$ nie jest przyczyną $y_1$ ⟺ $a_{12}=0$; predykcyjność, nie sprawczość.
  • Rozszerzenia: VECM, SVAR, BVAR; polska szkoła bayesowska (UEK Kraków, NBP).

Bibliografia

  • Granger, C. W. J. (1969). Investigating Causal Relations by Econometric Models and Cross-spectral Methods. Econometrica, 37(3), 424–438. https://doi.org/10.2307/1912791
  • Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.
  • Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer. https://doi.org/10.1007/978-3-540-27752-1
  • Sims, C. A. (1980). Macroeconomics and Reality. Econometrica, 48(1), 1–48. https://doi.org/10.2307/1912017
  • Osiewalski, J. (2001). Ekonometria bayesowska w zastosowaniach. Wydawnictwo Akademii Ekonomicznej w Krakowie.

Dalej: Szeregi czasowe · Autokorelacja · Dane panelowe

Literatura uzupełniająca
  • Podręcznik: H. Lütkepohl, New Introduction to Multiple Time Series Analysis, Springer, 2005
  • Podręcznik: J. D. Hamilton, Time Series Analysis, Princeton University Press, 1994, rozdz. 10–11
  • C. A. Sims, „Macroeconomics and Reality", Econometrica 48(1), 1980
  • J. Osiewalski, prace ze szkoły ekonometrii bayesowskiej (Uniwersytet Ekonomiczny w Krakowie)
  • Pakiety: R vars, tsDyn; Python statsmodels.tsa.api.VAR
Oprogramowanie

R:

library(vars)
## przyklad VAR(1): symulacja i estymacja
A <- matrix(c(0.5,0.1, 0.2,0.4), 2, 2); c0 <- c(1,2)
set.seed(1); T <- 300; Y <- matrix(0, T, 2)
for (t in 2:T) Y[t,] <- c0 + A %*% Y[t-1,] + rnorm(2)
m <- VAR(Y, p = 1, type = "const")
roots(m)                       # moduly < 1 -> stabilny
irf(m, n.ahead = 12)           # funkcje odpowiedzi na impuls
fevd(m, n.ahead = 12)          # dekompozycja wariancji prognozy
causality(m, cause = "y2")     # przyczynowosc Grangera

Python:

import numpy as np
from statsmodels.tsa.api import VAR
A = np.array([[0.5,0.2],[0.1,0.4]]); c0 = np.array([1,2])
rng = np.random.default_rng(1); T = 300; Y = np.zeros((T,2))
for t in range(1,T): Y[t] = c0 + A @ Y[t-1] + rng.standard_normal(2)
res = VAR(Y).fit(1)
print(np.abs(np.linalg.eigvals(A)))   # < 1 -> stabilny
print(res.irf(12).irfs)               # IRF
print(res.fevd(12).decomp)            # FEVD
print(res.test_causality('y1', 'y2', kind='f'))