rm(list=ls())
library(lmtest)   #  podpora regresie
library(outliers) # analyza odlahlych hodnot (outliers)
library(gptstudio)
library(kableExtra)
library(knitr)
library(dplyr)
library(broom)
library(corrplot)
library(sandwich)
library(lmtest)

1 Motivácia

V empirickej ekonómii často predpokladáme, že vysvetľujúca premenná nemá účinok iba v aktuálnom období, ale aj v budúcnosti. Typickým príkladom je choroba: vie nám ovplyvniť žiovt teraz o mesiac aj o rok.

Všeobecný model distribuovaných oneskorení môžeme zapísať ako

\[ y_t = \alpha + \beta_0 x_t + \beta_1 x_{t-1} + \beta_2 x_{t-2} + \cdots + u_t. \]

Problém je zrejmý: model má nekonečne veľa parametrov \(\beta_j\). Koyckova transformácia rieši tento problém tým, že predpokladá geometricky klesajúce oneskorené účinky.

2 Koyckov predpoklad

Koyckov model predpokladá, že koeficienty distribuovaných oneskorení majú tvar

\[ \beta_j = \beta \lambda^j, \qquad j = 0,1,2,\ldots, \]

kde

\[ |\lambda| < 1. \]

Parameter \(\lambda\) určuje rýchlosť doznievania účinku. Ak je \(\lambda\) blízko nuly, účinok \(x_t\) rýchlo mizne. Ak je \(\lambda\) blízko jednej, účinok je veľmi perzistentný.

Pôvodný model je potom

\[ y_t = \alpha + \beta x_t + \beta\lambda x_{t-1} + \beta\lambda^2 x_{t-2} + \cdots + u_t. \]

3 Odvodenie Koyckovej transformácie

Začnime s nekonečným modelom

\[ y_t = \alpha + \beta \sum_{j=0}^{\infty} \lambda^j x_{t-j} + u_t. \]

Posuňme rovnicu o jedno obdobie dozadu a vynásobme ju \(\lambda\):

\[ \lambda y_{t-1} = \lambda \alpha + \beta \sum_{j=1}^{\infty} \lambda^j x_{t-j} + \lambda u_{t-1}. \]

Odpočítaním druhej rovnice od prvej dostaneme

\[ y_t - \lambda y_{t-1} = \alpha(1-\lambda) + \beta x_t + u_t - \lambda u_{t-1}. \]

Teda

\[ y_t = \alpha(1-\lambda) + \beta x_t + \lambda y_{t-1} + v_t, \]

kde

\[ v_t = u_t - \lambda u_{t-1}. \]

Toto je Koyckova autoregresívna transformácia. Nekonečný model distribuovaných oneskorení sa zmenil na model s oneskorenou závislou premennou.

4 Interpretácia parametrov

Transformovaný model je

\[ y_t = c + \beta x_t + \lambda y_{t-1} + v_t, \]

kde

\[ c = \alpha(1-\lambda). \]

Z toho plynie:

\[ \alpha = \frac{c}{1-\lambda}. \]

Krátkodobý efekt \(x_t\) na \(y_t\) je \(\beta\). Dlhodobý efekt je súčet všetkých oneskorených účinkov:

\[ \beta + \beta\lambda + \beta\lambda^2 + \cdots = \frac{\beta}{1-\lambda}. \]

Mediánové oneskorenie, teda približný počet období, po ktorých sa realizuje polovica celkového dlhodobého účinku, možno vyjadriť ako

\[ m = \frac{\log(0.5)}{\log(\lambda)}. \]

Tento vzorec dáva zmysel pre \(0<\lambda<1\).

5 Ekonometrické upozornenie

Transformácia je algebraicky elegantná, ale prináša dôležitý problém. Nová chyba je

\[ v_t = u_t - \lambda u_{t-1}. \]

