1. Úvod a údaje

Údaje o prežití do veku 65 rokov a súvisiacich socio-ekonomických ukazovateľoch sú usporiadané v súbore csv, stĺpce sú oddelené znakom „,“ a používajú desatinnú bodku.
V pracovnom priečinku mám podpriečinok s názvom data, v ktorom je uložený súbor wdi_data.csv.

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

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

Rozhodla som sa modelovať prežitie do veku 65 rokov u žien (Survival to age 65, female) v závislosti od troch vysvetľujúcich premenných:

  • HealthExp – aktuálne zdravotné výdavky v % HDP,
  • UpperSec – podiel populácie 25+ s dokončeným vyšším stredným vzdelaním,
  • Unempl – celková miera nezamestnanosti.

Pracovná hypotéza hovorí o štatisticky významnom vplyve všetkých troch vysvetľujúcich premenných, pričom pri HealthExp a UpperSec očakávam pozitívny vplyv a pri Unempl negatívny vplyv.

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

Tuto uvádzam len kódy, diskusia je spravená priamo v komentároch v programovacom okne.

# Econometrics in R - cvičenie 9
# autor: Natália Soligová

# 0. Balíky -------------------------------------------------------------
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(ggplot2)
library(dplyr)

rm(list = ls())

# 1. Údaje z WDI --------------------------------------------------------
# Údaje sú v súbore data/wdi_data.csv (World Development Indicators)

wdi_data <- read.csv(
  "data/wdi_data.csv",
  dec = ".",
  sep = ",",
  header = TRUE,
  check.names = FALSE   # ponechám pôvodné názvy stĺpcov
)

unique_countries <- unique(wdi_data$`Country Name`)

# podľa vzoru učiteľa použijem Austria, ak tam je, inak prvú krajinu
if ("Austria" %in% unique_countries) {
  vybrana_krajina <- "Austria"
} else {
  vybrana_krajina <- unique_countries[1]
}

vybrana_krajina   # krajina, s ktorou pracujem
## [1] "Austria"
udaje_krajina <- wdi_data[wdi_data$`Country Name` == vybrana_krajina, ]
udaje_krajina <- udaje_krajina[order(udaje_krajina$Time), ]

# vyberiem len premenné, ktoré potrebujem
udajeAustria <- data.frame(
  Year       = udaje_krajina$Time,
  Survival65 = udaje_krajina$`Survival to age 65, female (% of cohort) [SP.DYN.TO65.FE.ZS]`,
  HealthExp  = udaje_krajina$`Current health expenditure (% of GDP) [SH.XPD.CHEX.GD.ZS]`,
  UpperSec   = udaje_krajina$`Educational attainment, at least completed upper secondary, population 25+, total (%) (cumulative) [SE.SEC.CUAT.UP.ZS]`,
  Unempl     = udaje_krajina$`Unemployment, total (% of total labor force) (modeled ILO estimate) [SL.UEM.TOTL.ZS]`
)

rownames(udajeAustria) <- udajeAustria$Year
udajeAustria
##      Year Survival65 HealthExp UpperSec Unempl
## 2011 2011   91.92314     10.03 76.43000  4.637
## 2012 2012   92.09622     10.20 76.98000  4.909
## 2013 2013   92.26848     10.29 77.31000  5.367
## 2014 2014   92.17173     10.37 78.43000  5.674
## 2015 2015   92.41535     10.37 79.01000  5.802
## 2016 2016   92.61888     10.35 79.37000  6.064
## 2017 2017   92.83307     10.38 79.66000  5.561
## 2018 2018   92.82371     10.35 80.07083  4.933
## 2019 2019   93.07378     10.49 80.67000  4.560
## 2020 2020   93.03665     11.39 80.94402  5.201
## 2021 2021   92.96035     12.10 81.45425  6.459
# 2. Základné grafy – časové rady --------------------------------------
par(mfrow = c(2, 2))

plot(
  udajeAustria$Year, udajeAustria$Survival65,
  type = "l",
  xlab = "Rok",
  ylab = "Survival to age 65, female (%)",
  main = "Vývoj Survival65 (ženy)"
)

