Práca s databázou - import údajov, grafy, štatistiky

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

#Import udajov

udaje <- read.csv("Employment_Unemployment_GDP_data.csv",
header = TRUE,
sep = ",",
dec = ".",
stringsAsFactors = FALSE)


udajeSlovakRepublic <- udaje[udaje$Country=="Slovak Republic",c("Unemployment.Rate","Employment.Sector..Agriculture", "Year","GDP..in.USD.", "Employment.Sector..Industry","Employment.Sector..Services")]
rownames(udajeSlovakRepublic)<-udajeSlovakRepublic$Year
udajeSlovakRepublic <- udajeSlovakRepublic[order(udajeSlovakRepublic$Year), ]
udajeSlovakRepublic$logGDP <- log(udajeSlovakRepublic$GDP..in.USD.)



udajeSlovakRepublic

NA

2. Lineárna regresia v základnom tvare

V tejto časti sme urobili základný regresný model, kde skúmame, čo ovplyvňuje mieru nezamestnanosti na Slovensku. Do modelu sme dali podiel zamestnanosti v jednotlivých sektoroch a logaritmus HDP, ktorý lepšie zachytáva ekonomický rast.

library(ggplot2)

# môj regresný logaritmovaný model
model_log <- lm(Unemployment.Rate ~ 
                  Employment.Sector..Agriculture +
                  Employment.Sector..Industry +
                  Employment.Sector..Services +
                  logGDP,
                data = udajeSlovakRepublic)

summary(model_log)

Call:
lm(formula = Unemployment.Rate ~ Employment.Sector..Agriculture + 
    Employment.Sector..Industry + Employment.Sector..Services + 
    logGDP, data = udajeSlovakRepublic)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3249 -2.3479  0.4943  1.6275  4.9565 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     2.531e+06  2.243e+06   1.128 0.269106    
Employment.Sector..Agriculture -2.531e+04  2.243e+04  -1.128 0.269112    
Employment.Sector..Industry    -2.530e+04  2.243e+04  -1.128 0.269190    
Employment.Sector..Services    -2.530e+04  2.243e+04  -1.128 0.269175    
logGDP                         -1.427e+01  3.662e+00  -3.896 0.000583 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.895 on 27 degrees of freedom
Multiple R-squared:  0.4891,    Adjusted R-squared:  0.4134 
F-statistic: 6.461 on 4 and 27 DF,  p-value: 0.0008785

Z výsledkov vyšlo, že premenné Agriculture, Industry a Services nie sú štatisticky významné. To znamená, že podľa tohto modelu nevieme potvrdiť, že tieto sektory majú priamy vplyv na nezamestnanosť. Na druhú stranu, premenná logGDP je významná a jej koeficient je záporný, takže keď HDP rastie, nezamestnanosť má tendenciu klesať. F-test ukazuje, že model je ako celok významný, takže dáva zmysel ho ďalej používať.

3. Detekcia autokorelácie

Grafická informácia

library(ggplot2)

udajeSlovakRepublic$fitted <- fitted(model_log)

ggplot(udajeSlovakRepublic, aes(x = Year, y = Unemployment.Rate)) +
  geom_point(color = "steelblue", size = 3) +
  geom_line(aes(y = fitted), color = "red", size = 2) +
  labs(
    title = "Unemployment rate in Slovak Republic: Empirical (blue) vs. Fitted (red)",
    x = "Year",
    y = "Unemployment rate (%)"
  ) +
  theme_minimal()

V grafe porovnávame reálne hodnoty nezamestnanosti s hodnotami, ktoré predpovedal model. Na grafe môžeme vidieť, že červená čiara nie presne lemuje modré body. Toto naznačuje, že model nemusí úplne zachytávať časový vývoj nezamestnanosti.

Zároveň je možné vidieť, že hodnoty idú v čase „vlnovito“, čo môže byť znak, že medzi rokmi existuje nejaká závislosť. Práve preto budeme ďalej testovať, či v dátach nie je autokorelácia.

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

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í")

Z ACF grafu vidíme, že reziduá nášho modelu sú výrazne autokorelované. Hodnota pri lag 1 je nad hranicou významnosti, čo znamená, že reziduá z jedného roka sú silne spojené s reziduami z predchádzajúceho roka. Podobne aj lag 2 a lag 3 sú stále relatívne vysoké.

Toto naznačuje, že model nezachytáva časovú závislosť v dátach a pravdepodobne v ňom vzniká autokorelácia, ktorú musíme ďalej otestovať aj formálnymi štatistickými testami.

Durbin–Watsonov test

Pri Durbin–Watsonovom teste pracujeme s týmito hypotézami:

• H₀ : v reziduách nie je prítomná autokorelácia prvého rádu. Chyby medzi rokmi sú náhodné a nesúvisia spolu.
• H₁ : v reziduách je prítomná autokorelácia prvého rádu. Chyby medzi rokmi sú navzájom závislé.

library(lmtest)
dwtest(model_log
       )

    Durbin-Watson test

