Ekonometria v R – cvičenie 9 (autokorelácia rezíduí)

1. Úvod a údaje

V tomto cvičení pracujeme s databázou powerconsumption.csv, ktorá obsahuje merania po 10 minútach. Cieľom je odhadnúť lineárny model spotreby elektriny a následne overiť predpoklad nezávislosti rezíduí (autokoreláciu).

Budeme modelovať spotrebu v zóne 1 (PowerConsumption_Zone1) v závislosti od: - teploty (Temperature), - vlhkosti (Humidity), - rýchlosti vetra (WindSpeed), - difúznych tokov (GeneralDiffuseFlows, DiffuseFlows).

1.1 Pracovná hypotéza

Predpokladáme, že meteorologické podmienky a intenzita difúznych tokov súvisia so spotrebou elektriny. Znamenko vplyvu (kladné/záporné) nefixujeme „nasilu“, pretože spotreba môže reagovať rozdielne v zime a v lete (kúrenie vs. chladenie). Dôležitá je pre nás najmä štatistická významnosť a autokorelácia rezíduí v časovom rade.

1.2 Načítanie a príprava dát

rm(list=ls())

library(dplyr)
library(lubridate)
library(ggplot2)
library(lmtest)
library(sandwich)

# nacitanie dat (desatinna bodka, oddelovac ciarka)
# predpoklad: subor je v priecinku "udaje/" pod nazvom "powerconsumption.csv"
udaje <- read.csv("powerconsumption.csv", sep = ",", dec=".", header = TRUE)

# parsovanie casu
udaje <- udaje %>%
  mutate(
    Datetime = mdy_hm(Datetime)
  ) %>%
  arrange(Datetime)

# kontrola
str(udaje)

'data.frame':   52416 obs. of  9 variables:
 $ Datetime              : POSIXct, format: "2017-01-01 00:00:00" "2017-01-01 00:10:00" ...
 $ Temperature           : num  6.56 6.41 6.31 6.12 5.92 ...
 $ Humidity              : num  73.8 74.5 74.5 75 75.7 76.9 77.7 78.2 78.1 77.3 ...
 $ WindSpeed             : num  0.083 0.083 0.08 0.083 0.081 0.081 0.08 0.085 0.081 0.082 ...
 $ GeneralDiffuseFlows   : num  0.051 0.07 0.062 0.091 0.048 0.059 0.048 0.055 0.066 0.062 ...
 $ DiffuseFlows          : num  0.119 0.085 0.1 0.096 0.085 0.108 0.096 0.093 0.141 0.111 ...
 $ PowerConsumption_Zone1: num  34056 29815 29128 28229 27336 ...
 $ PowerConsumption_Zone2: num  16129 19375 19007 18361 17872 ...
 $ PowerConsumption_Zone3: num  20241 20131 19668 18899 18442 ...

head(udaje)

Poznámka: dáta sú po 10 minútach, takže pri priamom OLS modeli bude autokorelácia takmer istá. Pre prehľadnosť si údaje zhrnieme na hodinové priemery (časový rad je stále dostatočne dlhý, ale čitateľnejší).

udaje_h <- udaje %>%
  mutate(Datetime_h = floor_date(Datetime, unit = "hour")) %>%
  group_by(Datetime_h) %>%
  summarise(
    Temperature = mean(Temperature),
    Humidity = mean(Humidity),
    WindSpeed = mean(WindSpeed),
    GeneralDiffuseFlows = mean(GeneralDiffuseFlows),
    DiffuseFlows = mean(DiffuseFlows),
    PC_Zone1 = mean(PowerConsumption_Zone1),
    PC_Zone2 = mean(PowerConsumption_Zone2),
    PC_Zone3 = mean(PowerConsumption_Zone3),
    .groups = "drop"
  )