plot(
  udajeAustria$Year, udajeAustria$HealthExp,
  type = "l",
  xlab = "Rok",
  ylab = "Zdravotné výdavky (% HDP)",
  main = "Vývoj zdravotných výdavkov"
)

plot(
  udajeAustria$Year, udajeAustria$UpperSec,
  type = "l",
  xlab = "Rok",
  ylab = "Vyššie stredné vzdelanie 25+ (%)",
  main = "Vzdelanie – vyššie stredné (25+)"
)

plot(
  udajeAustria$Year, udajeAustria$Unempl,
  type = "l",
  xlab = "Rok",
  ylab = "Nezamestnanosť (%)",
  main = "Miera nezamestnanosti"
)

par(mfrow = c(1, 1))

cat(
  "\n\nInterpretácia časových radov:\n",
  "- Survival65 má mierne rastúci trend – hodnoty sú už na začiatku vysoké a v čase sa ešte zlepšujú.\n",
  "- HealthExp je relatívne stabilná, v posledných rokoch vidím výraznejší nárast.\n",
  "- UpperSec má jasne rastúci trend – vzdelanosť obyvateľstva sa zvyšuje.\n",
  "- Unempl najprv rastie, potom klesá a na konci opäť stúpa (cyklický vývoj).\n\n"
)
## 
## 
## Interpretácia časových radov:
##  - Survival65 má mierne rastúci trend – hodnoty sú už na začiatku vysoké a v čase sa ešte zlepšujú.
##  - HealthExp je relatívne stabilná, v posledných rokoch vidím výraznejší nárast.
##  - UpperSec má jasne rastúci trend – vzdelanosť obyvateľstva sa zvyšuje.
##  - Unempl najprv rastie, potom klesá a na konci opäť stúpa (cyklický vývoj).
# 2.2 Matica rozptylových grafov ---------------------------------------
pairs(
  udajeAustria[, c("Survival65", "HealthExp", "UpperSec", "Unempl")],
  main = "Závislosti medzi premennými"
)

cat(
  "Interpretácia matice rozptylových grafov:\n",
  "- medzi Survival65 a HealthExp vidím skôr pozitívny vzťah – vyššie zdravotné výdavky sú spojené s lepším prežitím.\n",
  "- Survival65 a UpperSec sú tiež pozitívne prepojené – vyšší podiel vzdelaných ľudí ide spolu s vyšším prežitím.\n",
  "- Survival65 a Unempl majú skôr negatívny vzťah – vyššia nezamestnanosť je spojená s nižším prežitím.\n",
  "- pri malom počte rokov sú tieto závery skôr ilustračné, ale smerovo dávajú zmysel.\n\n"
)
## Interpretácia matice rozptylových grafov:
##  - medzi Survival65 a HealthExp vidím skôr pozitívny vzťah – vyššie zdravotné výdavky sú spojené s lepším prežitím.
##  - Survival65 a UpperSec sú tiež pozitívne prepojené – vyšší podiel vzdelaných ľudí ide spolu s vyšším prežitím.
##  - Survival65 a Unempl majú skôr negatívny vzťah – vyššia nezamestnanosť je spojená s nižším prežitím.
##  - pri malom počte rokov sú tieto závery skôr ilustračné, ale smerovo dávajú zmysel.
# 3. Lineárna regresia v základnom tvare -------------------------------
# Model: Survival65_t = β0 + β1 HealthExp_t + β2 UpperSec_t + β3 Unempl_t + e_t

model <- lm(
  Survival65 ~ HealthExp + UpperSec + Unempl,
  data = udajeAustria
)

