Autokorelacja składnika losowego — diagnoza i korekta

Streszczenie

Autokorelacja w modelach ekonometrycznych: przyczyny, test Durbina-Watsona, Breuscha-Godfreya, błędy Neweya-Westa (HAC) i modele z korektą — R i Python

Czym jest autokorelacja

Jedno z założeń KMNK mówi, że składniki losowe są nieskorelowane między obserwacjami:

$$ \text{Cov}(\varepsilon_i, \varepsilon_j) = 0 \quad \text{dla } i \neq j $$

Autokorelacja (korelacja szeregowa) występuje, gdy to założenie jest złamane — błąd w jednym okresie zależy od błędu w innym. Najczęściej w szeregach czasowych:

$$ \varepsilon_t = \rho\, \varepsilon_{t-1} + u_t, \qquad |\rho| < 1 $$

To autokorelacja pierwszego rzędu AR(1). Gdy $\rho > 0$ — dodatnia (typowa w ekonomii), gdy $\rho < 0$ — ujemna.

Źródła autokorelacji

Przyczyna	Mechanizm
Inercja zjawisk	PKB, inflacja, bezrobocie zmieniają się płynnie
Pominięta zmienna	Pominięty trend/cykl trafia do składnika losowego
Zła forma funkcyjna	Liniowy model dla nieliniowej zależności
Wygładzanie danych	Średnie ruchome, interpolacja

Dlaczego to problem

Estymator MNK pozostaje nieobciążony i zgodny, ALE:

❌ Błędy standardowe są błędne (zwykle zaniżone) → statystyki t zawyżone
❌ Testy istotności zawodne — „istotność" tam, gdzie jej nie ma
❌ MNK nieefektywne — istnieje estymator o mniejszej wariancji (GLS)

Najgroźniejsza konsekwencja: dodatnia autokorelacja sprawia, że zmienne wyglądają na bardziej istotne niż są. Łatwo o fałszywe wnioski.

Diagnostyka: wykres reszt

Pierwszy krok — narysuj reszty w czasie i reszty względem opóźnionych reszt:

model <- lm(y ~ x1 + x2, data = dane)
reszty <- residuals(model)

par(mfrow = c(1, 2))
plot(reszty, type = "l", main = "Reszty w czasie",
     xlab = "t", ylab = expression(hat(epsilon)[t]))
abline(h = 0, col = "red", lty = 2)

plot(reszty[-length(reszty)], reszty[-1],
     main = "Reszty: t-1 vs t",
     xlab = expression(hat(epsilon)[t-1]),
     ylab = expression(hat(epsilon)[t]))

Wyraźny wzorzec (fale, układanie się wzdłuż prostej) → autokorelacja.

Test Durbina-Watsona

Klasyczny test dla AR(1). Statystyka:

$$ DW = \frac{\sum_{t=2}^{n}(\hat{\varepsilon}_t - \hat{\varepsilon}_{t-1})^2}{\sum_{t=1}^{n}\hat{\varepsilon}_t^2} \approx 2(1 - \hat{\rho}) $$

Interpretacja statystyki $DW \in [0, 4]$:

Wartość DW	$\hat{\rho}$	Wniosek
≈ 2	≈ 0	Brak autokorelacji ✅
→ 0	→ +1	Silna dodatnia autokorelacja
→ 4	→ −1	Silna ujemna autokorelacja

library(lmtest)
dwtest(model)        # H0: brak autokorelacji (rho = 0)

Ograniczenia DW: wykrywa tylko AR(1); zawodzi gdy w modelu są opóźnienia zmiennej objaśnianej (lagged Y); ma „strefę nierozstrzygnięcia".

Test Breuscha-Godfreya (lepszy)

Ogólniejszy — wykrywa autokorelację dowolnego rzędu i działa z opóźnionymi regresorami. To preferowany test.

Regresja pomocnicza reszt na regresory modelu + opóźnione reszty:

$$ \hat{\varepsilon}_t = X_t\gamma + \rho_1\hat{\varepsilon}_{t-1} + \dots + \rho_p\hat{\varepsilon}_{t-p} + u_t $$

Statystyka $LM = n \cdot R^2 \sim \chi^2_p$ pod H0 (brak autokorelacji do rzędu p).

bgtest(model, order = 1)   # AR(1)
bgtest(model, order = 4)   # do AR(4) — np. dane kwartalne

from statsmodels.stats.diagnostic import acorr_breusch_godfrey
import statsmodels.api as sm

model = sm.OLS(y, sm.add_constant(X)).fit()
lm_stat, lm_pval, f_stat, f_pval = acorr_breusch_godfrey(model, nlags=4)
print(f"LM p-value: {lm_pval:.4f}")   # < 0.05 -> autokorelacja

Korekta 1: Błędy standardowe HAC (Newey-West)

Najprostsze i najczęstsze rozwiązanie — nie zmieniamy estymatora, tylko poprawiamy błędy standardowe na odporne na autokorelację i heteroskedastyczność (HAC — Heteroskedasticity and Autocorrelation Consistent):

library(sandwich)
coeftest(model, vcov = NeweyWest(model, lag = 4))

model_hac = sm.OLS(y, sm.add_constant(X)).fit(
    cov_type='HAC', cov_kwds={'maxlags': 4}
)
print(model_hac.summary())

Zalety: prostota, nie wymaga znajomości struktury autokorelacji. To domyślny wybór w praktyce empirycznej.

Korekta 2: Uogólniona MNK (GLS / Cochrane-Orcutt)

Gdy znamy strukturę (np. AR(1)), możemy estymować efektywniej przez transformację quasi-różnicową:

$$ y_t - \rho y_{t-1} = \beta_0(1-\rho) + \beta_1(x_t - \rho x_{t-1}) + u_t $$

library(orcutt)
co <- cochrane.orcutt(model)   # iteracyjna estymacja rho i beta
summary(co)

Korekta 3: Lepsza specyfikacja modelu

Często autokorelacja to objaw, nie choroba — sygnał błędnej specyfikacji:

Dodaj opóźnienia zmiennych (model dynamiczny ADL)
Dodaj trend lub zmienne sezonowe
Sprawdź formę funkcyjną (logarytmy, kwadraty)
Rozważ różnicowanie niestacjonarnego szeregu

Zanim „naprawisz" błędy standardowe, zadaj pytanie: czy model jest dobrze specyfikowany? Pominięty trend lub sezonowość to najczęstsze ukryte źródło autokorelacji.

Procedura postępowania

📊 Narysuj reszty w czasie
🔬 Wykonaj test Breuscha-Godfreya (rzędu dopasowanego do częstotliwości danych)
❓ Jeśli autokorelacja → sprawdź specyfikację (trend, sezonowość, opóźnienia)
🛠️ Jeśli zostaje → zastosuj błędy HAC (Newey-West)
⚡ Dla efektywności → rozważ GLS/Cochrane-Orcutt

Powiązane: Heteroskedastyczność · Szeregi czasowe