Prognozowanie PKB — tutorial krok po kroku (ARIMA)

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

Zbudujemy model prognozujący PKB na kilka kwartałów do przodu. Przejdziemy pełną ścieżkę: dane → diagnoza → model → walidacja → prognoza.

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

Krok 1: Wczytaj dane

# 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 i nieregularna reszta. Klasyczny szereg makroekonomiczny.

Krok 3: Test stacjonarności

ARIMA wymaga stacjonarności. PKB prawie zawsze jest niestacjonarne (ma trend). Testujemy formalnie:

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 pomagają dobrać rzędy 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 używamy automatycznego doboru minimalizującego 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

Dobry model ma reszty będące białym szumem (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 mają autokorelację → model niedopasowany, zwiększ rząd lub zmień specyfikację.

Krok 7: Walidacja na danych testowych

Podziel dane: trenuj na starszych, testuj na nowszych:

# 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ą z horyzontem — im dalej prognozujemy, tym większa niepewność. To naturalne i uczciwe.

Wnioski i dobre praktyki

Zawsze testuj stacjonarność przed ARIMA (ADF, KPSS)
Różnicuj szeregi z trendem; uwzględnij sezonowość (SARIMA)
auto.arima to dobry start, ale sprawdź reszty
Waliduj out-of-sample — dopasowanie in-sample myli
Raportuj przedziały ufności — punktowa prognoza to za mało
Uważaj na załamania strukturalne (kryzysy, pandemia) — ARIMA ich nie przewidzi

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

Powiązane: Szeregi czasowe — teoria · Modele VAR · Regresja od zera