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 reklama: výdavky na reklamu dnes môžu ovplyvniť predaje dnes, ale aj o mesiac alebo o štvrťrok neskôr.
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_austria <- udaje_svet[udaje_svet$Country == "Austria",]
udaje_austria
## Country Region Year Infant_deaths Under_five_deaths
## 48 Austria European Union 2014 3.1 3.8
## 338 Austria European Union 2012 3.3 4.0
## 458 Austria European Union 2006 3.9 4.7
## 467 Austria European Union 2003 4.3 5.2
## 661 Austria European Union 2013 3.2 3.9
## 1101 Austria European Union 2001 4.5 5.4
## 1136 Austria European Union 2004 4.2 5.1
## 1306 Austria European Union 2015 3.0 3.7
## 1455 Austria European Union 2002 4.4 5.3
## 1576 Austria European Union 2000 4.6 5.5
## 1606 Austria European Union 2009 3.7 4.5
## 1780 Austria European Union 2005 4.1 4.9
## 2517 Austria European Union 2008 3.8 4.5
## 2751 Austria European Union 2011 3.5 4.2
## 2782 Austria European Union 2007 3.8 4.6
## 2807 Austria European Union 2010 3.6 4.3
## Adult_mortality Alcohol_consumption Hepatitis_B Measles BMI Polio
## 48 65.8150 12.40 98 87 25.5 98
## 338 69.0215 12.10 92 78 25.5 92
## 458 78.3930 12.58 83 61 25.3 83
## 467 87.6535 12.31 83 46 25.1 84
## 661 67.3475 12.10 95 82 25.5 95
## 1101 91.5155 12.35 44 34 25.0 83
## 1136 84.8235 12.17 83 47 25.2 83
## 1306 64.6875 11.60 93 88 25.5 93
## 1455 89.4275 12.09 81 39 25.1 82
## 1576 94.8820 13.35 33 35 25.0 71
## 1606 76.5700 11.87 83 64 25.4 83
## 1780 83.4800 12.14 86 54 25.2 86
## 2517 74.9395 12.03 83 62 25.3 83
## 2751 71.8545 11.90 89 73 25.4 89
## 2782 78.0555 12.54 85 56 25.3 85
## 2807 74.5020 12.10 86 69 25.4 86
## Diphtheria Incidents_HIV GDP_per_capita Population_mln
## 48 98 0.08 44245 8.55
## 338 92 0.08 44550 8.43
## 458 83 0.08 42496 8.27
## 467 84 0.08 39815 8.12
## 661 95 0.08 44299 8.48
## 1101 84 0.08 39185 8.04
## 1136 83 0.08 40651 8.17
## 1306 93 0.08 44196 8.64
## 1455 83 0.08 39636 8.08
## 1576 81 0.08 38843 8.01
## 1606 83 0.08 42655 8.34
## 1780 86 0.08 41281 8.23
## 2517 83 0.08 44440 8.32
## 2751 89 0.08 44451 8.39
## 2782 85 0.08 43938 8.30
## 2807 86 0.08 43335 8.36
## Thinness_ten_nineteen_years Thinness_five_nine_years Schooling
## 48 1.8 2.0 12.1
## 338 1.8 2.0 11.8
## 458 1.7 1.9 10.1
## 467 1.7 1.9 9.6
## 661 1.8 2.0 11.9
## 1101 1.7 1.9 9.2
## 1136 1.7 1.9 9.8
## 1306 1.9 2.1 12.1
## 1455 1.7 1.9 9.4
## 1576 1.7 1.9 9.0
## 1606 1.7 1.9 11.7
## 1780 1.7 1.9 9.9
## 2517 1.7 1.9 11.5
## 2751 1.7 2.0 11.8
## 2782 1.7 1.9 11.6
## 2807 1.7 1.9 11.8
## Economy_status_Developed Economy_status_Developing Life_expectancy
## 48 1 0 81.5
## 338 1 0 80.9
## 458 1 0 79.9
## 467 1 0 78.6
## 661 1 0 81.1
## 1101 1 0 78.6
## 1136 1 0 79.2
## 1306 1 0 81.2
## 1455 1 0 78.7
## 1576 1 0 78.1
## 1606 1 0 80.3
## 1780 1 0 79.3
## 2517 1 0 80.4
## 2751 1 0 81.0
## 2782 1 0 80.2
## 2807 1 0 80.6
udaje_austria$y <- udaje_austria$Life_expectancy
udaje_austria$x <- udaje_austria$Schooling
#
fit_klasicky <- lm(y ~ x, data = udaje_austria)
summary(fit_klasicky)
##
## Call:
## lm(formula = y ~ x, data = udaje_austria)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45552 -0.23790 0.06823 0.15388 0.57232
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.38688 0.70784 99.44 < 2e-16 ***
## x 0.88523 0.06499 13.62 1.81e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2979 on 14 degrees of freedom
## Multiple R-squared: 0.9298, Adjusted R-squared: 0.9248
## F-statistic: 185.5 on 1 and 14 DF, p-value: 1.81e-09
dwtest(fit_klasicky)
##
## Durbin-Watson test
##
## data: fit_klasicky
## DW = 1.9684, p-value = 0.4796
## alternative hypothesis: true autocorrelation is greater than 0
udaje_austria$y_lag <- c(NA, head(udaje_austria$y, -1)) # posunie to o jedno obdobie dozadu
udaje2 <- na.omit(udaje_austria) # 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.47316 -0.19549 -0.00457 0.19436 0.52919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.91942 5.72774 11.509 7.69e-08 ***
## x 0.86434 0.06620 13.057 1.88e-08 ***
## y_lag 0.05836 0.06993 0.835 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2883 on 12 degrees of freedom
## Multiple R-squared: 0.9345, Adjusted R-squared: 0.9236
## F-statistic: 85.64 on 2 and 12 DF, p-value: 7.876e-08
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 65.91942331
## 2 beta 0.86434298
## 3 lambda 0.05836244
## 4 alpha = c/(1-lambda) 70.00509098
## 5 dlhodobý efekt beta/(1-lambda) 0.91791472
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) 65.919423 5.493975 11.9985 4.843e-08 ***
## x 0.864343 0.047819 18.0753 4.530e-10 ***
## y_lag 0.058362 0.067142 0.8692 0.4018
## ---
## 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 = 2.3363, p-value = 0.7469
## 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ú: