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)
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.
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. \]
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.
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\).
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.
# 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
#
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
dwtest(fit_klasicky)
##
## Durbin-Watson test
##
## data: fit_klasicky
## DW = 1.7021, p-value = 0.2213
## alternative hypothesis: true autocorrelation is greater than 0
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
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
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.
Koyckova transformácia je užitočná, keď chceme odhadnúť distribuované oneskorenia bez veľkého počtu parametrov. Jej hlavné výhody sú:
Jej hlavné obmedzenia sú: