Autokorelácia rezíduí

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

1. Načítanie dát a predspracovanie

eco <- read.csv("economics.csv", stringsAsFactors = FALSE)
eco$date <- as.Date(eco$date)
str(eco)

'data.frame':   574 obs. of  7 variables:
 $ X       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ date    : Date, format: "1967-07-01" "1967-08-01" "1967-09-01" ...
 $ pce     : num  507 510 516 512 517 ...
 $ pop     : num  198712 198911 199113 199311 199498 ...
 $ psavert : num  12.6 12.6 11.9 12.9 12.8 11.8 11.7 12.3 11.7 12.3 ...
 $ uempmed : num  4.5 4.7 4.6 4.9 4.7 4.8 5.1 4.5 4.1 4.6 ...
 $ unemploy: int  2944 2945 2958 3143 3066 3018 2878 3001 2877 2709 ...

summary(eco)

       X              date                 pce               pop            psavert      
 Min.   :  1.0   Min.   :1967-07-01   Min.   :  506.7   Min.   :198712   Min.   : 2.200  
 1st Qu.:144.2   1st Qu.:1979-06-08   1st Qu.: 1578.3   1st Qu.:224896   1st Qu.: 6.400  
 Median :287.5   Median :1991-05-16   Median : 3936.8   Median :253060   Median : 8.400  
 Mean   :287.5   Mean   :1991-05-17   Mean   : 4820.1   Mean   :257160   Mean   : 8.567  
 3rd Qu.:430.8   3rd Qu.:2003-04-23   3rd Qu.: 7626.3   3rd Qu.:290291   3rd Qu.:11.100  
 Max.   :574.0   Max.   :2015-04-01   Max.   :12193.8   Max.   :320402   Max.   :17.300  
    uempmed          unemploy    
 Min.   : 4.000   Min.   : 2685  
 1st Qu.: 6.000   1st Qu.: 6284  
 Median : 7.500   Median : 7494  
 Mean   : 8.609   Mean   : 7771  
 3rd Qu.: 9.100   3rd Qu.: 8686  
 Max.   :25.200   Max.   :15352  

pce_ts <- ts(eco$pce, start = c(as.numeric(format(min(eco$date), "%Y")), as.numeric(format(min(eco$date), "%m"))), frequency = 12)


plot(pce_ts, main = "PCE (Personal Consumption Expenditures) - časová rada", ylab = "pce")

2. Jednoduchá lineárna regresia

Ako príklad odhadnem model, kde závislá premenná je pce a vysvetľujúce premenné sú unemploy (počet nezamestnaných) a psavert (osobná miera úspor):

model <- lm(pce ~ unemploy + psavert, data = eco)
summary(model)

Call:
lm(formula = pce ~ unemploy + psavert, data = eco)

Residuals:
    Min      1Q  Median      3Q     Max 
-4908.2 -1051.0  -307.4   951.0  6577.6 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7400.27144  353.89141   20.91   <2e-16 ***
unemploy       0.54969    0.02783   19.75   <2e-16 ***
psavert     -799.79277   24.80336  -32.24   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1674 on 571 degrees of freedom
Multiple R-squared:  0.7794,    Adjusted R-squared:  0.7786 
F-statistic:  1009 on 2 and 571 DF,  p-value: < 2.2e-16

# uložím rezidua
res <- residuals(model)


# graf rezíduí
ts.plot(res, main = "Rezíduá z lineárneho modelu", ylab = "rezíduá")

3. Diagnostika rezíduí — autokorelácia

par(mfrow = c(1,2))
acf(res, main = "ACF rezíduí")
pacf(res, main = "PACF rezíduí")

par(mfrow = c(1,1))

Ljung–Box test

