Autokorelacja składnika losowego — diagnoza, dowody i korekta

Średniozaawansowany
Streszczenie

Autokorelacja od podstaw: model AR(1) i źródła, dowód, że estymator MNK traci efektywność i zaniża błędy standardowe (wariancja kanapkowa z macierzą Toeplitza), test Durbina-Watsona z wyprowadzeniem relacji DW≈2(1−ρ̂), test Breuscha-Godfreya, błędy HAC Neweya-Westa, uogólniona MNK i transformacja Cochrane'a-Orcutta jako BLUE — z wykresami w skali szarości, szczegółowo policzonymi przykładami i kodem R oraz Python.

Drugą połową założenia sferyczności klasycznego modelu MNK jest brak korelacji między składnikami losowymi różnych obserwacji. Jego naruszenie — autokorelacja, zwana też korelacją szeregową — jest typowe dla szeregów czasowych. Rozważymy jej skutki dla estymatora, sposoby wykrycia oraz korektę, dowodząc kluczowych własności i ilustrując je rachunkiem.

Definicja i model AR(1)

Definicja
Autokorelacja i proces AR(1)

Autokorelacja występuje, gdy składniki losowe różnych obserwacji są skorelowane:

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

Najczęstszą postacią jest autoregresja pierwszego rzędu AR(1):

$$ \varepsilon_t=\rho\,\varepsilon_{t-1}+u_t,\qquad |\rho|<1,\quad u_t\ \text{biały szum}. $$

Dla $\rho>0$ mówimy o autokorelacji dodatniej (inercyjnej, typowej w ekonomii), dla $\rho<0$ — ujemnej. Wówczas $\mathrm{Cov}(\varepsilon_t,\varepsilon_{t-h})=\dfrac{\sigma_u^2}{1-\rho^2}\,\rho^{|h|}$ — korelacja maleje geometrycznie z odległością w czasie.

Macierz kowariancji błędów AR(1) ma więc strukturę Toeplitza o wyrazach $\rho^{|i-j|}$: jedynki na przekątnej i geometrycznie zanikające wartości poza nią.

Powierzchnia macierzy kowariancji AR(1)
Macierz korelacji procesu AR(1): $\mathrm{Corr}(\varepsilon_i,\varepsilon_j)=\rho^{|i-j|}$ dla $\rho=0{,}6$. Na przekątnej wartość $1$; im dalej od przekątnej, tym słabsza korelacja. Przy homoskedastyczności bez autokorelacji powierzchnia byłaby zerowa wszędzie poza przekątną — tutaj jej grzbiet wzdłuż przekątnej rozlewa się na sąsiednie obserwacje.

Najczęstsze źródła autokorelacji to inercja zjawisk (PKB, inflacja zmieniają się płynnie), pominięta zmienna o przebiegu czasowym (trend, cykl trafia do składnika losowego), błędna forma funkcyjna oraz wygładzanie danych.

Konsekwencje dla estymatora MNK

Twierdzenie
Skutki autokorelacji

Przy autokorelacji estymator MNK $\hat{\boldsymbol\beta}$ pozostaje nieobciążony i zgodny (o ile regresory są egzogeniczne), lecz:

  1. traci efektywność — nie jest już BLUE;
  2. jego prawdziwa wariancja to macierz kanapkowa $(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega\,\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}$ z niediagonalną $\boldsymbol\Omega$;
  3. klasyczne błędy standardowe są błędne — przy autokorelacji dodatniej zwykle zaniżone, co zawyża statystyki $t$ i prowadzi do pozornej istotności.
Dowód
Dlaczego dodatnia autokorelacja zaniża błędy standardowe
  1. Postać wariancji. Jak w dowodzie dla [heteroskedastyczności](/ekonometria/heteroskedastycznosc/), z $\hat{\boldsymbol\beta}-\boldsymbol\beta=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\varepsilon$ i $\mathrm{Var}(\boldsymbol\varepsilon)=\boldsymbol\Omega$ otrzymujemy $$ \mathrm{Var}(\hat{\boldsymbol\beta}\mid\mathbf{X})=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega\,\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}. $$ Nieobciążoność zachowana, bo dowód korzysta tylko z $\mathbb{E}[\boldsymbol\varepsilon\mid\mathbf{X}]=\mathbf{0}$.
  2. Wyrazy pozadiagonalne nie znikają. Dla AR(1) $\boldsymbol\Omega_{ij}=\tfrac{\sigma_u^2}{1-\rho^2}\rho^{|i-j|}$. Środkowy czynnik $$ \mathbf{X}^\top\boldsymbol\Omega\mathbf{X}=\frac{\sigma_u^2}{1-\rho^2}\sum_i\sum_j \rho^{|i-j|}\,\mathbf{x}_i\mathbf{x}_j^\top $$ zawiera dodatnie składniki krzyżowe ($i\neq j$), których wzór klasyczny $\sigma^2\mathbf{X}^\top\mathbf{X}$ pomija.
  3. Znak obciążenia. Gdy regresory są same dodatnio autoskorelowane (np. trend), iloczyny $\mathbf{x}_i\mathbf{x}_j^\top$ dla bliskich $i,j$ są dodatnie, a $\rho^{|i-j|}>0$. Składniki krzyżowe **powiększają** prawdziwą wariancję ponad $\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$. Klasyczny estymator, ignorując je, **zaniża** wariancję — stąd zawyżone $t$ i pozorna istotność.
  4. Utrata efektywności. Skoro $\boldsymbol\Omega\neq\sigma^2\mathbf{I}$, [twierdzenie Gaussa-Markowa](/ekonometria/zalozenia-kmnk/) nie gwarantuje już optymalności MNK; estymator o mniejszej wariancji daje GLS.
∎

Diagnostyka graficzna

Pierwszym krokiem jest wykres reszt w czasie. Przy autokorelacji dodatniej reszty układają się w długie fale (po wartości dodatniej następuje kolejna dodatnia); przy braku autokorelacji zmieniają znak chaotycznie.

Reszty w czasie — autokorelacja dodatnia a brak autokorelacji
Reszty w funkcji czasu. Po lewej dodatnia autokorelacja: gładkie fale, długie serie tego samego znaku. Po prawej brak autokorelacji: szybka, nieregularna zmiana znaku. Wzrokowa ocena gładkości przebiegu reszt jest pierwszym sygnałem korelacji szeregowej.

Uzupełnieniem jest wykres reszty bieżącej względem opóźnionej; dodatnie nachylenie chmury punktów jest bezpośrednim obrazem współczynnika $\hat\rho$.

Reszta bieżąca względem opóźnionej — nachylenie rho
Reszta $\hat\varepsilon_t$ względem $\hat\varepsilon_{t-1}$. Punkty układające się wzdłuż prostej o dodatnim nachyleniu oznaczają, że duża reszta pociąga za sobą dużą resztę w następnym okresie — to próbkowy współczynnik autokorelacji $\hat\rho\approx 0{,}51$.

Test Durbina-Watsona

Twierdzenie
Statystyka Durbina-Watsona

Test klasyczny 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) $$

przyjmuje wartości w $[0,4]$: $DW\approx 2$ oznacza brak autokorelacji, $DW\to 0$ — silną dodatnią, $DW\to 4$ — silną ujemną. Rozkład zależy od $\mathbf{X}$, więc używa się tablicowych granic $d_L,d_U$: gdy $DWd_U$ nie odrzucamy, między nimi rozstrzygnięcie jest niemożliwe.

Dowód
Relacja DW ≈ 2(1 − ρ̂)
  1. Rozwinięcie licznika. $\sum_{t=2}^n(\hat\varepsilon_t-\hat\varepsilon_{t-1})^2=\sum_{t=2}^n\hat\varepsilon_t^2+\sum_{t=2}^n\hat\varepsilon_{t-1}^2-2\sum_{t=2}^n\hat\varepsilon_t\hat\varepsilon_{t-1}$.
  2. Przybliżenie sum brzegowych. Dla dużego $n$ obie sumy kwadratów różnią się tylko o jeden skrajny wyraz od pełnej sumy: $\sum_{t=2}^n\hat\varepsilon_t^2\approx\sum_{t=2}^n\hat\varepsilon_{t-1}^2\approx\sum_{t=1}^n\hat\varepsilon_t^2$.
  3. Podstawienie definicji $\hat\rho$. Ponieważ $\hat\rho=\dfrac{\sum_{t=2}^n\hat\varepsilon_t\hat\varepsilon_{t-1}}{\sum_{t=1}^n\hat\varepsilon_t^2}$, dzieląc licznik DW przez $\sum_t\hat\varepsilon_t^2$: $$ DW\approx\frac{\sum\hat\varepsilon_t^2+\sum\hat\varepsilon_t^2-2\sum\hat\varepsilon_t\hat\varepsilon_{t-1}}{\sum\hat\varepsilon_t^2}=2-2\hat\rho=2(1-\hat\rho). $$ Przybliżenie jest dokładne asymptotycznie; w małej próbie wyrazy brzegowe powodują rozbieżność między $DW$ a $2(1-\hat\rho)$.
∎
Strefy decyzyjne testu Durbina-Watsona
Strefy decyzyjne testu Durbina-Watsona dla $n=12$ i jednego regresora ($d_L\approx 0{,}97$, $d_U\approx 1{,}33$). Ciemne pasy to obszary odrzucenia (autokorelacja dodatnia po lewej, ujemna po prawej), jasne — strefy nierozstrzygnięcia, środek — brak autokorelacji. Pionowa linia zaznacza wartość $DW=0{,}64$ z przykładu, leżącą w obszarze odrzucenia.
Przykład
Durbin-Watson: autokorelacja dodatnia

Dla regresji $y$ na czas $t$ ($n=12$, $\hat y=1{,}0+0{,}5\,t$) reszty tworzą gładką falę: $\hat\varepsilon=(-1{,}44;\ -0{,}10;\ 0{,}90;\ 1{,}33;\ 1{,}12;\ 0{,}43;\ -0{,}43;\ -1{,}12;\ -1{,}33;\ -0{,}90;\ 0{,}10;\ 1{,}44)$. Liczymy

$$ DW=\frac{\sum_{t=2}^{12}(\hat\varepsilon_t-\hat\varepsilon_{t-1})^2}{\sum_{t=1}^{12}\hat\varepsilon_t^2}=0{,}64. $$

Granice tablicowe $d_L\approx0{,}97$, $d_U\approx1{,}33$. Ponieważ $0{,}64<0{,}97=d_L$, odrzucamy $H_0$ — występuje silna autokorelacja dodatnia. Z relacji $DW\approx2(1-\hat\rho)$ implikowane $\hat\rho\approx1-\tfrac{0{,}64}{2}=0{,}68$; bezpośrednia korelacja reszt z opóźnieniem wynosi $0{,}51$ (różnica wynika z asymptotycznego charakteru wzoru i wyrazów brzegowych przy $n=12$).

Przykład
Durbin-Watson: autokorelacja ujemna

Dla reszt zmieniających znak niemal naprzemiennie, $\hat\varepsilon=(0{,}9;\ -1{,}2;\ 0{,}6;\ -0{,}8;\ 1{,}1;\ -0{,}5;\ 0{,}7;\ -1{,}0;\ 0{,}4;\ -0{,}9;\ 1{,}0;\ -0{,}3)$, otrzymujemy

$$ DW=3{,}52,\qquad \hat\rho=-0{,}81. $$

Wartość $DW=3{,}52$ przekracza $4-d_L\approx 3{,}03$, więc odrzucamy $H_0$ na rzecz autokorelacji ujemnej — kolejne reszty systematycznie zmieniają znak.

Test Breuscha-Godfreya

Test Durbina-Watsona wykrywa wyłącznie AR(1) i zawodzi przy opóźnionej zmiennej objaśnianej. Ogólniejszy jest test Breuscha-Godfreya.

Twierdzenie
Test Breuscha-Godfreya (LM)

Estymujemy regresję pomocniczą reszt na regresory modelu oraz $p$ opóźnionych reszt:

$$ \hat\varepsilon_t=\mathbf{x}_t^\top\boldsymbol\gamma+\rho_1\hat\varepsilon_{t-1}+\dots+\rho_p\hat\varepsilon_{t-p}+u_t. $$

Statystyka $\mathrm{LM}=(n-p)\,R^2_{\text{aux}}\sim\chi^2_p$ przy $H_0$ (brak autokorelacji do rzędu $p$). Test obejmuje dowolny rząd i dopuszcza opóźnione regresory, dlatego jest preferowany.

Przykład
Breusch-Godfrey rzędu 1