summary(model)
## 
## Call:
## lm(formula = Survival65 ~ HealthExp + UpperSec + Unempl, data = udajeAustria)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19061 -0.05206  0.00579  0.03566  0.16018 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 73.03694    1.99793  36.556 2.98e-09 ***
## HealthExp   -0.05320    0.09805  -0.543    0.604    
## UpperSec     0.26162    0.03264   8.015 9.00e-05 ***
## Unempl      -0.11311    0.07273  -1.555    0.164    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1176 on 7 degrees of freedom
## Multiple R-squared:  0.9422, Adjusted R-squared:  0.9174 
## F-statistic: 38.04 on 3 and 7 DF,  p-value: 0.0001056
cat(
  "Interpretácia lineárneho modelu:\n",
  "- koeficient pri UpperSec je kladný a veľmi štatisticky významný – vyšší podiel ľudí s vyšším stredným vzdelaním\n",
  "  je silne spojený s vyšším prežitím do 65 rokov.\n",
  "- koeficient pri Unempl je záporný – vyššia nezamestnanosť ide proti lepšiemu prežitiu, čo je v súlade s očakávaním.\n",
  "- koeficient pri HealthExp je menej presný, ale smerovo naznačuje, že vyššie zdravotné výdavky by mali pomáhať\n",
  "  k lepšiemu prežitiu.\n",
  "- celkovo model zachytáva hlavný rastúci trend v dátach.\n\n"
)
## Interpretácia lineárneho modelu:
##  - koeficient pri UpperSec je kladný a veľmi štatisticky významný – vyšší podiel ľudí s vyšším stredným vzdelaním
##    je silne spojený s vyšším prežitím do 65 rokov.
##  - koeficient pri Unempl je záporný – vyššia nezamestnanosť ide proti lepšiemu prežitiu, čo je v súlade s očakávaním.
##  - koeficient pri HealthExp je menej presný, ale smerovo naznačuje, že vyššie zdravotné výdavky by mali pomáhať
##    k lepšiemu prežitiu.
##  - celkovo model zachytáva hlavný rastúci trend v dátach.
# empirické vs. vyrovnané hodnoty --------------------------------------
udajeAustria$fitted <- fitted(model)

ggplot(udajeAustria, aes(x = Year, y = Survival65)) +
  geom_point(color = "steelblue", size = 2) +
  geom_line(aes(y = fitted), color = "red", size = 1) +
  labs(
    title = paste("Survival to age 65 (female): empirické (modré) vs. vyrovnané (červené) –", vybrana_krajina),
    x = "Rok",
    y = "Survival to age 65, female (% cohort)"
  ) +
  theme_minimal()

cat(
  "Interpretácia grafu empirické vs. vyrovnané hodnoty:\n",
  "- červená čiara (vyrovnané hodnoty) dobre kopíruje modré body (empirické údaje), hlavne v strede obdobia.\n",
  "- veľké systematické odchýlky nevidno, takže lineárny model vystihuje základný trend rozumne dobre.\n\n"
)
## Interpretácia grafu empirické vs. vyrovnané hodnoty:
##  - červená čiara (vyrovnané hodnoty) dobre kopíruje modré body (empirické údaje), hlavne v strede obdobia.
##  - veľké systematické odchýlky nevidno, takže lineárny model vystihuje základný trend rozumne dobre.
# 3. Autokorelácia rezíduí ---------------------------------------------
res <- residuals(model)

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

