Vector Autoregression (VAR) to wielowymiarowe rozszerzenie modelu AR, które pozwala modelować kilka powiązanych szeregów czasowych jednocześnie.
Kluczowe cechy:
ARIMA: Model jednowymiarowy - jedna zmienna w czasie
VAR: Model wielowymiarowy - wiele powiązanych zmiennych w czasie
Przykłady zastosowań VAR:
# Model VAR(p) dla n zmiennych:
# Y_t = c + A₁·Y_{t-1} + A₂·Y_{t-2} + ... + A_p·Y_{t-p} + ε_t
# Gdzie:
# Y_t - wektor n zmiennych w czasie t: [y₁_t, y₂_t, ..., yₙ_t]ᵀ
# c - wektor stałych: [c₁, c₂, ..., cₙ]ᵀ
# Aᵢ - macierz współczynników (n × n) dla opóźnienia i
# ε_t - wektor błędów losowych: [ε₁_t, ε₂_t, ..., εₙ_t]ᵀ
# Przykład VAR(1) dla 2 zmiennych (n=2, p=1):
# [y₁_t] [c₁] [a₁₁ a₁₂] [y₁_{t-1}] [ε₁_t]
# [y₂_t] = [c₂] + [a₂₁ a₂₂]·[y₂_{t-1}] + [ε₂_t]
# W rozwinięciu:
# y₁_t = c₁ + a₁₁·y₁_{t-1} + a₁₂·y₂_{t-1} + ε₁_t
# y₂_t = c₂ + a₂₁·y₁_{t-1} + a₂₂·y₂_{t-1} + ε₂_t
# Zauważ: y₁ zależy od przeszłości y₁ i y₂
# y₂ zależy od przeszłości y₁ i y₂Aby model VAR był poprawny, muszą być spełnione następujące założenia:
Jeśli szeregi nie są stacjonarne, używamy VECM (Vector Error Correction Model) lub różnicujemy dane przed estymacją.
# Instalacja wymaganych bibliotek
pip install pandas numpy matplotlib statsmodels
# Biblioteki do zainstalowania:
# - pandas: obsługa danych
# - numpy: obliczenia numeryczne
# - matplotlib: wizualizacja
# - statsmodels: modele VARimport pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller
# Wczytanie danych (przykład: PKB i inflacja)
# Załóżmy, że mamy dane miesięczne
dane = pd.read_csv('dane_makro.csv', parse_dates=['data'], index_col='data')
# Wybór zmiennych do modelu VAR
zmienne = ['PKB', 'inflacja', 'stopa_procentowa']
df_var = dane[zmienne]
# Sprawdzenie danych
print(df_var.head())
print(f"\nOkres: {df_var.index.min()} do {df_var.index.max()}")
print(f"Liczba obserwacji: {len(df_var)}")
# Wizualizacja
df_var.plot(subplots=True, figsize=(12, 8))
plt.tight_layout()
plt.show()def test_stacjonarnosci_var(df):
"""
Testuje stacjonarność wszystkich zmiennych w DataFrame.
"""
print("=== TESTY STACJONARNOŚCI (ADF) ===\n")
wyniki = {}
for kolumna in df.columns:
wynik = adfuller(df[kolumna].dropna(), autolag='AIC')
p_value = wynik[1]
wyniki[kolumna] = p_value
print(f"{kolumna}:")
print(f" Statystyka ADF: {wynik[0]:.4f}")
print(f" p-value: {p_value:.4f}")
if p_value <= 0.05:
print(f" ✓ STACJONARNY")
else:
print(f" ✗ NIESTACJONARNY - wymagane różnicowanie")
print()
return wyniki
# Użycie
wyniki = test_stacjonarnosci_var(df_var)
# Jeśli zmienne niestacjonarne, różnicowanie:
df_var_diff = df_var.diff().dropna()# Tworzenie modelu VAR
model = VAR(df_var_diff)
# Automatyczny dobór optymalnego rzędu opóźnienia
# Testuje p od 0 do maxlags
wybor_p = model.select_order(maxlags=15)
print(wybor_p.summary())
# Wyświetl kryteria informacyjne dla różnych p
print("\nKryteria informacyjne:")
print(wybor_p)
# Najlepsze p według różnych kryteriów:
print(f"\nNajlepsze p według AIC: {wybor_p.aic}")
print(f"Najlepsze p według BIC: {wybor_p.bic}")
print(f"Najlepsze p według HQIC: {wybor_p.hqic}")
# Zazwyczaj wybieramy według BIC (bardziej konserwatywne)
optymalny_p = wybor_p.bic
print(f"\n✓ Wybrany rząd opóźnienia: p = {optymalny_p}")# Estymacja modelu VAR z wybranym p
model_var = model.fit(maxlags=optymalny_p)
# Podsumowanie modelu
print(model_var.summary())
# Współczynniki modelu
print("\n=== WSPÓŁCZYNNIKI ===")
print(model_var.params)
# Macierze współczynników A₁, A₂, ..., A_p
print("\n=== MACIERZE WSPÓŁCZYNNIKÓW ===")
for i in range(1, optymalny_p + 1):
print(f"\nMacierz A_{i}:")
print(model_var.coefs[i-1])
# R² dla każdej równania
print("\n=== R² DLA KAŻDEJ ZMIENNEJ ===")
for nazwa, r2 in zip(df_var_diff.columns, model_var.rsquared):
print(f"{nazwa}: R² = {r2:.4f}")# 1. Test autokorelacji reszt (Portmanteau)
# H0: Brak autokorelacji
test_auto = model_var.test_whiteness(nlags=10)
print("Test Portmanteau (autokorelacja):")
print(f"p-value: {test_auto.pvalue:.4f}")
if test_auto.pvalue > 0.05:
print("✓ Brak autokorelacji reszt\n")
else:
print("✗ Występuje autokorelacja\n")
# 2. Test normalności reszt (Jarque-Bera)
# H0: Reszty mają rozkład normalny
test_norm = model_var.test_normality()
print("Test normalności (Jarque-Bera):")
print(f"p-value: {test_norm.pvalue:.4f}")
if test_norm.pvalue > 0.05:
print("✓ Reszty są normalne\n")
else:
print("✗ Reszty nie są normalne\n")
# 3. Analiza reszt
reszty = model_var.resid
print(f"Średnia reszt: {reszty.mean().values}")
print(f"Odchylenie std reszt: {reszty.std().values}")# Wykresy reszt dla każdej zmiennej
fig, axes = plt.subplots(len(df_var_diff.columns), 1,
figsize=(12, 8))
for i, kolumna in enumerate(df_var_diff.columns):
axes[i].plot(reszty.index, reszty[kolumna])
axes[i].axhline(y=0, color='r', linestyle='--')
axes[i].set_title(f'Reszty dla {kolumna}')
axes[i].set_ylabel('Reszty')
axes[i].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Macierz korelacji reszt
print("\nMacierz korelacji reszt:")
print(reszty.corr())
# ACF reszt
from statsmodels.graphics.tsaplots import plot_acf
fig, axes = plt.subplots(len(df_var_diff.columns), 1,
figsize=(12, 8))
for i, kolumna in enumerate(df_var_diff.columns):
plot_acf(reszty[kolumna], ax=axes[i], lags=20)
axes[i].set_title(f'ACF reszt - {kolumna}')
plt.tight_layout()
plt.show()# Prognoza na n_steps okresów do przodu
n_steps = 12 # np. 12 miesięcy
# Model wymaga ostatnich p obserwacji jako dane wejściowe
lag_order = model_var.k_ar
dane_wejsciowe = df_var_diff.values[-lag_order:]
# Generowanie prognozy
prognoza = model_var.forecast(dane_wejsciowe, steps=n_steps)
# Konwersja do DataFrame
idx = pd.date_range(start=df_var_diff.index[-1],
periods=n_steps+1,
freq='M')[1:]
df_prognoza = pd.DataFrame(prognoza,
index=idx,
columns=df_var_diff.columns)
print("Prognoza:")
print(df_prognoza)
# UWAGA: Jeśli różnicowaliśmy dane, trzeba cofnąć różnicowanie!
# df_prognoza_oryginalna = df_prognoza.cumsum() + ostatnia_wartość# Wizualizacja danych historycznych i prognozy
fig, axes = plt.subplots(len(df_var_diff.columns), 1,
figsize=(12, 10))
for i, kolumna in enumerate(df_var_diff.columns):
# Dane historyczne
axes[i].plot(df_var_diff.index, df_var_diff[kolumna],
label='Dane historyczne', color='blue')
# Prognoza
axes[i].plot(df_prognoza.index, df_prognoza[kolumna],
label='Prognoza', color='red', linestyle='--')
axes[i].axvline(x=df_var_diff.index[-1],
color='green', linestyle=':',
label='Początek prognozy')
axes[i].set_title(f'{kolumna}')
axes[i].set_ylabel('Wartość')
axes[i].legend()
axes[i].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Impulse Response Function (IRF) pokazuje, jak szok (nagła zmiana) w jednej zmiennej wpływa na wszystkie zmienne w systemie w czasie.
Pytania na które odpowiada IRF:
# Obliczenie IRF dla n_periods okresów
irf = model_var.irf(periods=20)
# Wykres IRF - jak szok w każdej zmiennej wpływa na inne
irf.plot(orth=False)
plt.tight_layout()
plt.show()
# Ortogonalizowane IRF (Cholesky decomposition)
# Usuwa współczesną korelację między zmiennymi
irf.plot(orth=True)
plt.tight_layout()
plt.show()
# IRF dla konkretnej pary zmiennych
# Przykład: wpływ szoku w inflacji na PKB
irf_pkb_inflacja = irf.irfs[:, 0, 1] # [czas, zmienna_target, zmienna_impulse]
plt.figure(figsize=(10, 5))
plt.plot(irf_pkb_inflacja)
plt.axhline(y=0, color='r', linestyle='--')
plt.title('Wpływ szoku w inflacji na PKB')
plt.xlabel('Okresy po szoku')
plt.ylabel('Odpowiedź PKB')
plt.grid(True)
plt.show()FEVD pokazuje, jaka część wariancji błędu prognozy dla każdej zmiennej jest wyjaśniona przez szoki w innych zmiennych.
Odpowiada na pytanie: Które zmienne są najważniejsze dla wyjaśnienia zmienności danej zmiennej?
# Obliczenie FEVD
fevd = model_var.fevd(periods=20)
# Wykres FEVD
fevd.plot()
plt.tight_layout()
plt.show()
# Tabela FEVD dla wybranego okresu
print("\nFEVD dla okresu 10:")
print(fevd.summary().head(10))
# Interpretacja: jeśli FEVD pokazuje, że 60% wariancji
# PKB po 10 okresach wyjaśnia inflacja, to inflacja ma
# silny wpływ na PKBfrom statsmodels.tsa.stattools import grangercausalitytests
# Test przyczynowości Grangera
# H0: X nie powoduje Y w sensie Grangera
# Jeśli p < 0.05, odrzucamy H0 - X pomaga przewidywać Y
maxlag = 4
print("=== TEST PRZYCZYNOWOŚCI GRANGERA ===\n")
zmienne = df_var_diff.columns
for i, zmienna1 in enumerate(zmienne):
for zmienna2 in zmienne:
if zmienna1 != zmienna2:
print(f"\nCzy {zmienna1} powoduje {zmienna2}?")
dane_test = df_var_diff[[zmienna2, zmienna1]]
wynik = grangercausalitytests(dane_test, maxlag=maxlag,
verbose=False)
# Sprawdź p-value dla każdego opóźnienia
for lag in range(1, maxlag+1):
p_value = wynik[lag][0]['ssr_ftest'][1]
if p_value < 0.05:
print(f" Lag {lag}: p={p_value:.4f} ✓ TAK")
else:
print(f" Lag {lag}: p={p_value:.4f} ✗ NIE")import pandas as pd
import numpy as np
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller
# Generowanie przykładowych danych (w praktyce: wczytaj z pliku)
np.random.seed(42)
daty = pd.date_range(start='2010-01-01', end='2023-12-31', freq='Q')
# Symulacja danych makroekonomicznych
pkb = 100 + np.cumsum(np.random.normal(0.5, 1, len(daty)))
inflacja = 2 + np.cumsum(np.random.normal(0.1, 0.3, len(daty)))
bezrobocie = 5 + np.cumsum(np.random.normal(-0.05, 0.2, len(daty)))
dane = pd.DataFrame({
'PKB': pkb,
'inflacja': inflacja,
'bezrobocie': bezrobocie
}, index=daty)
print("Dane:")
print(dane.head(10))
# Wizualizacja
dane.plot(subplots=True, figsize=(12, 8))
plt.tight_layout()
plt.show()# Testowanie stacjonarności
print("=== TESTY STACJONARNOŚCI ===\n")
for kolumna in dane.columns:
wynik = adfuller(dane[kolumna])
print(f"{kolumna}: p-value = {wynik[1]:.4f}")
if wynik[1] > 0.05:
print(" ✗ Niestacjonarny\n")
else:
print(" ✓ Stacjonarny\n")
# Różnicowanie
dane_diff = dane.diff().dropna()
# Ponowny test
print("\n=== PO RÓŻNICOWANIU ===\n")
for kolumna in dane_diff.columns:
wynik = adfuller(dane_diff[kolumna])
print(f"{kolumna}: p-value = {wynik[1]:.4f}")
if wynik[1] <= 0.05:
print(" ✓ Stacjonarny\n")
# Wizualizacja po różnicowaniu
dane_diff.plot(subplots=True, figsize=(12, 8))
plt.suptitle('Dane po różnicowaniu')
plt.tight_layout()
plt.show()# Dobór rzędu opóźnienia
model = VAR(dane_diff)
wybor = model.select_order(maxlags=8)
print(wybor.summary())
optimal_p = wybor.bic
print(f"\nWybrany rząd: p = {optimal_p}")
# Estymacja modelu
model_fitted = model.fit(maxlags=optimal_p)
print("\n" + "="*60)
print(model_fitted.summary())
# Diagnostyka
print("\n=== DIAGNOSTYKA ===")
test_auto = model_fitted.test_whiteness(nlags=10)
print(f"Test autokorelacji p-value: {test_auto.pvalue:.4f}")
test_norm = model_fitted.test_normality()
print(f"Test normalności p-value: {test_norm.pvalue:.4f}")
# R² dla każdego równania
print("\nR² dla każdej zmiennej:")
for nazwa, r2 in zip(dane_diff.columns, model_fitted.rsquared):
print(f" {nazwa}: {r2:.4f}")# Prognoza na 8 kwartałów
n_forecast = 8
lag_order = model_fitted.k_ar
forecast_input = dane_diff.values[-lag_order:]
forecast = model_fitted.forecast(forecast_input, steps=n_forecast)
# DataFrame z prognozą
forecast_idx = pd.date_range(start=dane_diff.index[-1],
periods=n_forecast+1, freq='Q')[1:]
df_forecast = pd.DataFrame(forecast,
index=forecast_idx,
columns=dane_diff.columns)
print("Prognoza (różnice):")
print(df_forecast)
# IRF - Impulse Response Function
irf = model_fitted.irf(15)
irf.plot(orth=True)
plt.suptitle('Funkcje odpowiedzi na impuls (IRF)')
plt.tight_layout()
plt.show()
# FEVD
fevd = model_fitted.fevd(15)
fevd.plot()
plt.suptitle('Dekompozycja wariancji (FEVD)')
plt.tight_layout()
plt.show()VECM (Vector Error Correction Model) to rozszerzenie VAR dla niestacjonarnych zmiennych, które są skointegrowane.
Kointegracja oznacza, że mimo że zmienne są niestacjonarne, istnieje między nimi długoterminowa równowaga.
Przykłady kointegracji:
from statsmodels.tsa.vector_ar.vecm import coint_johansen
# Test Johansena dla kointegracji
# H0: brak kointegracji
def test_johansena(df, det_order=0, k_ar_diff=1):
"""
det_order: -1 (brak trendu), 0 (const), 1 (trend)
k_ar_diff: opóźnienia dla różnic
"""
wynik = coint_johansen(df, det_order=det_order, k_ar_diff=k_ar_diff)
print("=== TEST KOINTEGRACJI JOHANSENA ===\n")
print("Liczba zmiennych:", df.shape[1])
print("\nStatystyka śladu (Trace statistic):")
for i in range(len(wynik.lr1)):
print(f"\nr = {i}:")
print(f" Statystyka: {wynik.lr1[i]:.2f}")
print(f" Wartość krytyczna (5%): {wynik.cvt[i, 1]:.2f}")
if wynik.lr1[i] > wynik.cvt[i, 1]:
print(f" ✓ Odrzucamy H0 - jest kointegracja")
else:
print(f" ✗ Nie odrzucamy H0")
return wynik
# Użycie na niestacjonarnych danych
wynik_johansen = test_johansena(dane, det_order=0, k_ar_diff=2)from statsmodels.tsa.vector_ar.vecm import VECM
# Jeśli test Johansena pokazał kointegrację:
# Estymacja VECM z k_ar_diff opóźnieniami i coint_rank relacjami
coint_rank = 1 # liczba relacji kointegracyjnych (z testu Johansena)
k_ar_diff = 2 # opóźnienia dla różnic
model_vecm = VECM(dane,
k_ar_diff=k_ar_diff,
coint_rank=coint_rank,
deterministic='ci') # stała w relacji kointegracyjnej
vecm_fitted = model_vecm.fit()
# Podsumowanie
print(vecm_fitted.summary())
# Wektor kointegracyjny (długoterminowa równowaga)
print("\nWektor kointegracyjny (beta):")
print(vecm_fitted.beta)
# Macierz dostosowania (jak szybko system wraca do równowagi)
print("\nMacierz dostosowania (alpha):")
print(vecm_fitted.alpha)Przed estymacją: Testuj stacjonarność wszystkich zmiennych, wizualizuj dane, sprawdź korelacje.
Podczas modelowania: Używaj kryteriów BIC do wyboru p, sprawdzaj diagnostykę reszt, testuj stabilność modelu.
Po estymacji: Analizuj IRF i FEVD, przeprowadź test Grangera, waliduj prognozy out-of-sample.
Typowe błędy: Modelowanie niestacjonarnych szeregów bez różnicowania, ignorowanie kointegracji, zbyt duży rząd p (overfitting).
ARIMA: Jedna zmienna, nie wymaga innych szeregów, dobry dla izolowanych prognoz.
VAR: Wiele zmiennych stacjonarnych, modeluje wzajemne oddziaływania, pozwala na analizę IRF.
VECM: Wiele zmiennych niestacjonarnych ze związkami długoterminowymi, modeluje kointegrację i korekty błędów.
Wybór modelu: Jeśli jedna zmienna → ARIMA. Wiele zmiennych stacjonarnych → VAR. Wiele niestacjonarnych ze związkami → VECM.
VAR(p) - wielowymiarowy model szeregów czasowych modelujący wzajemne zależności między zmiennymi.
Wymagania: Stacjonarność, dobór p (BIC/AIC), diagnostyka reszt.
Narzędzia analizy: IRF (wpływ szoku), FEVD (dekompozycja wariancji), test Grangera (przyczynowość).
VECM: Rozszerzenie dla niestacjonarnych zmiennych z kointegracją.
Biblioteki: statsmodels.tsa.api (VAR, VECM), statsmodels.tsa.stattools
Książki: "Time Series Analysis" (Hamilton), "Applied Econometric Time Series" (Enders), "New Introduction to Multiple Time Series Analysis" (Lütkepohl)
Dokumentacja: statsmodels.org/stable/vector_ar.html, statsmodels VAR Tutorial
Kursy: Econometrics Academy (YouTube), Coursera Time Series Specialization
Dane: FRED Economic Data, Yahoo Finance, Kaggle Economics datasets