Špecifikácia modelu

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

Moja databáza obsahuje ročné makroekonomické ukazovatele Slovenskej republiky v období 2010–2022. Ide o kombináciu ukazovateľov inflácie, rastu HDP, verejných financií, bežného účtu a výkonnosti ekonomiky.

Pomocou viacnásobnej lineárnej regresie skúmame, do akej miery je miera nezamestnanosti ovplyvnená makroekonomickými ukazovateľmi, konkrétne ročným rastom HDP (GDP_Growth_Annual) a infláciou meranou indexom CPI (Inflation_CPI). Cieľom je zistiť, ktoré z týchto faktorov štatisticky významne prispievajú k vysvetleniu miery nezamestnanosti.

Výsledok tejto regresie je uvedený nižšie:

#######################################################################
# PRIPRAVA UDAJOV
#######################################################################
udaje <- read.csv("data.csv",dec=".",sep=";",header = TRUE)
# výber relevantných premenných
udaje <- udaje[ ,c("Unemployment_Rate","GDP_Growth_Annual", "Inflation_CPI")] 

################################################################################
# ZAKLADNA REGRESIA
################################################################################
attach(udaje)
model <- lm(Unemployment_Rate ~ +1 + GDP_Growth_Annual + Inflation_CPI, data = udaje)
summary(model)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0760 -2.4708 -0.2798  3.2210  4.5429 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)        10.2438     1.8215   5.624  0.00022
GDP_Growth_Annual   0.1789     0.4251   0.421  0.68284
Inflation_CPI      -0.3054     0.3087  -0.989  0.34578
                     
(Intercept)       ***
GDP_Growth_Annual    
Inflation_CPI        
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.478 on 10 degrees of freedom
Multiple R-squared:  0.1318,    Adjusted R-squared:  -0.04185 
F-statistic: 0.759 on 2 and 10 DF,  p-value: 0.4933

Testovať, či je model v správnej funkčnej forme (t. j. či je lineárna špecifikácia vhodná, alebo či by sme mali transformovať premenné, napríklad pomocou logaritmov alebo mocninami), možno vykonať viacerými spôsobmi.

1. Test RESET (test chyby špecifikácie Ramseyho regresnej rovnice - Ramsey Reset Test)

Ide o najštandardnejší formálny test nesprávnej špecifikácie funkčnej formy ale dá sa použiť aj pre prípady, ak sme nešpecifikovali všetky vysvetľujúce premenné.

Myšlienka: Nech pôvodný model má tvar \[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + u_t\] Ak je náš model správne špecifikovaný, potom pridaním mocnín vyrovnaných hodnôt (napr. \(\hat y_t^2\), \(\hat{y}_t^3\)) by sa pôvodný model nemal podstatne zlepšiť, teda budeme testovať pôvodný model uvedený vyššie a model

\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \gamma_2\hat y_t^2 + \gamma_3\hat{y}_t^3 + u_t\]

Budeme testovať hypotézu

\(H_0:\) model je správne špecifikovaný (\(\gamma_2 = \gamma_3 = 0\))

oproti

\(H_1:\) model je nesprávne špecifikovaný (\(\gamma_2 \ne 0 \quad \text{alebo} \quad \gamma_3 \ne 0\))

# Suppose your model is:
model <- lm(Unemployment_Rate ~ +1 + GDP_Growth_Annual + Inflation_CPI, data = udaje)

# RESET test from 'lmtest' package:
library(lmtest)
resettest(model)

    RESET test

data:  model
RESET = 0.70943, df1 = 2, df2 = 8, p-value = 0.5204

Interpretácia:

Na základe výsledkov RESET testu (p-hodnota = 0.5204) nezamietame nulovú hypotézu o správnej špecifikácii modelu. Model teda pravdepodobne neobsahuje štrukturálnu chybu a lineárna forma sa javí ako vhodná.

2. Grafická analýza

Graf Residuals vs. Fitted

Grafická analýza vzťahu medzi vyrovnanými hodnotami náhodnej premennej a rezíduami:

plot(model, which = 1)

Graf rezíduí voči vyrovnaným hodnotám ukazuje mierne systematické zakrivenie, čo naznačuje, že rezíduá nemajú úplne náhodný charakter. Tento vzor môže poukazovať na nelinearitu v modeli alebo na chýbajúcu vysvetľujúcu premennú.

Grafy C+R

Táto analýza nám môže pomôcť pri hľadaní odpovede na otázku, ktorú premennú by sme mali transfomovať pomocou nejakej známej funkcie. Vychádzajme z pôvodnej rovnice

