# 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₂

# Instalacja wymaganych bibliotek
pip install pandas numpy matplotlib statsmodels

# Biblioteki do zainstalowania:
# - pandas: obsługa danych
# - numpy: obliczenia numeryczne
# - matplotlib: wizualizacja
# - statsmodels: modele VAR

import 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()

# 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()

# 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 PKB

from 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()

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)