# test na náhodnosť rezíduí (Ljung-Box)
Box.test(res, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  res
X-squared = 4038.3, df = 12, p-value < 2.2e-16

Durbin–Watson test (prvá kontrola autokorelácie)

library(lmtest)
dwtest(model)

    Durbin-Watson test

data:  model
DW = 0.12773, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

Breusch–Godfrey test (viacnásobná autokorelácia)

# bgtest z balíka lmtest - test na autokorelaciu rádov 1..p
bgtest(model, order = 12) # testujeme autokoreláciu až do 12 periód

    Breusch-Godfrey test for serial correlation of order up to 12

data:  model
LM test = 509.81, df = 12, p-value < 2.2e-16

4. Robustné štandardné chyby — Newey–West

library(sandwich)
library(lmtest)
coeftest(model, vcov = NeweyWest(model, lag = 12, prewhite = FALSE))

t test of coefficients:

               Estimate  Std. Error  t value  Pr(>|t|)    
(Intercept) 7400.271439  906.339604   8.1650 2.068e-15 ***
unemploy       0.549694    0.093081   5.9056 6.038e-09 ***
psavert     -799.792766   57.529831 -13.9022 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5. Alternatíva: dynamický model (lagy závislej premennej)

# pridám lag 1 z pce
eco$pce_lag1 <- dplyr::lag(eco$pce, 1)
# odstránime prvé NA pozorovanie
eco_dyn <- na.omit(eco[, c("date", "pce", "pce_lag1", "unemploy", "psavert")])


model_dyn <- lm(pce ~ pce_lag1 + unemploy + psavert, data = eco_dyn)
summary(model_dyn)

Call:
lm(formula = pce ~ pce_lag1 + unemploy + psavert, data = eco_dyn)

Residuals:
     Min       1Q   Median       3Q      Max 
-169.204   -7.954    0.067    9.415  159.690 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 39.8107519  6.8857547    5.782 1.22e-08 ***
pce_lag1     1.0004423  0.0006150 1626.833  < 2e-16 ***
unemploy     0.0003373  0.0005312    0.635    0.526    
psavert     -2.8229677  0.6108563   -4.621 4.72e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 24.56 on 569 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 3.995e+06 on 3 and 569 DF,  p-value: < 2.2e-16

# test autokorelácie rezíduí modelu_dyn
res_dyn <- residuals(model_dyn)
acf(res_dyn, main = "ACF rezíduí - dynamický model")

Box.test(res_dyn, lag = 12, type = "Ljung-Box")

    Box-Ljung test

data:  res_dyn
X-squared = 32.551, df = 12, p-value = 0.001137

dwtest(model_dyn)

    Durbin-Watson test

data:  model_dyn
DW = 2.0407, p-value = 0.6428
alternative hypothesis: true autocorrelation is greater than 0

6. Záver a komentár

V modeli, kde bola závislou premennou spotreba pce a vysvetľujúcimi premennými unemploy a psavert, sme dostali veľmi vysokú štatistickú významnosť oboch koeficientov a relatívne vysoké R² (≈ 0.78). Samotné odhady však nie je možné považovať za spoľahlivé, pretože diagnostika rezíduí jednoznačne odhalila silnú pozitívnu autokoreláciu:

• ACF aj PACF grafy rezíduí ukázali dlhú pomalú decay (typický prejav trendovej časovej rady)
• Ljung–Box test (p < 2.2e-16) autokoreláciu potvrdil
• Durbin–Watson štatistika (DW = 0.1277) je extrémne vzdialená od 2 → výrazná pozitívna autokorelácia
• Breusch–Godfrey test pre 12 periód taktiež zamietol nezávislosť rezíduí (p < 2.2e-16)

Autokorelácia rezíduí znamená, že OLS odhady síce zostávajú neskreslené, ale ich štandardné chyby sú podhodnotené, čo môže viesť k mylným záverom o štatistickej významnosti koeficientov.

Robustné štandardné chyby – Newey–West

Po aplikácii Newey–West korekcie došlo k výraznému zvýšeniu štandardných chýb pri všetkých koeficientoch, čo je očakávané pri súčasnej autokorelácii aj heteroskedasticite.

• premenná unemploy zostala štatisticky významná, ale jej t-hodnota sa znížila,
• premenná psavert zostala vysoko významná,

To znamená, že model je aj po korekcii štatisticky silný, ale pôvodný bol príliš „optimistický“.

Dynamický model

Po pridaní jedného oneskorenia závislej premennej (pce_lag1) sa situácia zásadne zmenila:

• koeficient pri pce_lag1 ≈ 1.0004 → spotreba je silne autoregresívna,
• koeficient pri unemploy už nie je štatisticky významný (p = 0.526),
• psavert zostáva významná, ale jej efekt je menší než v prvom modeli,
• rezíduá sú dramaticky menšie (štandardná chyba prudko klesla z 1674 → 24.56).

Diagnostika rezíduí dynamického modelu:

• ACF ukazuje už len mierne korelácie,
• Ljung–Box test (p = 0.0011) naznačuje slabú autokoreláciu,
• Durbin–Watson (≈ 2.04, p = 0.64).

To znamená, že dynamický model veľmi dobre odstránil autokoreláciu.

Interpretácia výsledkov

• Pôvodný statický model bol síce štatisticky významný, ale trpel veľmi silnou autokoreláciou rezíduí, čo spôsobovalo skreslené závery.
• Newey–West korekcia poskytla realistickejšie štandardné chyby.
• Dynamický model s oneskorenou premennou pce_lag1:
  • výrazne zlepšil kvalitu modelu,
  • odstránil väčšinu autokorelácie,
  • ukázal, že spotreba má takmer jednotkovú autoregresiu,
  • zmenil významnosť premenných (unemploy prestal byť významný).

Z ekonomického hľadiska to znamená, že minulá spotreba je dominantným prediktorom dnešnej spotreby, čo je typické pre makroekonomické časové rady.

Zhrnutie

Autokorelácia bola v pôvodnom modeli extrémne silná, čo potvrdili všetky testy. Po dynamizácii modelu sa autokorelácia odstránila a nový model poskytuje spoľahlivejšie výsledky. Najlepší a najstabilnejší model je preto dynamický model s pce_lag1.