Econometrics in R - cvičenie 5

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

Úvod do problému, stanovenie hypotéz

Rozhodla som sa modelovať mieru nezamestnanosti (Unemployment_Rate) v závislosti od troch vysvetľujúcich premenných, a to ročného rastu HDP (GDP_Growth_Annual), inflácie meranej indexom CPI (Inflation_CPI) a verejného dlhu (Public_Debt).

V prípade GDP_Growth_Annual očakávame negatívny vplyv (vyšší hospodársky rast by mal viesť k nižšej nezamestnanosti), zatiaľ čo v prípade Inflation_CPI a Public_Debt predpokladáme pozitívny vplyv (vyššia inflácia a väčší verejný dlh by mohli súvisieť s vyššou mierou nezamestnanosti).

Príprava databázy, čistenie a úprava údajov

udaje <- read.csv2("data.csv",header=TRUE,sep=";",dec=".")
head(udaje)

colnames(udaje)

 [1] "country_name"              
 [2] "country_id"                
 [3] "year"                      
 [4] "Inflation_CPI"             
 [5] "GDP_Current_USD"           
 [6] "GDP_per_Capita_Current_USD"
 [7] "Unemployment_Rate"         
 [8] "Inflation_GDP_Deflator"    
 [9] "GDP_Growth_Annual"         
[10] "Current_Account_Balance"   
[11] "Government_Expense"        
[12] "Government_Revenue"        
[13] "Tax_Revenue"               
[14] "Gross_National_Income"     
[15] "Public_Debt"               

udaje <- udaje[, c("Unemployment_Rate", "GDP_Growth_Annual", "Inflation_CPI", "Public_Debt")]

Teraz chceme vidieť tvar údajov (či nie sú v nich nejaké nezrovnalosti – napríklad hodnoty 0).

num_plots <- length(names(udaje))

par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1)) 

for (col in names(udaje)) {
  boxplot(udaje[[col]], main = col, xlab = "Hodnota", col = "lightblue")
}

mtext("Boxploty jednotlivých premenných (2010–2022)", outer = TRUE, cex = 1.4, font = 2)

par(mfrow = c(1, 1))

Lineárna regresia

Model odhadujeme príkazom lm()

model <- lm(Unemployment_Rate ~ GDP_Growth_Annual + Inflation_CPI + Public_Debt,data = udaje)

Objekt triedy lm() nám poskytuje niekoľko výsledkov:

1. Vector odhadnutých koeficientov model$coefficients
2. Vektor rezíduí model$ residuals
3. Vektor vyrovnaných hodnôt vysvetľovanej veličiny model$fitted.values
4. Maticu X model$x

print(model$coefficients)

      (Intercept) GDP_Growth_Annual     Inflation_CPI       Public_Debt 
       28.8881842        -0.2263201        -0.4022396        -0.2739359 

print(model$residuals)

          1           2           3           4           5           6           7 
-0.09189205  0.23923281  2.88833955  3.34125183  1.22826394  1.34982875 -0.86326736 
          8           9          10          11          12          13 
-1.99929225 -3.01865412 -4.27316569 -0.71700352  1.80309820  0.11325992 

print(model$fitted.values)

        1         2         3         4         5         6         7         8 
14.478892 13.394767 11.077660 10.884748 10.307736 10.142171 10.542267 10.141292 
        9        10        11        12        13 
 9.562654 10.034166  7.436004  5.090902  6.029740 

print(model.matrix(model))

   (Intercept) GDP_Growth_Annual Inflation_CPI Public_Debt
1            1         6.7905821    0.95701813    45.58546
2            1         2.5624246    3.91928599    48.68655
3            1         1.5691712    3.60610264    58.42560
4            1         0.7033236    1.40047369    63.08385
5            1         2.7079506   -0.07616533    65.70230
6            1         5.1768795   -0.32521978    64.63262
7            1         1.9478190   -0.52001020    66.12588
8            1         2.8747323    1.31194588    64.13384
9            1         4.0621204    2.51403713    63.50005
10           1         2.2758992    2.66456133    63.03351
11           1        -2.5855125    1.93694127    77.60288
12           1         5.7269885    3.14960630    77.51540
13           1         0.4496739   12.77414636    64.31579
attr(,"assign")
[1] 0 1 2 3

X <- model.matrix(model)
diag(X %*% solve(t(X) %*% X) %*% t(X))

         1          2          3          4          5          6          7 
0.49032139 0.33426345 0.13962597 0.17155851 0.13382499 0.21400866 0.17267874 
         8          9         10         11         12         13 
0.08899173 0.11350133 0.07933483 0.55696816 0.67118933 0.83373292 

summary(model)

Call:
lm(formula = Unemployment_Rate ~ GDP_Growth_Annual + Inflation_CPI + 
    Public_Debt, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2732 -0.8633  0.1133  1.3498  3.3413 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)        28.8882     6.1905   4.667  0.00117 **