\[y_t = \beta_0 + \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] Túto rovnicu najprv odhadneme a potom vykresľujeme grafy, kde výraz component+residual (C+R) plot vykresľuje na zvislej osi \(\hat{\beta}_ix_{ti}+e_t\) a na vodorovnej osi vykresľuje hodnoty \(x_{ti}\)

Tieto grafy pomáhajú identifikovať nelineárne vzťahy pre každý regresor:

car::crPlots(model)

Component + Residual grafy naznačujú, že v mojom modeli sú prítomné nelineárne vzťahy. Výrazný odklon od lineárnej špecifikácie vidíme pri premenných Inflation_CPI a GDP_Growth_Annual, kde sa ružová krivka podstatne odchyľuje od modrej priamky. To naznačuje, že lineárny model nemusí správne zachytávať ich vplyv na nezamestnanosť a môže byť vhodné zvážiť transformáciu týchto premenných.

3. Nelineárna špecifikácia

Častokrát môžeme aj zložitejšie nelineárne vzťahy modelovať s pomocou ich aproximácie s polynómom, teda v prípade kvadratických členov

\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \dots + \gamma_i\hat x_{ik}^2 + \dots + \gamma_j\hat x_{jk}^2 + \dots + u_t\] Príklad na túto modifikáciu uvidíme nižšie.

3. Porovnanie základného a modifikovaného modelu

Predpokladajme, že sa pri nelineárnych úpravách pôvodnej rovnice dostaneme k zavedeniu kvadrátu premenných GDP_Growth_Annual a Inflation_CPI, nakoľko nás k tomu motivuje Component + Residual Plot uvedený vyššie. Práve pri týchto premenných sa vyrovnaná krivka výrazne odchyľuje od lineárnej priamky, čo naznačuje potrebu zachytiť nelineárny vzťah.

Ak má transformovaný model vyšší upravený koefcient determinácie \(R^2_{adj}\) a pri RESET test prijmeme alternatívnu hypotézu, odporúčame si výsledky potvrdiť s pomocou Anova testu oboch modelov a prípadne opakovaného Reset Testu uplatneneného na nelineárne transformovaný model

model <- lm(Unemployment_Rate ~ +1 + GDP_Growth_Annual + Inflation_CPI)
model_kvadr <- lm(Unemployment_Rate ~ +1 + GDP_Growth_Annual + I(GDP_Growth_Annual^2) + Inflation_CPI + I(Inflation_CPI^2))
summary(model_kvadr)

Call:
lm(formula = Unemployment_Rate ~ +1 + GDP_Growth_Annual + I(GDP_Growth_Annual^2) + 
    Inflation_CPI + I(Inflation_CPI^2))

Residuals:
   Min     1Q Median     3Q    Max 
-4.369 -2.450 -0.166  3.813  4.252 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            10.210105   2.300219   4.439  0.00217
GDP_Growth_Annual       0.368454   0.795945   0.463  0.65576
I(GDP_Growth_Annual^2) -0.041092   0.138746  -0.296  0.77464
Inflation_CPI          -0.250819   1.000185  -0.251  0.80831
I(Inflation_CPI^2)     -0.005236   0.076596  -0.068  0.94717
                         
(Intercept)            **
GDP_Growth_Annual        
I(GDP_Growth_Annual^2)   
Inflation_CPI            
I(Inflation_CPI^2)       
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.867 on 8 degrees of freedom
Multiple R-squared:  0.1416,    Adjusted R-squared:  -0.2877 
F-statistic: 0.3298 on 4 and 8 DF,  p-value: 0.8506

anova(model, model_kvadr)

Analysis of Variance Table

Model 1: Unemployment_Rate ~ +1 + GDP_Growth_Annual + Inflation_CPI
Model 2: Unemployment_Rate ~ +1 + GDP_Growth_Annual + I(GDP_Growth_Annual^2) + 
    Inflation_CPI + I(Inflation_CPI^2)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     10 120.96                           
2      8 119.60  2    1.3594 0.0455 0.9558

resettest(model_kvadr)

    RESET test

data:  model_kvadr
RESET = 2.7558, df1 = 2, df2 = 6, p-value = 0.1416

Po rozšírení pôvodného modelu o kvadratické členy premenných GDP_Growth_Annual a Inflation_CPI sa ukazuje, že ani jedna z novopridaných premenných nie je štatisticky významná. Rovnako však ostávajú nevýznamné aj pôvodné vysvetľujúce premenné.

Upravený koeficient determinácie \(R^2_{adj}\) sa navyše znížil (z pôvodného nízkeho čísla na ešte nižšiu hodnotu – približne –0.29), čo znamená, že rozšírený model vysvetľuje variabilitu nezávislej premennej ešte horšie než základný model.

Porovnanie modelov pomocou ANOVA testu ukazuje, že pridanie kvadratických členov neviedlo k zlepšeniu modelu (p-hodnota ≈ 0.956). Rozšírený model teda nie je štatisticky odlišný od základného a neponúka lepší opis dát.

Aj podľa RESET testu (p ≈ 0.142) nemáme dôkazy o tom, že by model trpel nesprávnou funkčnou formou po pridaní kvadratických členov. To znamená, že rozšírenie modelu o druhé mocniny nezlepšilo špecifikáciu, ale zároveň ani nepoukazuje na novú nelinearitu.

Tento modifikovaný model teda nepriniesol žiadne zlepšenie – kvadratické členy sú štatisticky nevýznamné a hodnota \(R^2_{adj}\) klesla. Na rozdiel od príkladu v šablóne, kde sa po nelineárnej úprave model zlepšil, v mojom prípade doplnenie nelineárnych premenných nebolo prínosné.

Z tohto dôvodu nemá zmysel kvadratické členy ďalej používať a vhodnejšie je ostať pri pôvodnom lineárnom modeli.

5. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie

Predpokladajme, že máme dummy premennú \(D_t\), ktorá obsahuje len nuly a jedničky. Takáto premenná sa dá využiť modelovanie zlomov, a to napr.

1. zlom v autonómnom člene \(\beta_0\) a to nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_D D+ \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] čo interpretujeme ako posun regresnej priamky (regresnej nadroviny) o \(\beta_D\) jednotiek pozdĺž zvislej osi a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\)
2. zlom v sklone regresnej priamky (nadroviny) a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\), čo dosiahneme nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_1 x_{t1} + \dots + \beta_{i}x_{ti} + \beta_{Di}D_tx_{ti}+ \dots + \beta_k x_{tk} + u_t\] kde teda sklon priamky pozdĺž premenne \(x_{ti}\) je \(\beta_i\) ale len v prípade \(D_t=0\), inak je ten sklon rovný \(\beta_i+\beta_{D_i}\).

V mojom prípade hľadáme zlom v jednej z vysvetľujúcich premenných GDP_Growth_Annual alebo Inflation_CPI.

Podľa Component + Residual grafov sa najvýraznejší odklon od linearity objavuje pri premenných Inflation_CPI, kde krivka naznačuje zmenu smerovania približne pri hodnotách okolo 5–6 % inflácie.

Preto zavediem dummy premennú DUM, ktorá vyjadruje, či inflácia prekročila túto hodnotu.

Definícia dummy premennej:

ak je Inflation_CPI < 6, potom DUM = 0

ak je Inflation_CPI ≥ 6, potom DUM = 1

Takto dokážeme preskúmať, či sa intercept (autonómny člen) modelu mení pri vyššej úrovni inflácie.

udaje$DUM <- ifelse(udaje$Inflation_CPI < 6, 0, 1)
modelD_auto <- lm(Unemployment_Rate ~ +1 + DUM + GDP_Growth_Annual + Inflation_CPI,data=udaje ) 
summary(modelD_auto)

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

Residuals:
   Min     1Q Median     3Q    Max 
-4.217 -2.530  0.000  3.296  4.295 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)        10.0627     2.1491   4.682  0.00115 **
DUM                -1.6582     8.9072  -0.186  0.85644   
GDP_Growth_Annual   0.1772     0.4473   0.396  0.70122   
Inflation_CPI      -0.1833     0.7320  -0.250  0.80793   
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.659 on 9 degrees of freedom
Multiple R-squared:  0.1351,    Adjusted R-squared:  -0.1532 
F-statistic: 0.4687 on 3 and 9 DF,  p-value: 0.7114

Na základe odhadu modelu môžeme konštatovať, že zavedenie dummy premennej DUM, ktorá reprezentuje možné zmeny v úrovni nezamestnanosti pri vyšších hodnotách inflácie, neprinieslo žiadny štatisticky významný efekt. Koeficient pri premennej DUM je síce záporný (–1.66), no jeho p-hodnota (0.856) je vysoko nad bežnými hladinami významnosti, čo znamená, že rozdiel v autonómnom člene medzi oboma skupinami (Inflation_CPI < 6 a ≥ 6) nie je štatisticky preukázateľný.