summary(udaje_h)

   Datetime_h                   Temperature        Humidity       WindSpeed      
 Min.   :2017-01-01 00:00:00   Min.   : 3.602   Min.   :12.71   Min.   :0.05467  
 1st Qu.:2017-04-01 23:45:00   1st Qu.:14.404   1st Qu.:58.32   1st Qu.:0.07817  
 Median :2017-07-01 23:30:00   Median :18.759   Median :69.82   Median :0.08550  
 Mean   :2017-07-01 23:30:00   Mean   :18.810   Mean   :68.26   Mean   :1.95949  
 3rd Qu.:2017-09-30 23:15:00   3rd Qu.:22.867   3rd Qu.:81.35   3rd Qu.:4.91533  
 Max.   :2017-12-30 23:00:00   Max.   :39.695   Max.   :94.75   Max.   :5.93367  
 GeneralDiffuseFlows  DiffuseFlows         PC_Zone1        PC_Zone2        PC_Zone3    
 Min.   :  0.019     Min.   :  0.0400   Min.   :14329   Min.   : 8686   Min.   : 6191  
 1st Qu.:  0.064     1st Qu.:  0.1242   1st Qu.:26293   1st Qu.:17017   1st Qu.:13148  
 Median :  9.947     Median :  8.2413   Median :32342   Median :20787   Median :16428  
 Mean   :182.697     Mean   : 75.0280   Mean   :32345   Mean   :21043   Mean   :17835  
 3rd Qu.:326.488     3rd Qu.:105.8833   3rd Qu.:37318   3rd Qu.:24678   3rd Qu.:21598  
 Max.   :953.350     Max.   :861.0000   Max.   :51844   Max.   :36255   Max.   :47224  

---

2. Lineárna regresia v základnom tvare

Odhadujeme model (pre zónu 1):

\[ PC\_Zone1_t = \beta_0 + \beta_1 Temperature_t + \beta_2 Humidity_t + \beta_3 WindSpeed_t + \beta_4 GeneralDiffuseFlows_t + \beta_5 DiffuseFlows_t + e_t \]

model <- lm(PC_Zone1 ~ Temperature + Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows,
            data = udaje_h)

summary(model)

Call:
lm(formula = PC_Zone1 ~ Temperature + Humidity + WindSpeed + 
    GeneralDiffuseFlows + DiffuseFlows, data = udaje_h)

Residuals:
     Min       1Q   Median       3Q      Max 
-14359.7  -4907.0   -581.5   4124.6  17464.5 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)         26818.9010   519.9802  51.577  < 2e-16 ***
Temperature           535.9826    15.6051  34.347  < 2e-16 ***
Humidity              -57.7531     5.1953 -11.116  < 2e-16 ***
WindSpeed            -149.8401    33.0769  -4.530 5.98e-06 ***
GeneralDiffuseFlows    -1.8396     0.3638  -5.057 4.34e-07 ***
DiffuseFlows            0.2148     0.6912   0.311    0.756    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6280 on 8730 degrees of freedom
Multiple R-squared:  0.2111,    Adjusted R-squared:  0.2107 
F-statistic: 467.2 on 5 and 8730 DF,  p-value: < 2.2e-16

Graf: empíria vs. fitted

udaje_h$fitted <- fitted(model)

ggplot(udaje_h, aes(x = Datetime_h, y = PC_Zone1)) +
  geom_line() +
  geom_line(aes(y = fitted), linewidth = 0.8) +
  labs(
    title = "Spotreba elektriny (Zone1): skutočné hodnoty vs. vyrovnané hodnoty",
    x = "Čas (hodina)",
    y = "Power consumption (Zone1)"
  ) +
  theme_minimal()

---

3. Autokorelácia rezíduí

V časových radoch často platí, že rezíduá nie sú nezávislé. Ak je chyba v čase \(t\) prepojená s chybou v čase \(t-1\), hovoríme o autokorelácii.

3.1 Rezíduá a ACF

res <- residuals(model)

acf(res, lag.max = 48, main = "ACF rezíduí (do 48 hodín)")

3.2 Durbin–Watsonov test

DW test je klasická kontrola autokorelácie 1. rádu.

dwtest(model)

    Durbin-Watson test

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