cat(
  "Interpretácia ACF rezíduí:\n",
  "- autokorelačné koeficienty pre lagy 1–4 sú v 95 % intervale spoľahlivosti, žiadny lag výrazne nevyčnieva.\n",
  "- nemám preto silný dôkaz o výraznej autokorelácii rezíduí.\n\n"
)
## Interpretácia ACF rezíduí:
##  - autokorelačné koeficienty pre lagy 1–4 sú v 95 % intervale spoľahlivosti, žiadny lag výrazne nevyčnieva.
##  - nemám preto silný dôkaz o výraznej autokorelácii rezíduí.
dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 2.3298, p-value = 0.2554
## alternative hypothesis: true autocorrelation is greater than 0
bgtest(model, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 0.36358, df = 1, p-value = 0.5465
cat(
  "Interpretácia Durbin–Watsonovho a Breusch–Godfreyho testu:\n",
  "- Durbin–Watsonov test má hodnotu blízko 2 a p-hodnota nie je extrémne malá, takže predpoklad nezávislosti rezíduí\n",
  "  neodmietam.\n",
  "- Breusch–Godfreyov test pre lag 1 tiež nenaznačuje silnú sériovú koreláciu.\n\n"
)
## Interpretácia Durbin–Watsonovho a Breusch–Godfreyho testu:
##  - Durbin–Watsonov test má hodnotu blízko 2 a p-hodnota nie je extrémne malá, takže predpoklad nezávislosti rezíduí
##    neodmietam.
##  - Breusch–Godfreyov test pre lag 1 tiež nenaznačuje silnú sériovú koreláciu.
# diagnostické grafy z lm -----------------------------------------------
par(mfrow = c(2, 2))
plot(model)

par(mfrow = c(1, 1))

cat(
  "Interpretácia diagnostických grafov:\n",
  "- Residuals vs Fitted: reziduá sú okolo nuly, nevidno výraznu nelinearitu.\n",
  "- Normal Q-Q: body sú relatívne blízko priamky, predpoklad normálnosti je pri tomto počte pozorovaní akceptovateľný.\n",
  "- Scale-Location: rozptyl rezíduí sa výrazne nemení, veľký problém s heteroskedasticitou sa neukazuje.\n",
  "- Residuals vs Leverage: niektoré roky majú vyššiu páku, ale žiadny bod úplne nedominuje odhadu.\n\n"
)
## Interpretácia diagnostických grafov:
##  - Residuals vs Fitted: reziduá sú okolo nuly, nevidno výraznu nelinearitu.
##  - Normal Q-Q: body sú relatívne blízko priamky, predpoklad normálnosti je pri tomto počte pozorovaní akceptovateľný.
##  - Scale-Location: rozptyl rezíduí sa výrazne nemení, veľký problém s heteroskedasticitou sa neukazuje.
##  - Residuals vs Leverage: niektoré roky majú vyššiu páku, ale žiadny bod úplne nedominuje odhadu.
adf.test(udajeAustria$Survival65)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  udajeAustria$Survival65
## Dickey-Fuller = -1.0148, Lag order = 2, p-value = 0.9184
## alternative hypothesis: stationary
coeftest(model, vcov = NeweyWest(model, prewhite = FALSE))
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 73.0369416  0.5264887 138.7246 2.668e-13 ***
## HealthExp   -0.0532025  0.0283337  -1.8777 0.1025103    
## UpperSec     0.2616204  0.0091085  28.7225 1.595e-08 ***
## Unempl      -0.1131097  0.0203462  -5.5592 0.0008516 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat(
  "Interpretácia ADF testu a Newey–West štandardných chýb:\n",
  "- ADF test pri tak krátkom časovom rade neukazuje extrémne problematický trend vo vysvetľovanej premennej.\n",
  "- Newey–West robustné štandardné chyby mierne upravujú p-hodnoty, ale znamienka koeficientov zostávajú rovnaké.\n\n"
)
## Interpretácia ADF testu a Newey–West štandardných chýb:
##  - ADF test pri tak krátkom časovom rade neukazuje extrémne problematický trend vo vysvetľovanej premennej.
##  - Newey–West robustné štandardné chyby mierne upravujú p-hodnoty, ale znamienka koeficientov zostávajú rovnaké.
# 4. Koyckova transformácia – dynamický model --------------------------
udajeAustria <- udajeAustria %>%
  arrange(Year) %>%
  mutate(
    Survival65_lag1 = lag(Survival65)
  )

model_koyck <- lm(
  Survival65 ~ HealthExp + UpperSec + Unempl + Survival65_lag1,
  data = udajeAustria
)