Rovnako aj ostatné premenné (GDP_Growth_Annual a Inflation_CPI) ostávajú v tomto modeli štatisticky nevýznamné a celková vysvetľovacia schopnosť modelu je nízka (R² = 0.1351, upravené R² dokonca negatívne).

V tomto prípade sa teda nepotvrdilo, že by existoval signifikantný zlom v autonómnom člene regresnej rovnice pri oddeľovaní krajín podľa úrovne inflácie. Zavedenie dummy premennej teda nezlepšuje model a nenaznačuje odlišné správanie nezamestnanosti v závislosti od toho, či inflácia prekročila stanovenú hranicu.

modelD_sklon <- lm(Unemployment_Rate ~ +1  + Inflation_CPI + I(DUM*Inflation_CPI) + GDP_Growth_Annual, data=udaje) 
summary(modelD_sklon)

Call:
lm(formula = Unemployment_Rate ~ +1 + Inflation_CPI + I(DUM * 
    Inflation_CPI) + GDP_Growth_Annual, data = udaje)

Residuals:
   Min     1Q Median     3Q    Max 
-4.217 -2.530  0.000  3.296  4.295 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)             10.0627     2.1491   4.682  0.00115
Inflation_CPI           -0.1833     0.7320  -0.250  0.80793
I(DUM * Inflation_CPI)  -0.1298     0.6973  -0.186  0.85644
GDP_Growth_Annual        0.1772     0.4473   0.396  0.70122
                         
(Intercept)            **
Inflation_CPI            
I(DUM * Inflation_CPI)   
GDP_Growth_Annual        
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.659 on 9 degrees of freedom
Multiple R-squared:  0.1351,    Adjusted R-squared:  -0.1532 
F-statistic: 0.4687 on 3 and 9 DF,  p-value: 0.7114

V tomto modeli sme zaviedli interakčný člen I(DUM * Inflation_CPI), ktorý umožňuje, aby sa vplyv inflácie na mieru nezamestnanosti líšil medzi obdobiami s nižšou infláciou (DUM = 0) a obdobiami s vyššou infláciou (DUM = 1).

Z výsledkov však vyplýva, že:

koeficient pri Inflation_CPI je štatisticky nevýznamný (p = 0.80793),

interakčný koeficient I(DUM * Inflation_CPI) je taktiež úplne nevýznamný (p = 0.85644).

To znamená, že neexistuje dôkaz o tom, že by sa sklon regresnej priamky vzhľadom na infláciu líšil medzi obdobiami s nižšou a vyššou infláciou. Inými slovami, inflácia nemá štatisticky významny vplyv na mieru nezamestnanosti ani v jednej z týchto dvoch skupín období a rozdiel v sklonoch sa nepotvrdil.

Okrem toho:

upravený koeficient determinácie (−0.1532) naznačuje, že pridaním interakčného termu sa model nezlepšil, ale dokonca mierne zhoršil,

jediný štatisticky významný parameter je intercept,

celkový F-test (p = 0.7114) ukazuje, že model ako celok nie je štatisticky významný.

V tomto prípade sa teda nepotvrdilo, že by existoval štatisticky významný zlom v sklone regresnej funkcie. Interakcia medzi infláciou a dummy premennou je nevýznamná, takže vplyv inflácie na nezamestnanosť sa medzi obdobiami s nižšou a vyššou infláciou nelíši. Zavedením tejto úpravy sa model nezlepšil a jeho vysvetľovacia sila dokonca klesla.

anova(model, modelD_sklon)

Analysis of Variance Table

Model 1: Unemployment_Rate ~ +1 + GDP_Growth_Annual + Inflation_CPI
Model 2: Unemployment_Rate ~ +1 + Inflation_CPI + I(DUM * Inflation_CPI) + 
    GDP_Growth_Annual
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     10 120.96                           
2      9 120.50  1   0.46403 0.0347 0.8564

resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 3.7203, df1 = 2, df2 = 7, p-value = 0.0793

ANOVA analýza ukázala, že pridanie interakčného členu I(DUM * Inflation_CPI) do modelu nezlepšilo jeho vysvetľovaciu schopnosť (p = 0.8564). RESET test naznačuje miernu možnosť nesprávnej špecifikácie modelu (p = 0.0793), preto je vhodné zvážiť doplnenie ďalších premenných alebo nelineárnych členov.