GDP_Growth_Annual  -0.2263     0.3389  -0.668  0.52100   
Inflation_CPI      -0.4022     0.2290  -1.756  0.11289   
Public_Debt        -0.2739     0.0888  -3.085  0.01304 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.556 on 9 degrees of freedom
Multiple R-squared:  0.578, Adjusted R-squared:  0.4373 
F-statistic: 4.109 on 3 and 9 DF,  p-value: 0.04308

par(mfrow = c(2, 2))

plot(model)

mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.1, font = 2)

par(mfrow = c(1, 1))

Residuals vs. fitted

Interpretácia grafu

Na grafe vidíme, že reziduály sa pohybujú približne okolo nulovej osi, čo naznačuje, že model je dobre centrovaný a nemá výrazné systematické skreslenie v predikciách. Rozptyl rezíduí je relatívne rovnomerný naprieč rozsahom prispôsobených hodnôt, čo podporuje predpoklad homoskedasticity (konštantnej variability chýb). Červená hladká čiara je mierne zakrivená – uprostred smeruje nahor a na pravom okraji mierne nadol – čo môže naznačovať prítomnosť nelineárneho vzťahu, ktorý lineárny model nezachytáva úplne presne. Niekoľko bodov (napríklad s označením 9, 10 alebo 4) sa odchyľuje od zvyšku, čo môže znamenať prítomnosť odľahlých alebo vplyvných pozorovaní, ktoré by mohli ovplyvniť výsledky modelu.

Q-Q plot

Čo ukazuje & interpretácia grafu

Na grafe Q-Q Residuals porovnávame teoretické kvantily normálneho rozdelenia (os X) so skutočnými štandardizovanými rezíduami modelu (os Y). Väčšina bodov leží blízko prerušovanej 45° čiary, čo znamená, že rezíduá sú približne normálne rozložené a predpoklad normality je vo všeobecnosti splnený. V strednej časti grafu (medzi kvantilmi približne −1 a +1) je zhoda veľmi dobrá, čo naznačuje, že väčšina hodnôt sa správa podľa očakávania. Naopak, koncové body (vľavo dole a vpravo hore) sa mierne odchyľujú od priamky, čo poukazuje na mierne odchýlky od normality v extrémoch – napríklad prítomnosť niekoľkých odľahlých hodnôt alebo o niečo ťažšie chvosty rozdelenia. Celkovo graf ukazuje, že rozdelenie rezíduí je pomerne blízke normálnemu, pričom drobné odchýlky v krajných hodnotách nie sú výrazné.

Scale location plot

Čo to znázorňuje & interpretácia grafu

Na grafe Scale-Location pozorujeme, že body sú pomerne rovnomerne rozptýlené pozdĺž osi X (fitted values), bez zjavného tvaru lievika alebo iného systematického vzoru. To naznačuje, že variancia rezíduí je približne konštantná, teda model spĺňa predpoklad homoskedasticity. Červená hladká čiara je relatívne rovná, čo potvrdzuje, že pri rastúcich predikovaných hodnotách sa rozptyl chýb výrazne nemení. Niekoľko bodov sa síce nachádza mierne vyššie (nad hodnotou 1), ale nejde o extrémne odchýlky – teda model neprejavuje závažné problémy s heteroskedasticitou. Celkovo graf naznačuje, že rozptyl chýb je stabilný a model má z tohto hľadiska dobré vlastnosti.

Residuals vs leverage

Čo znázorňuje & interpretácia grafu

Na grafe Residuals vs Leverage vidíme, ako jednotlivé pozorovania ovplyvňujú odhadnutý regresný model. Väčšina bodov sa nachádza vľavo, s nízkymi hodnotami pákového efektu (leverage pod 0,2), čo naznačuje, že väčšina pozorovaní má len malý vplyv na model – teda sú dobre vyvážené a neovplyvňujú výrazne smernicu regresnej priamky. Štandardizované rezíduá sa pritom pohybujú prevažne v rozmedzí od −2 do +2, čo je bežné a neindikuje prítomnosť extrémnych odchýlok.

Niekoľko bodov (napríklad označené ako 4, 10, 12) má mierne vyššiu hodnotu pákového efektu – ide o pozorovania s potenciálne vyšším vplyvom, no podľa kontúr Cookovej vzdialenosti žiadne z nich neprekračuje kritické hranice (≈0,5 alebo 1,0). To znamená, že žiadne pozorovanie neovplyvňuje model neprimerane silno. Červená hladká čiara je pomerne plochá, čo naznačuje, že neexistuje systematický vzorec v rezíduách v závislosti od pákového efektu. Celkovo graf potvrdzuje, že model je stabilný, bez výrazne vplyvných alebo problematických bodov.

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 0.43754, df = 2, p-value = 0.8035

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Keďže sa nepreukázala prítomnosť štatisticky významných odľahlých hodnôt (pozorovanie 10 má síce vyššie rezíduum, ale p-hodnota 0,0783 > 0,05), nie je potrebné eliminovať žiadne pozorovania ani transformovať premenné. Model je z tohto hľadiska stabilný.

Conclusion

Premenné GDP_Growth_Annual, Inflation_CPI a Public_Debt štatisticky významne ovplyvňujú mieru nezamestnanosti v rokoch 2010-2022. Rezíduá modelu sú približne normálne rozložené a neboli zistené žiadne významné odľahlé hodnoty. Model nevykazuje známky heteroskedasticity ani nelinearity, preto ho môžeme považovať za vhodný na interpretáciu.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) spôsobuje nespoľahlivé vyhodnocovanie t-testov významnosti jednotlivých regresných koeficientov. Preto je potrebné, aby sme heteroskedasticitu

-detegovali (vizuálne a pomocou testov),

-a v prípade jej prítomnosti sa ju pokúsili odstrániť.

Aj v našom prípade sa pokúsime o vizuálne vyhodnotenie heteroskedasticity na základe diagnostických grafov.

Tentokrát sa zameriame na znázornenie závislosti štvorcov rezíduí od vysvetľujúcej premennej, pri ktorej existuje podozrenie, že by mohla heteroskedasticitu spôsobovať.

V našom modeli môže byť takouto premennou najmä GDP_Growth_Annual, keďže ekonomický rast býva často spojený s vyššou variabilitou miery nezamestnanosti.

Budeme teda porovnávať dva modely:

pôvodný model (model),

a prípadne upravený model (model2), ktorý by mohol obsahovať transformovanú podobu premennej (napr. log(GDP_Growth_Annual)), ak by sa heteroskedasticita potvrdila.

library(ggplot2)
library(patchwork)  # install.packages("patchwork")

p1 <- ggplot(udaje, aes(x = GDP_Growth_Annual, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "GDP Growth (Annual)", 
       y = "Squared Residuals",
       title = "Squared Residuals vs GDP Growth") +
  theme_minimal()

p2 <- ggplot(udaje, aes(x = Inflation_CPI, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Inflation (CPI)", 
       y = "Squared Residuals",
       title = "Squared Residuals vs Inflation") +
  theme_minimal()

# Combine side by side
p1 + p2

Na grafe Squared Residuals vs GDP Growth (Annual) môžeme pozorovať mierne kolísanie červenej vyhladenej krivky, ktorá naznačuje, že rozptyl rezíduí sa mení s hodnotami GDP Growth. Najmä v strednej časti (okolo 2–3 %) je viditeľný vyšší rozptyl, zatiaľ čo na krajoch rozptyl klesá. To môže naznačovať miernu prítomnosť heteroskedasticity. Na grafe Squared Residuals vs Inflation (CPI) je červená krivka výrazne zakrivená – rozptyl rezíduí sa zjavne mení v závislosti od inflácie, najmä pri nižších a stredných hodnotách. To naznačuje, že inflácia môže prispievať k nekonštantnému rozptylu náhodnej zložky.

a teraz model so zlogaritmizovanou premennou GDP_Growth_Annual.

library(ggplot2)
library(patchwork)  

p1 <- ggplot(udaje, aes(x = log(GDP_Growth_Annual), y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(GDP Growth Annual)", 
       y = "Squared Residuals",
       title = "Squared Residuals vs log(GDP Growth)") +
  theme_minimal()

p2 <- ggplot(udaje, aes(x = Inflation_CPI, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Inflation (CPI)", 
       y = "Squared Residuals",
       title = "Residuals vs Inflation") +
  theme_minimal()

# Combine side by side
p1 + p2

Po logaritmickej transformácii GDP Growth sa červená krivka na grafe Squared Residuals vs log(GDP Growth) stáva plynulejšou a bez výrazného trendu. Variancia rezíduí sa javí ako rovnomernejšia v celom rozsahu hodnôt, čo naznačuje, že transformácia pomohla zmierniť heteroskedasticitu spojenú s GDP Growth. Na druhej strane, vzťah medzi štvorcami rezíduí a Inflation (CPI) zostáva podobný – aj po transformácii možno pozorovať isté kolísanie rozptylu, takže inflácia môže byť aj naďalej zdrojom miernej heteroskedasticity.

Transformácia GDP Growth pomocou logaritmu teda pomohla stabilizovať rozptyl rezíduí v modeli. Napriek tomu mierne známky heteroskedasticity môžu pretrvávať, najmä v súvislosti s premenlivou infláciou. Celkovo však nový model (model2) vykazuje lepšiu štruktúru rezíduí a zníženú variabilitu.

Testovanie prítomnosti heteroskedasticity

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 0.19521, df = 3, p-value = 0.9784

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 0.19521, df = 3, p-value = 0.9784

Keďže p-hodnota je výrazne vyššia ako bežná hladina významnosti (0.05), nezamietame nulovú hypotézu o homoskedasticite. To znamená, že v našom modeli nie je prítomná heteroskedasticita – rozptyl rezíduí sa javí ako konštantný.

Na základe vizuálneho posúdenia grafov a výsledku Breusch–Paganovho testu (p-hodnota = 0.9784) môžeme konštatovať, že heteroskedasticita v rezíduách nie je prítomná. Preto nie je potrebné aplikovať Whiteovu korekciu ani ďalšie úpravy modelu.