Ak je pôvodná chyba \(u_t\) biely šum, potom \(v_t\) je MA(1) proces. Navyše v regresii vystupuje \(y_{t-1}\), ktorý obsahuje \(u_{t-1}\). Preto môže byť \(y_{t-1}\) korelované s \(v_t\). To znamená, že obyčajná metóda najmenších štvorcov nemusí byť všeobecne konzistentná bez dodatočných predpokladov.

V praxi sa preto často používajú robustné štandardné chyby pre časové rady, prípadne inštrumentálne premenné alebo dynamické metódy odhadovania. V kontexte časových radov je dôležité pracovať s pojmami slabej závislosti, miešania a dlhodobej variance. Slabá závislosť znamená, že korelácia medzi vzdialenými pozorovaniami postupne mizne; v asymptotickej teórii potom štandardná variance často nestačí a nahrádza ju dlhodobá variance.

6 Príklad

# import the dataset and create a data.frame udaje
#
udaje_svet <- read.csv("udaje/Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)
##      Country                        Region Year Infant_deaths Under_five_deaths
## 1    Turkiye                   Middle East 2015          11.1              13.0
## 2      Spain                European Union 2015           2.7               3.3
## 3      India                          Asia 2007          51.5              67.9
## 4     Guyana                 South America 2006          32.8              40.5
## 5     Israel                   Middle East 2012           3.4               4.3
## 6 Costa Rica Central America and Caribbean 2006           9.8              11.2
##   Adult_mortality Alcohol_consumption Hepatitis_B Measles  BMI Polio Diphtheria
## 1        105.8240                1.32          97      65 27.8    97         97
## 2         57.9025               10.35          97      94 26.0    97         97
## 3        201.0765                1.57          60      35 21.2    67         64
## 4        222.1965                5.68          93      74 25.3    92         93
## 5         57.9510                2.89          97      89 27.0    94         94
## 6         95.2200                4.19          88      86 26.4    89         89
##   Incidents_HIV GDP_per_capita Population_mln Thinness_ten_nineteen_years
## 1          0.08          11006          78.53                         4.9
## 2          0.09          25742          46.44                         0.6
## 3          0.13           1076        1183.21                        27.1
## 4          0.79           4146           0.75                         5.7
## 5          0.08          33995           7.91                         1.2
## 6          0.16           9110           4.35                         2.0
##   Thinness_five_nine_years Schooling Economy_status_Developed
## 1                      4.8       7.8                        0
## 2                      0.5       9.7                        1
## 3                     28.0       5.0                        0
## 4                      5.5       7.9                        0
## 5                      1.1      12.8                        1
## 6                      1.9       7.9                        0
##   Economy_status_Developing Life_expectancy
## 1                         1            76.5
## 2                         0            82.8
## 3                         1            65.4
## 4                         1            67.0
## 5                         0            81.7
## 6                         1            78.2
udaje_india <- udaje_svet[udaje_svet$Country == "India",]
udaje_india
##      Country Region Year Infant_deaths Under_five_deaths Adult_mortality
## 3      India   Asia 2007          51.5              67.9        201.0765
## 296    India   Asia 2010          45.1              58.2        191.5935
## 528    India   Asia 2012          40.9              52.0        185.2720
## 640    India   Asia 2015          34.9              43.5        177.9465
## 805    India   Asia 2014          36.9              46.2        180.3880
## 1048   India   Asia 2013          38.8              49.0        182.8295
## 1463   India   Asia 2000          66.7              91.7        221.6305
## 1661   India   Asia 2008          49.4              64.6        197.9155
## 1906   India   Asia 2004          57.8              77.8        210.2435
## 2087   India   Asia 2006          53.6              71.1        204.1320
## 2204   India   Asia 2005          55.7              74.4        207.1880
## 2263   India   Asia 2009          47.2              61.4        194.7545
## 2276   India   Asia 2003          60.0              81.1        213.2990
## 2298   India   Asia 2002          62.1              84.5        216.3545
## 2484   India   Asia 2001          64.4              88.1        218.9925
## 2625   India   Asia 2011          43.0              55.0        188.4325
##      Alcohol_consumption Hepatitis_B Measles  BMI Polio Diphtheria
## 3                   1.57          60      35 21.2    67         64
## 296                 2.73          38      38 21.4    76         79
## 528                 3.04          73      36 21.5    79         82
## 640                 3.00          87      69 21.7    86         87
## 805                 3.00          79      60 21.7    84         85
## 1048                3.03          70      51 21.6    82         83
## 1463                0.93          67      34 20.7    57         58
## 1661                1.92          29      33 21.2    69         70
## 1906                1.18          60      34 21.0    58         63
## 2087                1.35          60      34 21.1    66         65
## 2204                1.25          80      34 21.0    65         65
## 2263                2.46          37      34 21.3    73         74
## 2276                1.17          67      34 20.9    57         61
## 2298                1.08          67      34 20.8    58         59
## 2484                0.98          67      34 20.7    58         59
## 2625                2.95          44      27 21.4    79         82
##      Incidents_HIV GDP_per_capita Population_mln Thinness_ten_nineteen_years
## 3             0.13           1076        1183.21                        27.1
## 296           0.10           1244        1234.28                        27.0
## 528           0.08           1347        1265.78                        26.9
## 640           0.06           1606        1310.15                        26.7
## 805           0.07           1503        1295.60                        26.8
## 1048          0.08           1416        1280.84                        26.8
## 1463          0.34            758        1056.58                        27.7
## 1661          0.12           1093        1200.67                        27.0
## 1906          0.17            898        1129.62                        27.2
## 2087          0.14           1015        1165.49                        27.1
## 2204          0.16            954        1147.61                        27.2
## 2263          0.11           1162        1217.73                        27.0
## 2276          0.19            845        1111.52                        27.3
## 2298          0.22            797        1093.32                        27.4
## 2484          0.27            781        1075.00                        27.5
## 2625          0.09           1293        1250.29                        26.9
##      Thinness_five_nine_years Schooling Economy_status_Developed
## 3                        28.0       5.0                        0
## 296                      27.8       5.4                        0
## 528                      27.6       5.6                        0
## 640                      27.3       6.3                        0
## 805                      27.4       6.1                        0
## 1048                     27.5       5.8                        0
## 1463                     28.6       4.4                        0
## 1661                     27.9       5.2                        0
## 1906                     28.2       4.7                        0
## 2087                     28.0       4.9                        0
## 2204                     28.1       4.8                        0
## 2263                     27.8       5.3                        0
## 2276                     28.3       4.7                        0
## 2298                     28.4       4.6                        0
## 2484                     28.5       4.5                        0
## 2625                     27.7       5.4                        0
##      Economy_status_Developing Life_expectancy
## 3                            1            65.4
## 296                          1            66.7
## 528                          1            67.5
## 640                          1            68.6
## 805                          1            68.3
## 1048                         1            67.9
## 1463                         1            62.5
## 1661                         1            65.8
## 1906                         1            64.1
## 2087                         1            64.9
## 2204                         1            64.5
## 2263                         1            66.2
## 2276                         1            63.7
## 2298                         1            63.3
## 2484                         1            62.9
## 2625                         1            67.1
udaje_india$y <- udaje_india$Life_expectancy
udaje_india$x <- udaje_india$Schooling

#

6.1 Klasický odhad

fit_klasicky <- lm(y ~ x, data = udaje_india)
summary(fit_klasicky)
## 
## Call:
## lm(formula = y ~ x, data = udaje_india)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8206 -0.3757  0.1343  0.2496  0.7289 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  48.0741     1.0216   47.06  < 2e-16 ***
## x             3.3883     0.1965   17.24 7.96e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4348 on 14 degrees of freedom
## Multiple R-squared:  0.955,  Adjusted R-squared:  0.9518 
## F-statistic: 297.3 on 1 and 14 DF,  p-value: 7.957e-11

6.1.1 Test na autokoreláciu rezíduí (Durbin-Watson test)

dwtest(fit_klasicky)
## 
##  Durbin-Watson test
## 
## data:  fit_klasicky
## DW = 1.7021, p-value = 0.2213
## alternative hypothesis: true autocorrelation is greater than 0

6.2 Odhad Koyckovej transformácie

udaje_india$y_lag <- c(NA, head(udaje_india$y, -1))   # posunie to o jedno obdobie dozadu
udaje2 <- na.omit(udaje_india)                  # vynecha sa prve pozorovanie (NA)

fit <- lm(y ~ x + y_lag, data = udaje2)
summary(fit)
## 
## Call:
## lm(formula = y ~ x + y_lag, data = udaje2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77802 -0.21632  0.07075  0.27444  0.53069 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 53.94152    3.63189  14.852 4.35e-09 ***
## x            3.55931    0.20683  17.209 8.00e-10 ***
## y_lag       -0.10351    0.06074  -1.704    0.114    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4086 on 12 degrees of freedom
## Multiple R-squared:  0.9659, Adjusted R-squared:  0.9603 
## F-statistic: 170.1 on 2 and 12 DF,  p-value: 1.563e-09

Odhadované parametre sú:

coef_fit <- coef(fit)

c_hat <- coef_fit["(Intercept)"]
beta_hat <- coef_fit["x"]
lambda_hat <- coef_fit["y_lag"]
alpha_hat <- c_hat / (1 - lambda_hat)
long_run_hat <- beta_hat / (1 - lambda_hat)

estimates <- data.frame(
  parameter = c("c", "beta", "lambda", "alpha = c/(1-lambda)", "dlhodobý efekt beta/(1-lambda)"),
  odhad = c(c_hat, beta_hat, lambda_hat, alpha_hat, long_run_hat)
)

estimates
##                        parameter      odhad
## 1                              c 53.9415153
## 2                           beta  3.5593089
## 3                         lambda -0.1035132
## 4           alpha = c/(1-lambda) 48.8816216
## 5 dlhodobý efekt beta/(1-lambda)  3.2254338

7 Robustné štandardné chyby typu Newey-West

Pre dynamické modely časových radov je rozumné ukázať aj HAC štandardné chyby. Balíky sandwich a lmtest nie sú súčasťou úplne základnej inštalácie R, preto ich načítame len vtedy, ak sú dostupné.

  nw_vcov <- sandwich::NeweyWest(fit, lag = 4, prewhite = FALSE, adjust = TRUE)
  lmtest::coeftest(fit, vcov. = nw_vcov)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 53.941515   1.775194  30.386 1.012e-12 ***
## x            3.559309   0.235468  15.116 3.559e-09 ***
## y_lag       -0.103513   0.027782  -3.726  0.002895 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8 Diagnostika rezíduí

res <- residuals(fit)

plot(res, type = "l", main = "Rezíduá z transformovaného modelu",
     xlab = "čas", ylab = "rezíduum")
abline(h = 0, lty = 2)

dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 1.954, p-value = 0.3889
## alternative hypothesis: true autocorrelation is greater than 0

Ak rezíduá vykazujú autokoreláciu, klasické OLS štandardné chyby môžu byť zavádzajúce. To je jeden z dôvodov, prečo sme vyššie uviedli HAC štandardné chyby.

9 Praktické zhrnutie

Koyckova transformácia je užitočná, keď chceme odhadnúť distribuované oneskorenia bez veľkého počtu parametrov. Jej hlavné výhody sú:

  1. nekonečný lagový model sa zredukuje na jednoduchý dynamický regresný model;
  2. parameter \(\lambda\) má jasnú interpretáciu zotrvačnosti;
  3. dlhodobý efekt má jednoduchý tvar \(\beta/(1-\lambda)\).

Jej hlavné obmedzenia sú:

  1. geometrický tvar oneskorení môže byť príliš reštriktívny;
  2. transformovaná chyba má dynamickú štruktúru;
  3. oneskorená závislá premenná môže byť korelovaná s chybou;
  4. pre inferenciu v časových radoch je často potrebné používať robustné štandardné chyby alebo vhodné inštrumenty.