3.3 Breusch–Godfrey test

BG test umožňuje testovať autokoreláciu aj pre viac lagov. Pri hodinových dátach dáva zmysel skúsiť napr. 24 hodín (1 deň).

bgtest(model, order = 24)

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

data:  model
LM test = 8283.8, df = 24, p-value < 2.2e-16

---

4. Ako riešiť autokoreláciu (prakticky)

Ukážeme dve bežné možnosti: 1) dynamický model s oneskorenou závislou premennou, 2) robustné štandardné chyby (Newey–West).

4.1 Dynamizácia: lag závislej premennej (Koyckov typ uvažovania)

Pridáme \(PC\_Zone1_{t-1}\) do regresie. Tým zachytíme zotrvačnosť spotreby.

udaje_h <- udaje_h %>%
  arrange(Datetime_h) %>%
  mutate(PC_Zone1_lag1 = lag(PC_Zone1))

model_dyn <- lm(PC_Zone1 ~ Temperature + Humidity + WindSpeed + GeneralDiffuseFlows + DiffuseFlows + PC_Zone1_lag1,
                data = udaje_h)

summary(model_dyn)

Call:
lm(formula = PC_Zone1 ~ Temperature + Humidity + WindSpeed + 
    GeneralDiffuseFlows + DiffuseFlows + PC_Zone1_lag1, data = udaje_h)

Residuals:
     Min       1Q   Median       3Q      Max 
-10137.3  -1555.1   -320.4   1366.6  10766.9 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          2.927e+03  2.458e+02  11.907  < 2e-16 ***
Temperature          4.413e+01  6.923e+00   6.373 1.94e-10 ***
Humidity            -1.301e+01  2.172e+00  -5.991 2.17e-09 ***
WindSpeed            5.274e+00  1.378e+01   0.383    0.702    
GeneralDiffuseFlows  8.011e-01  1.518e-01   5.276 1.35e-07 ***
DiffuseFlows         3.796e+00  2.880e-01  13.182  < 2e-16 ***
PC_Zone1_lag1        8.976e-01  4.393e-03 204.335  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2612 on 8728 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.8636,    Adjusted R-squared:  0.8635 
F-statistic:  9210 on 6 and 8728 DF,  p-value: < 2.2e-16

dwtest(model_dyn)

    Durbin-Watson test

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

bgtest(model_dyn, order = 24)

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

data:  model_dyn
LM test = 7770.8, df = 24, p-value < 2.2e-16

Poznámka: pri dynamickom modeli sa často zmení významnosť „počasných“ premenných, lebo časť variability vysvetlí samotná zotrvačnosť spotreby.

4.2 Newey–West robustné štandardné chyby

Ak chceme ponechať pôvodný OLS odhad koeficientov, ale opraviť štandardné chyby pre autokoreláciu a heteroskedasticitu, použijeme Newey–West.

Pri hodinových dátach je bežná voľba napr. lag 24 (1 deň).

coeftest(model, vcov = NeweyWest(model, lag = 24, prewhite = FALSE))

t test of coefficients:

                       Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)         26818.90096   885.12630 30.2995 < 2.2e-16 ***
Temperature           535.98259    25.81092 20.7657 < 2.2e-16 ***
Humidity              -57.75314     9.35556 -6.1731 6.993e-10 ***
WindSpeed            -149.84009    56.57687 -2.6484  0.008101 ** 
GeneralDiffuseFlows    -1.83957     0.40275 -4.5675 5.003e-06 ***
DiffuseFlows            0.21485     0.79960  0.2687  0.788169    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

5. Záver

Na dátach spotreby elektriny (časový rad) je autokorelácia rezíduí prirodzená a očakávaná. Formálne testy (DW, BG) a ACF graf zvyčajne potvrdia, že predpoklad nezávislosti rezíduí v základnom OLS modeli nie je splnený.

Riešenia: - dynamický model (zachytí zotrvačnosť), - alebo robustné odhady štandardných chýb (Newey–West), ak nám ide hlavne o korektnú inferenciu.