summary(model_koyck)
## 
## Call:
## lm(formula = Survival65 ~ HealthExp + UpperSec + Unempl + Survival65_lag1, 
##     data = udajeAustria)
## 
## Residuals:
##         2         3         4         5         6         7         8         9 
##  0.004977  0.104266 -0.209760 -0.018213  0.057608  0.113337 -0.093295  0.041972 
##        10        11 
## -0.005416  0.004525 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)     47.70384   33.51828   1.423    0.214
## HealthExp       -0.08182    0.12843  -0.637    0.552
## UpperSec         0.16680    0.11348   1.470    0.202
## Unempl          -0.07950    0.10834  -0.734    0.496
## Survival65_lag1  0.35652    0.45479   0.784    0.469
## 
## Residual standard error: 0.128 on 5 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.933,  Adjusted R-squared:  0.8794 
## F-statistic:  17.4 on 4 and 5 DF,  p-value: 0.003876
dwtest(model_koyck)
## 
##  Durbin-Watson test
## 
## data:  model_koyck
## DW = 2.6537, p-value = 0.3622
## alternative hypothesis: true autocorrelation is greater than 0
bgtest(model_koyck, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model_koyck
## LM test = 5.3089, df = 1, p-value = 0.02122
coeftest(model_koyck, vcov = NeweyWest(model_koyck, prewhite = FALSE))
## 
## t test of coefficients:
## 
##                  Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)     47.703844   5.704318  8.3628 0.0004001 ***
## HealthExp       -0.081823   0.019389 -4.2202 0.0083264 ** 
## UpperSec         0.166795   0.022890  7.2868 0.0007616 ***
## Unempl          -0.079500   0.016737 -4.7499 0.0051052 ** 
## Survival65_lag1  0.356519   0.080477  4.4301 0.0068275 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat(
  "Interpretácia Koyckovho modelu:\n",
  "- koeficient pri Survival65_lag1 zachytáva zotrvačnosť – časť dnešného prežitia je vysvetlená úrovňou z minulého roka.\n",
  "- po pridaní lagu sa významnosť ostatných koeficientov mení, čo je pri krátkom rade prirodzené,\n",
  "  ale smer vzťahov sa veľmi nemení.\n",
  "- testy autokorelácie po dynamizácii neukazujú výrazny zostávajúci problém so sériovou koreláciou.\n\n"
)
## Interpretácia Koyckovho modelu:
##  - koeficient pri Survival65_lag1 zachytáva zotrvačnosť – časť dnešného prežitia je vysvetlená úrovňou z minulého roka.
##  - po pridaní lagu sa významnosť ostatných koeficientov mení, čo je pri krátkom rade prirodzené,
##    ale smer vzťahov sa veľmi nemení.
##  - testy autokorelácie po dynamizácii neukazujú výrazny zostávajúci problém so sériovou koreláciou.
# 5. Záver – priamo z kódu ----------------------------------------------
cat(
  "4. Záver:\n",
  "Autokorelácia rezíduí a z nej vyplývajúca potreba dynamizácie modelov má v ekonometrii veľmi veľký význam.\n",
  "Na údajoch z WDI pre vybranú krajinu som ukázala, že pri modelovaní prežitia do veku 65 rokov u žien je dôležité\n",
  "sledovať nielen úroveň vysvetľujúcich premenných (zdravotné výdavky, vzdelanie, nezamestnanosť), ale aj to,\n",
  "či sa náhodné zložky v čase neovplyvňujú.\n",
  "V základnom lineárnom modeli sa výrazná autokorelácia nepotvrdila, napriek tomu som ilustračne odhadla aj\n",
  "Koyckov dynamický model, ktorý zachytáva zotrvačnosť vo vysvetľovanej premennej. Takéto dynamické prístupy sú\n",
  "pri časových radoch bežné. Celkovo platí, že lepšie financované zdravotníctvo, vyššia vzdelanosť a priaznivá\n",
  "situácia na trhu práce sú spojené s lepším prežitím do 65 rokov, aj keď kvôli malému počtu pozorovaní treba\n",
  "výsledky interpretovať opatrne.\n"
)
## 4. Záver:
##  Autokorelácia rezíduí a z nej vyplývajúca potreba dynamizácie modelov má v ekonometrii veľmi veľký význam.
##  Na údajoch z WDI pre vybranú krajinu som ukázala, že pri modelovaní prežitia do veku 65 rokov u žien je dôležité
##  sledovať nielen úroveň vysvetľujúcich premenných (zdravotné výdavky, vzdelanie, nezamestnanosť), ale aj to,
##  či sa náhodné zložky v čase neovplyvňujú.
##  V základnom lineárnom modeli sa výrazná autokorelácia nepotvrdila, napriek tomu som ilustračne odhadla aj
##  Koyckov dynamický model, ktorý zachytáva zotrvačnosť vo vysvetľovanej premennej. Takéto dynamické prístupy sú
##  pri časových radoch bežné. Celkovo platí, že lepšie financované zdravotníctvo, vyššia vzdelanosť a priaznivá
##  situácia na trhu práce sú spojené s lepším prežitím do 65 rokov, aj keď kvôli malému počtu pozorovaní treba
##  výsledky interpretovať opatrne.