Prognozowanie PKB — modele ARIMA w praktyce

Streszczenie

Praktyczny tutorial prognozowania szeregu czasowego: od wczytania danych PKB przez test stacjonarności i dobór ARIMA do prognozy z przedziałami ufności — w R i Python

Cel tutorialu

Celem jest zbudowanie modelu prognozującego PKB na kilka kwartałów do przodu. Analizę przeprowadzimy w pełnym zakresie: od wczytania danych, przez diagnozę i estymację modelu, po walidację i prognozę.

To kontynuacja przewodnika o regresji, teraz dla szeregów czasowych.

Krok 1: Wczytanie danych

# R — pobierz realne PKB z FRED (lub wczytaj CSV)
library(forecast)

# Przykład: wbudowane dane lub z pliku
# dane <- read.csv("pkb.csv")
# pkb <- ts(dane$pkb, start = c(2000, 1), frequency = 4)  # kwartalne

# Dla ilustracji — symulujemy realistyczny szereg PKB
set.seed(42)
n <- 80
trend <- 100 * (1.005)^(1:n)               # wzrost ~2% rocznie
sezon <- rep(c(-1.5, 0.5, 1.0, 0), length.out = n)
pkb <- ts(trend + sezon + arima.sim(list(ar = 0.6), n) * 2,
          start = c(2004, 1), frequency = 4)

# Python
import pandas as pd
import statsmodels.api as sm

# Przykład: pobierz z FRED przez pandas_datareader
# import pandas_datareader as pdr
# pkb = pdr.get_data_fred('GDPC1', start='2000-01-01')

# lub wczytaj wbudowane dane makro
dane = sm.datasets.macrodata.load_pandas().data
pkb = dane['realgdp']

Krok 2: Wizualizacja i dekompozycja

# Wykres szeregu
autoplot(pkb) +
  labs(title = "Realne PKB (kwartalne)", y = "PKB", x = "Rok")

# Dekompozycja na trend + sezonowość + reszta
dekompozycja <- stl(pkb, s.window = "periodic")
autoplot(dekompozycja)

Co widzimy: rosnący trend, regularna sezonowość kwartalna oraz nieregularna składowa resztkowa — struktura typowa dla szeregów makroekonomicznych.

Krok 3: Test stacjonarności

ARIMA wymaga stacjonarności szeregu. PKB jest niemal zawsze niestacjonarne (zawiera trend). Weryfikacji stacjonarności dokonuje się formalnie za pomocą testów statystycznych:

library(tseries)
adf.test(pkb)        # Augmented Dickey-Fuller
# H0: szereg NIEstacjonarny (jest pierwiastek jednostkowy)
# p > 0.05 -> nie odrzucamy -> niestacjonarny -> trzeba różnicować

from statsmodels.tsa.stattools import adfuller

wynik = adfuller(pkb)
print(f'ADF: {wynik[0]:.3f}, p-value: {wynik[1]:.3f}')
# p > 0.05 -> niestacjonarny

Różnicowanie

# Pierwsza różnica (zmiana kwartał do kwartału)
pkb_diff <- diff(pkb)
adf.test(pkb_diff)   # teraz p < 0.05 -> stacjonarny

# Ile różnicowań potrzeba? (automatycznie)
library(forecast)
ndiffs(pkb)          # zwykle 1
nsdiffs(pkb)         # różnicowanie sezonowe

Krok 4: Identyfikacja modelu (ACF/PACF)

Wykresy autokorelacji służą do identyfikacji rzędów AR i MA:

ggtsdisplay(pkb_diff)   # szereg + ACF + PACF razem

# Reguły:
# - ACF urywa się po q, PACF wygasa -> MA(q)
# - PACF urywa się po p, ACF wygasa -> AR(p)
# - oba wygasają -> ARMA(p,q)

Wzorzec	Model
ACF urywa się po lag q	MA(q)
PACF urywa się po lag p	AR(p)
Oba stopniowo wygasają	ARMA(p,q)

Krok 5: Dobór modelu — auto.arima

