Econometrics in R - cvičenie 9

S využitím databázy Databaza.csv, ktorá obsahuje počty obyvateľov slovenských obcí podľa okresu a mesiaca v roku 2024.

Pri ďalšej práci budeme používať knižnice

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

---

1. Úvod a údaje

Údaje o počte obyvateľov sú uložené v súbore Databaza.csv.
Každý riadok predstavuje jednu obec, stĺpce:

• Okres – názov okresu,
• Obec – názov obce,
• 2024M01 … 2024M12 – počet obyvateľov v danom mesiaci roku 2024.

Súbor predpokladáme uložený v podpriečinku udaje v pracovnom priečinku (aby boli údaje oddelené od zvyšku projektu).

1.1 Úvod do problému, stanovenie hypotéz

Rozhodli sme sa sledovať vývoj počtu obyvateľov v čase vo vybranej obci.
Konkrétne budeme modelovať mesačný počet obyvateľov v obci
„Bratislava - mestská časť Staré Mesto“ v závislosti od času (index mesiaca).

Označme:

• \(Y_t\): počet obyvateľov v mesiaci \(t\),
• \(t\): časový index (1 = január 2024, 2 = február 2024, …, 12 = december 2024).

Naša pracovná hypotéza:

• predpokladáme, že počet obyvateľov sa v priebehu roka mierne mení (napr. vplyvom migrácie),
• očakávame hladký vývoj v čase, teda silnú autokoreláciu, t. j. počet obyvateľov v mesiaci \(t\) je veľmi podobný počtu v mesiaci \(t-1\),
• očakávame kladný trendový koeficient (\(\beta_1 > 0\)), ak počet obyvateľov má mierne rastúcu tendenciu.

1.2 Príprava databázy, úprava údajov

Databáza je v tzv. „wide“ tvare (mesiace v stĺpcoch). Pre časovú analýzu jednej obce je výhodnejšie mať údaje v „long“ tvare: každý riadok = jedna obec v jednom mesiaci.

Najskôr načítame údaje, potom vyberieme jednu obec a preklopíme mesačné stĺpce do časového radu.

# načítanie databázy
udaje <- read.csv("Databaza.csv",
                  sep = ";",
                  dec = ",",
                  header = TRUE,
                  check.names = FALSE)

# vyberieme jednu konkrétnu obec (možeš zmeniť na inú)
udajeMesto <- udaje[udaje$Obec == "Bratislava - mestská časť Staré Mesto", ]

# preklopenie z wide na long (mesiace do riadkov)
udajeMesto_long <- udajeMesto %>%
  pivot_longer(
    cols = dplyr::starts_with("2024M"),
    names_to = "RokMesiac",
    values_to = "Pocet_obyvatelov"
  ) %>%
  arrange(RokMesiac) %>% 
  mutate(
    Rok = as.numeric(substr(RokMesiac, 1, 4)),
    Mes = as.numeric(substr(RokMesiac, 6, 7)),
    Time = dplyr::row_number()  # 1 = január, 12 = december
  )

udajeMesto_long

---

2. Lineárna regresia v základnom tvare

Ide o odhad rovnice

\[ \text{Pocet\_obyvatelov}_t = \beta_0 + \beta_1 \text{Time}_t + e_t, \]

kde:

• \(\text{Pocet\_obyvatelov}_t\) – počet obyvateľov vo vybranom mesiaci,
• \(\text{Time}_t\) – index mesiaca (1–12),
• \(e_t\) – náhodná zložka (rezíduum).

library(ggplot2)

# lineárny trendový model
model <- lm(Pocet_obyvatelov ~ Time, data = udajeMesto_long)
summary(model)

Call:
lm(formula = Pocet_obyvatelov ~ Time, data = udajeMesto_long)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.232 -12.214   1.946  12.539  23.684 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 47339.47      10.96 4318.14  < 2e-16 ***
Time           24.44       1.49   16.41 1.47e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.81 on 10 degrees of freedom
Multiple R-squared:  0.9642,    Adjusted R-squared:  0.9606 
F-statistic: 269.2 on 1 and 10 DF,  p-value: 1.472e-08

---

3. Autokorelácia rezíduí

V tejto časti sa pozrieme na ďalší dôležitý predpoklad klasického lineárneho regresného modelu – nezávislosť rezíduí. V časových radoch (ako je náš mesačný počet obyvateľov) sa často stáva, že chyba v čase \(t\) je systematicky spätá s chybou v čase \(t-1\), čo nazývame autokoreláciou rezíduí.

3.1 Čo je autokorelácia rezíduí?

Autokorelácia rezíduí je situácia, keď platí:

\[ \operatorname{Corr}(e_t, e_{t-k}) \neq 0 \quad \text{pre niektoré } k \neq 0. \]

Najčastejšie sa skúma autokorelácia prvého rádu:

\[ e_t = \rho e_{t-1} + \nu_t,\quad |\rho| < 1. \]

3.2 Dôsledky autokorelácie rezíduí

• odhady koeficientov \(\hat{\beta}\) sú nestranné,
• ale neefektívne,
• štandardné chyby sú skreslené (často podhodnotené), teda p-hodnoty sa javia menšie,
• t- a F-testy sú potom skreslené.

---

3.3 Detekcia autokorelácie

Grafická informácia

# pridáme do dát fitted values z modelu
udajeMesto_long$fitted <- fitted(model)

# scatterplot + vyrovnaná trendová čiara
ggplot(udajeMesto_long, aes(x = Time, y = Pocet_obyvatelov)) +
  geom_point(color = "steelblue", size = 2) +
  geom_line(aes(y = fitted), color = "red", size = 1) +
  labs(
    title = "Počet obyvateľov: empirické údaje (modrá) vs. odhadnutý trend (červená)",
    x = "Mesiac (2024)",
    y = "Počet obyvateľov"
  ) +
  theme_minimal()

Analýzou obrázka vidíme, že mesačné hodnoty počtu obyvateľov na seba plynulo nadväzujú – v susedných mesiacoch sú si veľmi podobné. Ak by sme si zobrazili rezíduá (rozdiel medzi empirickými a vyrovnanými hodnotami), je prirodzené očakávať, že počas niekoľkých po sebe idúcich mesiacov budú mať podobné znamienko. To je typický vizuálny signál možnej autokorelácie.

# uložíme si reziduá z pôvodného modelu
res <- residuals(model)

---

ACF graf (Autocorrelation Function)

Táto funkcia priraďuje odhad korelácie, ktorá je medzi jednotlivými rezíduami v aktuálnom období a období posunutom (Lag) o \(k\) období späť.

acf(res, lag.max = 4, main = "Autokorelačná funkcia rezíduí")

Na tomto grafe je modrou prerušovanou čiarou vyjadrený aj 95 % interval spoľahlivosti pre hodnotu autokorelačného koeficientu s príslušným posunom. Ak sa niektorý stĺpec nachádza mimo tohto pásma, naznačuje to štatisticky významnú autokoreláciu pri danom posune.

---

Durbin–Watsonov test

Durbin–Watsonov koeficient (DW) je vypočítaný z rezíduí podľa vzorca

\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^{2}}{\sum_{t=1}^{n} e_t^{2}} \]

kde \(n\) je počet pozorovaní. Medzi koeficientom autokorelácie dvoch susedných rezíduí a DW platí približný vzťah

\[ \hat{\rho} \approx 1 - \frac{DW}{2}. \]

Hodnoty:

• blízke nule → silná pozitívna autokorelácia,
• blízke 4 → silná negatívna autokorelácia,
• okolo 2 → žiadny výrazný problém s autokoreláciou.

V praxi sa často používa intuitívne pravidlo: ak sa DW nachádza v intervale približne 1,8 až 2,2, problém autokorelácie väčšinou nepovažujeme za vážny.

library(lmtest)
dwtest(model)

    Durbin-Watson test

data:  model
DW = 2.4468, p-value = 0.6766
alternative hypothesis: true autocorrelation is greater than 0

DW test má isté obmedzenia – regresory nesmú byť časovo posunuté a nesmú obsahovať oneskorené pozorovania vysvetľovanej veličiny ako regresory. Tieto obmedzenia neplatia pre Breusch–Godfreyov test.

---

Breusch–Godfreyov test (BG test)

Čo test testuje

BG test je formálnym testom autokorelácie (sériovej korelácie) rezíduí:

\[ u_t = \rho_1 u_{t-1} + \rho_2 u_{t-2} + \cdots + \rho_p u_{t-p} + \varepsilon_t. \]

Testuje, či sú rezíduá korelované v čase až do rádu \(p\), čo je typické pre časové rady.

Na rozdiel od DW testu:

• testuje autokoreláciu s vyšším posunom (lag 1, 2, 3, …),
• funguje aj pri modeloch s oneskorenou vysvetľovanou premennou ako regresorom,
• pracuje s nekonštantnými regresormi (trend, dummy premenné atď.).

---

Hypotézy

Pre test autokorelácie až do rádu \(p\):

• Nulová hypotéza \(H_0\): žiadna sériová korelácia
  \[ \rho_1 = \rho_2 = \cdots = \rho_p = 0 \]

• Alternatívna hypotéza \(H_1\): séria je autokorelovaná
  Aspoň jedna z \(\rho_j\) je odlišná od nuly.

---

Ako funguje BG test

1. Odhadneme pôvodnú regresiu (u nás trendový model) a získame rezíduá \(e_t\).
2. Spustíme pomocnú regresiu: \[ e_t = \alpha_0 + \alpha_1 x_{1t} + \cdots + \alpha_k x_{kt} + \rho_1 e_{t-1} + \cdots + \rho_p e_{t-p} + v_t, \] kde \(x_{jt}\) sú pôvodné regresory.
3. Z tejto regresie vypočítame: \[ \text{BG} = (n - p) R^2_{\text{aux}} \]
4. Pod \(H_0\) platí približne: \[ \text{BG} \sim \chi^2_p. \]

---

Interpretácia

• Veľká hodnota testovej štatistiky BG (malá p-hodnota) → zamietame \(H_0\) → autokorelácia prítomná.
• Malá hodnota BG (veľká p-hodnota) → nezamietneme \(H_0\) → nie je dôkaz autokorelácie.

---

Praktický výpočet v R

bgtest(model, order = 1)

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

data:  model
LM test = 0.70926, df = 1, p-value = 0.3997

V našich údajoch (mesačný počet obyvateľov vo vybranej obci v roku 2024) sa ukáže, či BG test potvrdí alebo vyvráti prítomnosť autokorelácie prvého rádu. Aj keď pri krátkom časovom rade (len 12 mesiacov) nebýva test príliš silný, pre demonštračné účely použijeme aj postup odstránenia dôsledkov autokorelácie pomocou dynamizácie modelu.

---

Ako riešiť autokoreláciu

Koyckova transformácia a Koyckova rovnica

V predchádzajúcej časti sme sa venovali autokorelácii rezíduí. Teraz rozšírime model o dynamickú štruktúru založenú na tzv. geometricky klesajúcom rozložení oneskorených efektov. Takýto model sa nazýva Koyckov model a jeho ústredným prvkom je Koyckova transformácia.

---

Východiskový model s distribuovaným oneskorením

Uvažujme model:

\[ y_t = \alpha + \beta_0 x_t + \beta_1 x_{t-1} + \beta_2 x_{t-2} + \cdots + u_t. \]

---

Koyckova štruktúra koeficientov

Koyck navrhol, že oneskorené efekty majú geometricky klesajúcu podobu:

\[ \beta_k = \lambda^k \beta_0, \qquad 0 < \lambda < 1. \]

Tým sa distribučné oneskorenie zjednoduší tak, že namiesto nekonečného počtu parametrov odhadujeme len:

• \(\beta_0\) – okamžitý efekt,
• \(\lambda\) – rýchlosť tlmenia efektov oneskorenia.

---

Koyckova transformácia

Pôvodný model:

\[ y_t = \alpha + \beta_0 x_t + \lambda \beta_0 x_{t-1} + \lambda^2 \beta_0 x_{t-2} + \cdots + u_t. \]

Vynásobíme model konštantou \(\lambda\) a posunieme časový index o 1 dozadu:

\[ \lambda y_{t-1} = \lambda\alpha + \beta_0 x_{t-1} + \lambda \beta_0 x_{t-2} + \cdots + \lambda u_{t-1}. \]

Odčítaním oboch výrazov dostaneme:

\[ y_t - \lambda y_{t-1} = \alpha(1-\lambda) + \beta_0 (x_t - \lambda x_{t-1}) + (u_t - \lambda u_{t-1}). \]

To je Koyckova transformácia.

---

Koyckova rovnica (autoregresívny model)

Po algebraickej úprave:

\[ y_t = \alpha(1-\lambda) + \beta_0 x_t + \lambda y_{t-1} + v_t, \]

kde:

\[ v_t = u_t - \lambda u_{t-1}. \]

Toto je Koyckova rovnica – dynamický autoregresívny model so závislosťou na minulosti.

V našom kontexte:

• \(y_t\) = počet obyvateľov v mesiaci \(t\),
• \(x_t\) – napríklad časový trend (Time),
• \(y_{t-1}\) – počet obyvateľov v predchádzajúcom mesiaci.

---

Odstraňovanie problému autokorelácie rezíduí

Odhad Koyckovho modelu v R

Najjednoduchšia implementácia pre naše údaje:

udajeMesto_long <- udajeMesto_long %>%
  arrange(Time) %>%
  mutate(
    Pocet_obyvatelov_lag1 = dplyr::lag(Pocet_obyvatelov)
  )

model_koyck <- lm(Pocet_obyvatelov ~ Time + Pocet_obyvatelov_lag1, 
                  data = udajeMesto_long)

summary(model_koyck)

Call:
lm(formula = Pocet_obyvatelov ~ Time + Pocet_obyvatelov_lag1, 
    data = udajeMesto_long)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.591 -16.219   2.003  13.554  23.918 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)   
(Intercept)           58792.7569 15755.3637   3.732  0.00577 **
Time                     30.9479     8.3103   3.724  0.00584 **
Pocet_obyvatelov_lag1    -0.2422     0.3330  -0.727  0.48779   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.75 on 8 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.9609,    Adjusted R-squared:  0.9512 
F-statistic:  98.4 on 2 and 8 DF,  p-value: 2.329e-06

dwtest(model_koyck)

    Durbin-Watson test

data:  model_koyck
DW = 2.2885, p-value = 0.5587
alternative hypothesis: true autocorrelation is greater than 0

V dynamizovanom modeli vystupuje aj oneskorená hodnota vysvetľovanej premennej (Pocet_obyvatelov_lag1). Ak je jej koeficient kladný a menší ako 1, interpretujeme to tak, že časť úrovne z minulého mesiaca sa prenáša do aktuálneho mesiaca (zotrvačnosť). Pri hodnotení kvality modelu je vhodné porovnať napríklad Adjusted R-squared pôvodného a dynamického modelu, prípadne znova overiť autokoreláciu rezíduí (DW, BG test).

---

Newey–West robustné štandardné chyby

Ešte jednou možnosťou, ako reagovať na prítomnosť autokorelácie (a heteroskedasticity), je použiť robustné štandardné chyby podľa Newey–Westa. Tie nemenia samotné odhady koeficientov, ale korigujú odhady smerodajných odchýlok.

library(sandwich)
library(lmtest)

coeftest(model, vcov = NeweyWest(model))

t test of coefficients:

              Estimate Std. Error  t value  Pr(>|t|)    
(Intercept) 4.7339e+04 6.5101e+00 7271.681 < 2.2e-16 ***
Time        2.4441e+01 7.7275e-01   31.628 2.347e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

4. Záver

Na jednoduchom príklade mesačného počtu obyvateľov vo vybranej slovenskej obci v roku 2024 sme ukázali:

• ako odhadnúť trendový model v čase,
• ako graficky a štatisticky testovať autokoreláciu rezíduí (ACF, Durbin–Watson, BG test),
• ako možno model dynamizovať pomocou Koyckovej rovnice,
• a ako využiť Newey–Westove robustné štandardné chyby.

Autokorelácia rezíduí a z nej vyplývajúca potreba dynamizácie modelov má v ekonometrii veľký význam. Tu sme uviedli len základné prístupy – k zložitejším patrí napríklad Almonov model distribuovaného oneskorenia, metóda Cochran–Orcutt alebo všeobecnejšie ARIMA modely, ktoré sa používajú pri dlhších časových radoch.