data:  model_log
DW = 0.54927, p-value = 9.366e-09
alternative hypothesis: true autocorrelation is greater than 0

Výsledok testu nám v našom prípade ukazuje hodnotu DW = 0.549, čo je výrazne nižšie ako hodnota okolo 2, ktorá by znamenala žiadnu autokoreláciu. P-hodnota je extrémne malá (9.366e-09), takže nulovú hypotézu zamietame. Z toho vyplýva, že v našom modeli je prítomná pozitívna autokorelácia rezíduí.

Breusch–Godfreyov test (BG test)

Durbin–Watsonov test nám ukázal, že v našom modeli je problém s autokoreláciou. Na potvrdenie tohto výsledku použijeme Breusch–Godfreyov test, ktorý je všeobecnejší a vie zachytiť aj vyšší rád autokorelácie.

Pri BG teste pracujeme s týmito hypotézami:

• H₀ : v reziduách nášho modelu nie je autokorelácia daného rádu. Chyby sú náhodné a medzi rokmi nie sú systematicky prepojené.
• H₁ : v reziduách je prítomná autokorelácia daného rádu. Chyby medzi rokmi spolu súvisia a model nezachytáva časový priebeh úplne správne.

bgtest(model_log, order = 1)

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

data:  model_log
LM test = 15.78, df = 1, p-value = 7.114e-05

Test nám vyšiel s p-hodnotou 7.114e-05, čo je oveľa menej ako hranica 0,05. Preto nulovú hypotézu zamietame a berieme to tak, že v našom modeli je prítomná autokorelácia prvého rádu.

Tento výsledok potvrdzuje to, čo sme videli aj v ACF grafe a v Durbin–Watsonovom teste. Chyby modelu sa v čase opakujú a nadväzujú na seba, takže model nezachytáva časový vývoj nezamestnanosti úplne správne. Ide o problém v špecifikácii modelu, nie o problém v samotných dátach.

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

Keďže náš pôvodný model mal problém s autokoreláciou, pokúšame sa ho odstrániť pridaním oneskorenej závislej premennej. Ide o tzv. Koyckov model. Do modelu sme pridali nezamestnanosť z minulého roka (lag 1), keďže predpokladáme, že vývoj nezamestnanosti sa mení postupne a môže byť ovplyvnený hodnotami z predchádzajúceho obdobia.

Po pridaní tejto premennej odhadneme nový model a následne vykonáme Durbin–Watsonov test, aby sme skontrolovali, či autokorelácia zmizla. Ak je hodnota DW bližšie k číslu 2 a p-hodnota je väčšia ako 0,05, znamená to, že autokoreláciu sa podarilo odstrániť.

library(dplyr)

udajeSlovakRepublic <- udajeSlovakRepublic %>%
  arrange(Year) %>%
  mutate(
    Unemployment_lag1 = lag(Unemployment.Rate)
  )

model_koyck <- lm(Unemployment.Rate ~ 
                    Employment.Sector..Agriculture +
                    Employment.Sector..Industry +
                    Employment.Sector..Services +
                    logGDP +
                    Unemployment_lag1,
                  data = udajeSlovakRepublic)

summary(model_koyck)

Call:
lm(formula = Unemployment.Rate ~ Employment.Sector..Agriculture + 
    Employment.Sector..Industry + Employment.Sector..Services + 
    logGDP + Unemployment_lag1, data = udajeSlovakRepublic)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8401 -0.9695 -0.2207  0.9802  2.8456 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     2.047e+06  1.098e+06   1.864   0.0741 .  
Employment.Sector..Agriculture -2.047e+04  1.098e+04  -1.864   0.0741 .  
Employment.Sector..Industry    -2.047e+04  1.098e+04  -1.864   0.0741 .  
Employment.Sector..Services    -2.047e+04  1.098e+04  -1.864   0.0741 .  
logGDP                         -3.746e+00  2.237e+00  -1.674   0.1065    
Unemployment_lag1               8.215e-01  9.572e-02   8.582 6.37e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.414 on 25 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.8872,    Adjusted R-squared:  0.8646 
F-statistic: 39.31 on 5 and 25 DF,  p-value: 4.654e-11

Výsledky Koyckovho modelu

Keď sme do modelu pridali oneskorenú nezamestnanosť z minulého roka, tak sa výsledky zlepšili. Premenná Unemployment_lag1 vyšla veľmi významná, čo znamená, že nezamestnanosť z predchádzajúceho roka silno ovplyvňuje nezamestnanosť v aktuálnom roku. To dáva zmysel, keďže nezamestnanosť sa väčšinou nemení skokovo, ale postupne.

Ostatné premenné už nie sú také dôležité ako predtým, pretože veľkú časť zmeny v nezamestnanosti vysvetlí práve hodnota z minulého roka. Tým, že sme túto premennú pridali, model sa lepšie trafí do toho, ako sa nezamestnanosť správa v čase.

Hodnota R² sa zvýšila na približne 0.89, čo je veľmi dobré. Znamená to, že náš nový model už vie vysvetliť skoro celú zmenu v nezamestnanosti. Celkovo tento model pôsobí stabilnejšie a presnejšie než ten pôvodný.

Durbin–Watsonov test po odhade Koyckovho modelu

Na novom modeli sme spravili opäť Durbin–Watsonov test, aby sme zistili, či sa autokorelácia podarila odstrániť. V našom prípade sme dostali hodnotu DW = 1.3468, čo je síce bližšie k hodnote 2 ako predtým, ale stále to nie je úplne ideálne. P-hodnota je 0.003, takže nulovú hypotézu opäť zamietame.

To znamená, že aj keď sa autokorelácia znížila, stále tam je a model ju úplne neodstránil. Koyckov model síce pomohol zlepšiť výsledky a spravil model presnejší, ale autokorelácia sa nestratila úplne.

dwtest(model_koyck)

    Durbin-Watson test

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

Na novom modeli sme spravili opäť Durbin–Watsonov test, aby sme zistili, či sa autokorelácia podarila odstrániť. V našom prípade sme dostali hodnotu DW = 1.3468, čo je síce bližšie k hodnote 2 ako predtým, ale stále to nie je úplne ideálne. P-hodnota je 0.003, takže nulovú hypotézu opäť zamietame.

Koyckov model síce pomohol zlepšiť výsledky a spravil model presnejší, ale autokorelácia sa nestratila úplne.

Newey–West robustné štandardné chyby

Keďže náš model má problémy s autokoreláciou a niektoré predpoklady klasickej regresie nie sú splnené, použili sme tzv. robustné smerodajné chyby. Tie nemenia samotné koeficienty, ale upravujú ich smerodajné chyby tak, aby boli výsledky spoľahlivejšie aj vtedy, keď sa model nespráva úplne ideálne.

coeftest(model_koyck, vcov = vcovHC(model_koyck, type = "HC3"))

t test of coefficients:

                                  Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)                     2.0472e+06  9.4600e+05  2.1640   0.04022 *  
Employment.Sector..Agriculture -2.0471e+04  9.4859e+03 -2.1580   0.04073 *  
Employment.Sector..Industry    -2.0471e+04  9.4780e+03 -2.1598   0.04058 *  
Employment.Sector..Services    -2.0470e+04  9.4689e+03 -2.1619   0.04040 *  
logGDP                         -3.7461e+00  2.4370e+00 -1.5372   0.13681    
Unemployment_lag1               8.2149e-01  1.0563e-01  7.7771 3.918e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Výsledky po použití robustných smerodajných chýb

Po použití robustných smerodajných chýb sa ukázalo, že náš model je stabilnejší a výsledky sú spoľahlivejšie. V porovnaní s predchádzajúcimi výsledkami vidíme, že oneskorená nezamestnanosť (Unemployment_lag1) zostáva výrazne významná, takže naozaj platí, že nezamestnanosť z minulého roka silno ovplyvňuje tú aktuálnu.

Zaujímavé je, že po úprave sa významnými stali aj premenné Agriculture, Industry a Services. Všetky tri majú veľmi podobné koeficienty aj p-hodnoty. Znamená to, že tieto premenné majú určitý vplyv na nezamestnanosť, ale tento vplyv sme videli až po tom, čo boli chyby „opravené“ robustným prístupom.

Premenná logGDP zostala nevýznamná, takže HDP samotné nemá v našom modeli až taký veľký vplyv, keď už berieme do úvahy oneskorenú nezamestnanosť.

Celkovo môžeme povedať, že robustné smerodajné chyby nám pomohli získať presnejší obraz o tom, ktoré premenné sú skutočne dôležité, aj keď model nemá úplne ideálne vlastnosti.

Zhrnutie

V našej analýze sme najskôr odhadli základný regresný model, aby sme zistili, čo ovplyvňuje mieru nezamestnanosti na Slovensku. Pôvodný model mal však problémy s autokoreláciou, čo sme videli v ACF grafe a potvrdili pomocou Durbin–Watsonovho a Breusch–Godfreyovho testu. To znamenalo, že model nezachytával časový vývoj dostatočne dobre a chyby sa v čase opakovali.

Preto sme vyskúšali Koyckov model, kde sme pridali nezamestnanosť z predchádzajúceho roka. Tento krok výrazne zlepšil kvalitu modelu a zvýšil R², takže model dokázal vysvetliť oveľa väčšiu časť zmien v nezamestnanosti. Autokorelácia sa síce znížila, ale úplne nezmizla.

Nakoniec sme použili robustné smerodajné chyby, ktoré upravujú výsledky tak, aby boli spoľahlivejšie aj vtedy, keď model nie je úplne ideálny. Po tejto úprave sa ukázalo, ktoré premenné sú naozaj dôležité a výsledky pôsobili stabilnejšie.

Celkovo môžeme povedať, že náš model po úpravách lepšie zachytáva vývoj nezamestnanosti v čase. Stále však platí, že časové rady sú špecifické a často si vyžadujú pokročilejšie modely, aby sme úplne odstránili autokoreláciu.