Dla danych z przykładu 1 regresja pomocnicza $\hat\varepsilon_t$ na $\{1,\ t,\ \hat\varepsilon_{t-1}\}$ (dla $t=2,\dots,12$, czyli $n-1=11$ obserwacji) daje $R^2_{\text{aux}}=0{,}429$, zatem

$$ \mathrm{LM}=(n-1)\,R^2_{\text{aux}}=11\cdot 0{,}429=4{,}72,\qquad \chi^2_{1;\,0{,}05}=3{,}84. $$

Ponieważ $4{,}72>3{,}84$ ($p=0{,}030$), odrzucamy brak autokorelacji — zgodnie z testem Durbina-Watsona.

Korekta I: błędy standardowe HAC (Newey-West)

Twierdzenie
Estymator HAC

Estymator Neweya-Westa szacuje środek macierzy kanapkowej, uwzględniając autokowariancje reszt do rzędu $L$ z malejącymi wagami $w_\ell=1-\dfrac{\ell}{L+1}$ (Bartletta):

$$ \widehat{\mathbf{S}}=\sum_t\hat\varepsilon_t^2\mathbf{x}_t\mathbf{x}_t^\top+\sum_{\ell=1}^{L}w_\ell\sum_{t=\ell+1}^{n}\hat\varepsilon_t\hat\varepsilon_{t-\ell}\big(\mathbf{x}_t\mathbf{x}_{t-\ell}^\top+\mathbf{x}_{t-\ell}\mathbf{x}_t^\top\big). $$

Wagi Bartletta zapewniają dodatnią półokreśloność macierzy. Błędy HAC są odporne jednocześnie na autokorelację i heteroskedastyczność i nie wymagają znajomości jej struktury — stąd ich rola domyślnego narzędzia w analizie empirycznej.

Podejście to jest autokorelacyjnym uogólnieniem odpornych błędów White’a: nie zmienia estymatora $\hat{\boldsymbol\beta}$, lecz koryguje jego macierz wariancji, naprawiając wnioskowanie.

Korekta II: GLS i transformacja Cochrane’a-Orcutta

Gdy struktura jest znana (AR(1)), efektywność odzyskuje się przez transformację quasi-różnicową, która „wybiela" błąd.

Twierdzenie
Transformacja Cochrane'a-Orcutta

Dla modelu $y_t=\beta_0+\beta_1 x_t+\varepsilon_t$ z błędem AR(1) ($\varepsilon_t=\rho\varepsilon_{t-1}+u_t$) przekształcony model

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

ma już biały składnik losowy $u_t$, więc MNK na danych przekształconych jest GLS, a zatem BLUE. W praktyce $\rho$ estymuje się iteracyjnie (Cochrane-Orcutt) albo łącznie z $\boldsymbol\beta$ (Prais-Winsten zachowuje pierwszą obserwację).

Dowód
Quasi-różnicowanie usuwa autokorelację
  1. Zapis modelu w dwóch okresach. Mamy $y_t=\beta_0+\beta_1 x_t+\varepsilon_t$ oraz $y_{t-1}=\beta_0+\beta_1 x_{t-1}+\varepsilon_{t-1}$.
  2. Odejmujemy $\rho$-krotność opóźnienia. $$ y_t-\rho y_{t-1}=\beta_0(1-\rho)+\beta_1(x_t-\rho x_{t-1})+(\varepsilon_t-\rho\varepsilon_{t-1}). $$
  3. Błąd staje się biały. Z definicji AR(1) $\varepsilon_t-\rho\varepsilon_{t-1}=u_t$ — szum o stałej wariancji i bez autokorelacji. Przekształcony model spełnia więc założenia klasyczne, a MNK na nim jest BLUE. To szczególny przypadek [GLS](/ekonometria/heteroskedastycznosc/) z $\boldsymbol\Omega$ o strukturze AR(1).
∎
Przykład
Jedna iteracja Cochrane'a-Orcutta

Dla danych z przykładu 1 estymowane $\hat\rho=0{,}51$ (z reszt MNK). Quasi-różnicujemy zmienne, $\tilde y_t=y_t-0{,}51\,y_{t-1}$ oraz $\tilde x_t=t-0{,}51(t-1)$, i estymujemy MNK na przekształconych danych:

$$ \hat\beta_1^{\text{CO}}=0{,}42\quad\text{wobec}\quad \hat\beta_1^{\text{MNK}}=0{,}50. $$

Oszacowanie nachylenia zmienia się, a — co istotniejsze — jego błąd standardowy zostaje policzony przy już niezautokorelowanym błędzie, więc jest poprawny. Procedurę powtarza się (nowe reszty → nowe $\hat\rho$ → ponowna transformacja) aż do zbieżności.

Najpierw specyfikacja, potem korekta.

Autokorelacja bywa objawem złej specyfikacji, a nie samodzielnym problemem. Pominięty trend lub sezonowość trafiają do składnika losowego, tworząc pozorną korelację szeregową. Zanim sięgniemy po HAC czy GLS, należy sprawdzić, czy modelu nie naprawia dodanie trendu, zmiennych sezonowych, opóźnień (model dynamiczny ADL) lub różnicowanie niestacjonarnego szeregu. Korekta błędów standardowych nie usuwa obciążenia ze złej specyfikacji.

Podsumowanie

Definicja
Autokorelacja w pigułce
  • Definicja: $\mathrm{Cov}(\varepsilon_i,\varepsilon_j)\neq0$; typowo AR(1) z kowariancją $\propto\rho^{|i-j|}$.
  • Skutki: MNK nieobciążony i zgodny, lecz nie-BLUE; klasyczne SE błędne (dla $\rho>0$ zaniżone → zawyżone $t$).
  • Diagnoza: wykres reszt w czasie, test Durbina-Watsona ($DW\approx2(1-\hat\rho)$) z granicami $d_L,d_U$, ogólniejszy test Breuscha-Godfreya ($\mathrm{LM}=(n-p)R^2_{\text{aux}}$).
  • Leczenie: błędy HAC Neweya-Westa (bez założeń o strukturze) albo GLS / Cochrane-Orcutt (efektywne, gdy znamy AR(1)).
  • Granica: najpierw poprawna specyfikacja (trend, sezonowość, opóźnienia), dopiero potem korekta.

Dalej: Heteroskedastyczność · Szeregi czasowe · Założenia KMNK · Testy hipotez

Literatura uzupełniająca
  • Podręcznik: J. M. Wooldridge, Introductory Econometrics: A Modern Approach, rozdz. 12
  • Podręcznik: W. H. Greene, Econometric Analysis, rozdz. 20 (modele szeregów czasowych)
  • J. Durbin, G. S. Watson, „Testing for Serial Correlation in Least Squares Regression", Biometrika 37–38, 1950–1951
  • W. K. Newey, K. D. West, „A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix", Econometrica 55(3), 1987, s. 703–708
  • Pakiety R: lmtest (dwtest, bgtest), sandwich (NeweyWest), orcutt
Oprogramowanie

R:

library(lmtest); library(sandwich)
t <- 1:12
y <- c(0.059, 1.902, 3.402, 4.325, 4.619, 4.432,
       4.068, 3.881, 4.175, 5.098, 6.598, 8.441)
m <- lm(y ~ t)

dwtest(m)                 # DW = 0,64; silna dodatnia autokorelacja
bgtest(m, order = 1)      # LM = 4,72; df = 1; p = 0,030

## Bledy HAC (Newey-West) — odporne na autokorelacje i heteroskedastycznosc
coeftest(m, vcov = NeweyWest(m, lag = 2, prewhite = FALSE))

## GLS / Cochrane-Orcutt (iteracyjna estymacja rho)
library(orcutt); summary(cochrane.orcutt(m))   # beta_t ~ 0,42

Python:

import numpy as np, statsmodels.api as sm
from statsmodels.stats.stattools import durbin_watson
from statsmodels.stats.diagnostic import acorr_breusch_godfrey

t = np.arange(1, 13)
y = np.array([0.059, 1.902, 3.402, 4.325, 4.619, 4.432,
              4.068, 3.881, 4.175, 5.098, 6.598, 8.441])
m = sm.OLS(y, sm.add_constant(t)).fit()

print(durbin_watson(m.resid))                 # 0,64
print(acorr_breusch_godfrey(m, nlags=1)[:2])  # LM = 4,72; p = 0,030
print(sm.OLS(y, sm.add_constant(t)).fit(cov_type='HAC',
      cov_kwds={'maxlags': 2}).bse)           # bledy HAC