W praktyce stosuje się automatyczny dobór modelu minimalizujący kryterium AIC:

model <- auto.arima(pkb,
                    seasonal = TRUE,    # uwzględnij sezonowość (SARIMA)
                    stepwise = FALSE,   # pełne przeszukanie
                    approximation = FALSE)
summary(model)
# Np. ARIMA(1,1,1)(0,1,1)[4] — model SARIMA dla danych kwartalnych

import pmdarima as pm

model = pm.auto_arima(pkb, seasonal=True, m=4,
                      stepwise=True, suppress_warnings=True,
                      trace=True)
print(model.summary())

Krok 6: Diagnostyka reszt

Poprawnie dobrany model charakteryzuje się resztami stanowiącymi biały szum (brak autokorelacji):

checkresiduals(model)
# Sprawdza: wykres reszt, ACF reszt, test Ljunga-Boxa
# H0 (Ljung-Box): reszty to biały szum -> chcemy p > 0.05

model.plot_diagnostics(figsize=(10, 8))
# Q-Q plot, ACF reszt, histogram — reszty powinny być ~normalne i nieskorelowane

Jeśli reszty wykazują autokorelację → model jest niedopasowany; należy zwiększyć rząd lub zmienić specyfikację.

Krok 7: Walidacja na danych testowych

Dane dzielimy na dwie części — estymacja na obserwacjach wcześniejszych, ocena na późniejszych:

# Podział: ostatnie 8 kwartałów jako test
train <- window(pkb, end = c(2021, 4))
test  <- window(pkb, start = c(2022, 1))

model_tr <- auto.arima(train)
prognoza_test <- forecast(model_tr, h = length(test))

# Miary dokładności
accuracy(prognoza_test, test)
# RMSE, MAE, MAPE — im mniejsze tym lepiej

Miara	Interpretacja
RMSE	Pierwiastek błędu średniokwadratowego (kara za duże błędy)
MAE	Średni błąd bezwzględny
MAPE	Średni błąd procentowy (porównywalny między szeregami)

Krok 8: Prognoza końcowa

Po walidacji estymujemy na wszystkich danych i prognozujemy w przyszłość:

# Prognoza na 8 kwartałów (2 lata)
model_final <- auto.arima(pkb)
prognoza <- forecast(model_final, h = 8, level = c(80, 95))

autoplot(prognoza) +
  labs(title = "Prognoza PKB na 8 kwartałów",
       y = "PKB", x = "Rok") +
  theme_minimal()

print(prognoza)

# Python — prognoza z przedziałami
prognoza, ci = model.predict(n_periods=8, return_conf_int=True)
print(prognoza)
print(ci)   # przedziały ufności

        Point Forecast   Lo 80   Hi 80   Lo 95   Hi 95
2024 Q1         128.4   126.1   130.7   124.9   131.9
2024 Q2         129.7   126.8   132.6   125.2   134.2
...

Przedziały ufności rosną wraz z horyzontem — im dalej sięga prognoza, tym większa niepewność. Jest to zjawisko naturalne i oczekiwane.

Wnioski i dobre praktyki

Stacjonarność należy zawsze testować przed dopasowaniem ARIMA (ADF, KPSS)
Szeregi z trendem wymagają różnicowania; należy uwzględnić sezonowość (SARIMA)
auto.arima stanowi dobry punkt wyjścia, lecz reszty wymagają sprawdzenia
Walidacja out-of-sample jest konieczna — dopasowanie in-sample bywa mylące
Przedziały ufności należy raportować — prognoza punktowa jest niewystarczająca
Załamania strukturalne (kryzysy, pandemia) wymagają ostrożności — ARIMA ich nie przewiduje

Ograniczenie: modele jednorównaniowe (ARIMA) nie uwzględniają zależności od innych zmiennych. Gdy PKB zależy od stóp procentowych, inflacji itp. → warto rozważyć modele VAR.

Powiązane: Szeregi czasowe — teoria · Modele VAR · Regresja w praktyce