# NENAČÍTAVAM žiadne balíky cez library(), aby to nespadlo.
# Len si overím, či náhodou nie sú dostupné.
have_lmtest <- requireNamespace("lmtest", quietly = TRUE)
have_sandwich <- requireNamespace("sandwich", quietly = TRUE)
have_lmtest
## [1] FALSE
have_sandwich
## [1] FALSE
V tomto cvičení používam údaje z Eurostatu
súbor: une_rt_a__custom_18708117_linear.csv.
Zameriavam sa na ročnú mieru nezamestnanosti na
Slovensku v rokoch 2014–2024.
# načítanie údajov
udaje <- read.csv("une_rt_a__custom_18708117_linear.csv",
stringsAsFactors = FALSE)
str(udaje)
## 'data.frame': 33 obs. of 11 variables:
## $ DATAFLOW : chr "ESTAT:UNE_RT_A(1.0)" "ESTAT:UNE_RT_A(1.0)" "ESTAT:UNE_RT_A(1.0)" "ESTAT:UNE_RT_A(1.0)" ...
## $ LAST.UPDATE: chr "11/09/25 23:00:00" "11/09/25 23:00:00" "11/09/25 23:00:00" "11/09/25 23:00:00" ...
## $ freq : chr "Annual" "Annual" "Annual" "Annual" ...
## $ age : chr "From 15 to 74 years" "From 15 to 74 years" "From 15 to 74 years" "From 15 to 74 years" ...
## $ unit : chr "Percentage of population in the labour force" "Percentage of population in the labour force" "Percentage of population in the labour force" "Percentage of population in the labour force" ...
## $ sex : chr "Total" "Total" "Total" "Total" ...
## $ geo : chr "Czechia" "Czechia" "Czechia" "Czechia" ...
## $ TIME_PERIOD: int 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 ...
## $ OBS_VALUE : num 6.1 5.1 4 2.9 2.2 2 2.6 2.8 2.2 2.6 ...
## $ OBS_FLAG : chr "" "" "" "" ...
## $ CONF_STATUS: logi NA NA NA NA NA NA ...
head(udaje)
## DATAFLOW LAST.UPDATE freq age
## 1 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## 2 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## 3 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## 4 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## 5 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## 6 ESTAT:UNE_RT_A(1.0) 11/09/25 23:00:00 Annual From 15 to 74 years
## unit sex geo TIME_PERIOD
## 1 Percentage of population in the labour force Total Czechia 2014
## 2 Percentage of population in the labour force Total Czechia 2015
## 3 Percentage of population in the labour force Total Czechia 2016
## 4 Percentage of population in the labour force Total Czechia 2017
## 5 Percentage of population in the labour force Total Czechia 2018
## 6 Percentage of population in the labour force Total Czechia 2019
## OBS_VALUE OBS_FLAG CONF_STATUS
## 1 6.1 NA
## 2 5.1 NA
## 3 4.0 NA
## 4 2.9 NA
## 5 2.2 NA
## 6 2.0 NA
Vyberiem len Slovensko a vytvorím časový rad nezamestnanosti:
# filtrovanie na Slovensko a výber stĺpcov
udaje_svk <- udaje[udaje$geo == "Slovakia",
c("TIME_PERIOD", "OBS_VALUE")]
# premenujem stĺpce
names(udaje_svk) <- c("Year", "Unemployment")
# zoradím podľa rokov
udaje_svk <- udaje_svk[order(udaje_svk$Year), ]
udaje_svk
## Year Unemployment
## 23 2014 13.1
## 24 2015 11.5
## 25 2016 9.6
## 26 2017 8.1
## 27 2018 6.5
## 28 2019 5.7
## 29 2020 6.7
## 30 2021 6.8
## 31 2022 6.1
## 32 2023 5.8
## 33 2024 5.3
Základné štatistiky:
summary(udaje_svk$Unemployment)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.300 5.950 6.700 7.745 8.850 13.100
A jednoduchý graf vývoja nezamestnanosti v čase (base R):
plot(udaje_svk$Year, udaje_svk$Unemployment,
type = "b",
xlab = "Rok",
ylab = "Nezamestnanosť (%)",
main = "Miera nezamestnanosti na Slovensku")
Ako jednoduchý základný model používam regresiu miery nezamestnanosti na časový trend:
\[ Unemployment_t = \beta_0 + \beta_1 \, Year_t + e_t. \]
model <- lm(Unemployment ~ Year, data = udaje_svk)
summary(model)
##
## Call:
## lm(formula = Unemployment ~ Year, data = udaje_svk)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0455 -0.6850 0.3918 0.8605 1.9591
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1378.8300 255.0566 5.406 0.000430 ***
## Year -0.6791 0.1263 -5.376 0.000447 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.325 on 9 degrees of freedom
## Multiple R-squared: 0.7625, Adjusted R-squared: 0.7361
## F-statistic: 28.9 on 1 and 9 DF, p-value: 0.0004471
Do grafu doplním regresnú priamku:
plot(udaje_svk$Year, udaje_svk$Unemployment,
xlab = "Rok",
ylab = "Nezamestnanosť (%)",
main = "Nezamestnanosť na Slovensku a lineárny trend")
points(udaje_svk$Year, udaje_svk$Unemployment, pch = 16)
abline(model, lty = 2, col = "blue")
V texte viem stručne opísať:
Year (klesajúci / rastúci
trend),Pri časových radoch ma zaujíma, či sú reziduá nezávislé v čase.
udaje_svk$residuals <- residuals(model)
udaje_svk$fitted <- fitted(model)
head(udaje_svk)
## Year Unemployment residuals fitted
## 23 2014 13.1 1.9590909 11.140909
## 24 2015 11.5 1.0381818 10.461818
## 25 2016 9.6 -0.1827273 9.782727
## 26 2017 8.1 -1.0036364 9.103636
## 27 2018 6.5 -1.9245455 8.424545
## 28 2019 5.7 -2.0454545 7.745455
# reziduá v čase
plot(udaje_svk$Year, udaje_svk$residuals,
type = "b",
xlab = "Rok",
ylab = "Reziduá",
main = "Reziduá základného modelu v čase")
abline(h = 0, col = "red")
# ACF rezíduí (funkcia acf je v base/stats)
acf(udaje_svk$residuals,
main = "ACF rezíduí – základný model")
Z grafov viem vizuálne posúdiť, či reziduá naznačujú autokoreláciu.
Keďže balík lmtest nemusí byť dostupný, vypočítam si
Durbin–Watsonov štatistický ukazovateľ sama z
definície:
\[ DW = \frac{\sum_{t=2}^{T} (e_t - e_{t-1})^2}{\sum_{t=1}^{T} e_t^2}. \]
e <- udaje_svk$residuals
num <- sum((e[2:length(e)] - e[1:(length(e)-1)])^2)
den <- sum(e^2)
DW_stat <- num / den
DW_stat
## [1] 0.4733047
Hodnotu DW môžem slovne okomentovať (≈ 2 bez autokorelácie, < 2 pozitívna, > 2 negatívna autokorelácia).
Ak by bol balík lmtest náhodou k dispozícii, viem použiť
aj klasický test:
lmtest::dwtest(model)
Ak autokorelácia vyzerá ako problém, jednou z možností je pridať do
modelu
oneskorenú hodnotu nezamestnanosti. Vytvorím lag o 1
rok:
\[ Unemployment_t = \alpha_0 + \alpha_1 \, Unemployment_{t-1} + \alpha_2 \, Year_t + u_t. \]
# oneskorená nezamestnanosť – jednoduchý lag o 1 rok
Unemp_lag1 <- c(NA, head(udaje_svk$Unemployment, -1))
udaje_dyn <- data.frame(
Year = udaje_svk$Year,
Unemployment = udaje_svk$Unemployment,
Unemp_lag1 = Unemp_lag1
)
# odstránim prvý riadok s NA
udaje_dyn <- udaje_dyn[!is.na(udaje_dyn$Unemp_lag1), ]
udaje_dyn
## Year Unemployment Unemp_lag1
## 2 2015 11.5 13.1
## 3 2016 9.6 11.5
## 4 2017 8.1 9.6
## 5 2018 6.5 8.1
## 6 2019 5.7 6.5
## 7 2020 6.7 5.7
## 8 2021 6.8 6.7
## 9 2022 6.1 6.8
## 10 2023 5.8 6.1
## 11 2024 5.3 5.8
Odhadnem dynamický model:
model_dyn <- lm(Unemployment ~ Unemp_lag1 + Year,
data = udaje_dyn)
summary(model_dyn)
##
## Call:
## lm(formula = Unemployment ~ Unemp_lag1 + Year, data = udaje_dyn)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8265 -0.2940 -0.1893 0.4054 1.1145
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 54.18923 297.51125 0.182 0.86063
## Unemp_lag1 0.70369 0.17212 4.088 0.00464 **
## Year -0.02605 0.14673 -0.178 0.86413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6549 on 7 degrees of freedom
## Multiple R-squared: 0.9142, Adjusted R-squared: 0.8897
## F-statistic: 37.29 on 2 and 7 DF, p-value: 0.000185
Skontrolujem ACF rezíduí aj pri dynamickom modeli:
acf(residuals(model_dyn),
main = "ACF rezíduí – dynamický model")
Ak je lmtest dostupný, môžem znovu urobiť testy:
lmtest::dwtest(model_dyn)
lmtest::bgtest(model_dyn, order = 1)
Ak sú nainštalované balíky lmtest a
sandwich, viem získať
Newey–West robustné štandardné chyby:
# robustné štandardné chyby pre základný model
lmtest::coeftest(model, vcov. = sandwich::NeweyWest(model))
# robustné štandardné chyby pre dynamický model
lmtest::coeftest(model_dyn, vcov. = sandwich::NeweyWest(model_dyn))
Ak balíky nie sú dostupné, tento chunk sa automaticky preskočí a knit nespadne.
V cvičení som:
Na záver viem stručne zhodnotiť: