Econometrics in R - cvičenie 9

S využitím databázy WHO Life Expactancy Data database.

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

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

1. Úvod a údaje

Údaje o očakávanej dĺžke života sú usporiadané v súbore csv, stĺpce sú oddelené znakom “,” a používajú desatinnú čiarku. V pracovnom priečinku som vytvoril podpriečinok s názvom udaje, aby som oddelil údaje od zvyšku projektu. Môj priečinok sa tiež nazýva udaje.

Nie všetky údaje budú použité, preto som vybral len niektoré stĺpce pre neskoršie použitie.

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

Rozhodol som sa modelovať strednú dĺžku života Life.expectancy v závislosti od troch vysvetľujúcich premenných a to BMI, HDP na obyvateľa GDP a stredný počet rokov štúdia Schooling.

Naša pracovná hypotéza hovorí o štatisticky významnom vplyve všetkých troch vysvetľujúcich premenných, pričom u premenných GDP a Schooling by malo ísť o pozitívny vplyv (očakávame kladné znamienko odhadovaného regresného koeficienta) a v prípade BMI by malo ísť of negatívny vplyv (so záporným znamienkom)

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

Tuto uvádzam len kódy, diskusia je urobená na predchádzajúcich cvičeniach.

udaje <- read.csv("udaje/Life_Expectancy_Data.csv",dec=".",sep=",",header = TRUE)
# select just the record from 2015
udajeAustria <- udaje[udaje$Country=="Austria",c("Life.expectancy","BMI","GDP","Schooling","Year")]
rownames(udajeAustria)<-udajeAustria$Year
udajeAustria <- udajeAustria[order(udajeAustria$Year), ]


udajeAustria

NA

2. Lineárna regresia v základnom tvare

Ide o odhad rovnice

\[Life.expectancy_t = \beta_0 + \beta_1 BMI_t + \beta_2 GDP_t + \beta_3Schooling_t + e_t\]

library(ggplot2)

# your regression model
model <- lm(Life.expectancy ~ BMI + GDP + Schooling, data = udajeAustria)
summary(model)

Call:
lm(formula = Life.expectancy ~ BMI + GDP + Schooling, data = udajeAustria)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6260 -1.0164 -0.3103  0.5488  4.7682 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) 51.2975638 25.8395778   1.985   0.0705 .
BMI         -0.0709884  0.7207499  -0.098   0.9232  
GDP          0.0002143  0.0001517   1.413   0.1831  
Schooling    1.6279409  1.9582731   0.831   0.4220  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.425 on 12 degrees of freedom
Multiple R-squared:  0.5132,    Adjusted R-squared:  0.3915 
F-statistic: 4.217 on 3 and 12 DF,  p-value: 0.02976

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 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ú podhodnotené teda p hodnoty sa javia menšie, než by mali byť - dochádza k falošnému zdaniu o štatistickej významnosti.
• teda t- a F-testy sú skreslené.

---

3.3 Detekcia autokorelácie

Grafická informácia

library(ggplot2)

# add fitted values to the dataframe
udajeAustria$fitted <- fitted(model)

# scatterplot + regression line + spline smoother
ggplot(udajeAustria, aes(x = Year, y = Life.expectancy)) +
  geom_point(color = "steelblue", size = 2) +
  
  # regression fitted line
  geom_line(aes(y = fitted), color = "red", size = 1) +
  

  labs(
    title = "Life Expectancy in Austria: Empirical Data (blue) vs. Fitted Data (Red)",
    x = "Year",
    y = "Life Expectancy"
  ) +
  theme_minimal()

Analýzou obrázka vidíme že empirické hodnoty sa počas niekoľkých susedných kompaktných období stále nachádzajú pod vyrovnanou hodnotou ako najmä v rokoch 2005 až 2007 respektíve 2013 až 2015. Vždy, ak sa dajú identifikovaťkompaktné úseky v čase, keď rezíduá majú rovnaké znamienko, signalizuje to možnú prítomnosť autokorelácie náhodných zložiek.

# ulosime si rezidua z povodneho modelu - nazvaneho model
res <- residuals(model)

---

ACF graf (Autocorrelation Function)

Táto funkcia priradzuje 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 vujadrený aj 95 % interval spoľahlivosti pre hodnotu autokorelačného koeficientu s príslušným posunom. Podľa výsledkov sa zdá že žiaden koeficient autokorelácie nie je 3tatistický významný, keďže pre posun väčší alebo rovný 1 je odhadnutý autokorelačný koeficient stále v hraniciach konfidenčného intervalu (intervalu spoľahlivosti).

---

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} \]

Durbin Watsonov koeficient nie je v pôvodnej forme tabelovaný. My ho ale v našich výstupoch tabelovaný máme a poskytuje teda aj p-hodnotu. Tento koeficient nadobúda hodnoty od nula po 4. Hodnoty blízke nule signalizujú významnú pozitívnu autokoreláciu ktorá je typická v mnohých časových radoch a hodnoty nad 2 vyjadrujú negatívnu autokoreláciu. V praxi poväčšine používame intuitívne pravidlo a to, či Durbin-Watsonov koeficient sa nachádza niekde v intervale 1,8 až 2,2. Problém autokoreácie potom nemusíme riešiť. V našom prípade je ten koeficient veľmi nízky a teda nám to signalizuje prítomnosť významnej autokorelácie pozitívnej autokorelácie rezíduí prvého rádu.

library(lmtest)
dwtest(model)

    Durbin-Watson test

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

DW test má svoje obmedzenie použitia - regresory nesmú byť časovo posunuté a nesmú obsahovať časovo oneskorené pozorovania vysvetľovanej veličiny ako regresory. Tieto ohraničenia neplatia pre Breush-Godfreyov test.

---

Breusch–Godfreyov test (BG test)

Čo test testuje

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

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

Testuje, či rezíduá sú koelované v čase a to aj pre \(p\) posunov (lagov), čo je typické pre časové rady (nie pre prierezové údaje!).

Na rozdiel od DW testu, BG test testuje:

• autokorelácie s vyšším posunom (lag 1, 2, 3, …),
• pracuje aj s modelmi, ktoré obsahujú časovo posunuté vysvetľované premenné ako regresory. (BG test),
• pracuje s nekonštantnými regresormi.

---

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ériová korelácia prítomná

Aspoň jedna z \(\rho_j\) je odlišná od nuly.

---

Ako funguje Breuschov–Godfreyho test

1. Odhadnite regresiu pomocou OLS a získajte rezíduá \(e_t\).

2. Spustite 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 \]

Regresory z pôvodného modelu musia byť zahrnuté.

3. Vypočítajte testovú štatistiku:

\[ \text{BG} = (n - p) R^2_{\text{aux}} \]

4. Pod nulovou hypotézou:

\[ \text{BG} \sim \chi^2_p \]

---

Interpretácia

• Veľká štatistika BG (malá hodnota p) → zamietnuť \(H_0\) → zistená autokorelácia.
  
• Malá štatistika BG (veľká hodnota p) → nezamietnuť \(H_0\) → žiadny dôkaz autokorelácie.

---

Intuícia

Ak sú rezíduá korelované v čase, oneskorené rezíduá \(e_{t-1}, e_{t-2}, \ldots\) pomôžu vysvetliť variabilitu v \(e_t\).

---

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 = 1.674, df = 1, p-value = 0.1957

V našich údajoch sme identifikovali týmto testom nízku úroveň autokorelácie s posunom o 1, teda závery DW testu a BG testu sú rozporuplné. Preto aspoň z demonštratívneho dôvodu pristúpime ku odstraňovaniu dôsledkov autokorelácie.

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. \]

Takýto model obsahuje veľký počet parametrov \(\beta_k\). Problémy:

• vysoká multikolinearita medzi \(x_t, x_{t-1}, x_{t-2}, \dots\)
• veľká variancia odhadov
• nestabilita parametrov pri malých výberoch

---

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

Začnime pôvodným modelom:

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

Vynásobme celý model konštantou \(\lambda\) a posuňme časový indexdozadu o 1:

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

Teraz odčítame oba výrazy:

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

Toto je Koyckova transformácia.

---

Koyckova rovnica (autoregresívny model)

Po algebraickej úprave získame:

\[ 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.

Interpretácia:

• \(\lambda\) – rýchlosť prispôsobenia (čím bližšie 1, tým pomalšia reakcia)
• \(\beta_0\) – okamžitý efekt \(x_t\)
• dlhodobý efekt:

\[ \frac{\beta_0}{1-\lambda} \]

---

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

Odhad Koyckovho modelu v R

Najjednoduchšia implementácia:

library(dplyr)

udajeAustria <- udajeAustria %>%
  arrange(Year) %>%
  mutate(
    Life.expectancy_lag1 = lag(Life.expectancy)
  )

model_koyck <- lm(Life.expectancy ~ GDP + BMI + Schooling +  Life.expectancy_lag1, 
                  data = udajeAustria)

summary(model_koyck)

Call:
lm(formula = Life.expectancy ~ GDP + BMI + Schooling + Life.expectancy_lag1, 
    data = udajeAustria)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9183 -0.7982 -0.2828  0.9154  4.1247 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)          35.6572901 28.7114197   1.242    0.243
GDP                   0.0001254  0.0001675   0.749    0.471
BMI                  -0.1107985  0.7403988  -0.150    0.884
Schooling             0.7848417  2.1058277   0.373    0.717
Life.expectancy_lag1  0.4245017  0.3200273   1.326    0.214

Residual standard error: 2.449 on 10 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.548, Adjusted R-squared:  0.3672 
F-statistic: 3.031 on 4 and 10 DF,  p-value: 0.07057

dwtest(model_koyck)

    Durbin-Watson test

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

Dynamizovaný model nám neidentifikuje žiaden regresor ako štatisicky významný. Koeficient pri y_lag1, čo je očakávaná dĺžka dožitia v minulom období, je kladný a menší ako 1. Hodnota 0.42 nám hovorí o zotrvačnom vplyve vysvetľovanej premennej z minulosti (niečo ako 42 percent očakávanej dĺžky života sa prenesie z minulého obdobia - čo je interpretáucia v tomto konkrétnom prípade dosť kostrbatá).

Pri hodnotení kvality nového modelu porovnajme len Adjusted R-squared. Z neho vyplýva, že pôvodný model nám poskytuje lepšie výsledky.

Newey–West robustné štandardné chyby

library(sandwich)
library(lmtest)

coeftest(model, vcov = NeweyWest(model))

t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)   
(Intercept)  5.1298e+01  1.4030e+01  3.6563 0.003288 **
BMI         -7.0988e-02  1.7758e-01 -0.3998 0.696364   
GDP          2.1430e-04  5.2481e-05  4.0834 0.001517 **
Schooling    1.6279e+00  7.9867e-01  2.0383 0.064182 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4. Záver

Autokorelácia rezíduí a z nej vyplývajúca potreba dynamizácie modelov má ekonómometrii veľmi veľký význam. My sme tu uviedli len najzákladnejšie problémy. K trochu zložitejším patrí Almonovej model riešenia dynamizácie alebo metóda Cochran-Orcutt. Samostatnou metodológiou je ale metodika arima modelov ktorá je dnes populárna a pri dlhých časových radoch sa používa najčastejšie.