The U.S. Gasoline Market, 52 Yearly Observations, 1953-2004

\(\ln (G/Pop)_t=\beta_0+\beta_1 \ln (Income/Pop)_t+\beta_2 \ln (GasP)_t+\beta_3 \ln (Pnc_t)+\beta_4 \ln (Puc)_t+\beta_5 t+\epsilon_t\)

Nahrání dat

getwd()
## [1] "C:/Users/Lukas/Desktop/Ekonometrie/AKM/cviceni 2017/cviceni 4"
Oil=read.csv2("Oil_AKM.csv", header = TRUE, sep = ";")
class(Oil$GasExp)                                       # zjistíme typ promenné, nekdy vám bude R hlásit, ze je daná metoda pouzitelná napríklad jen pro "numeric" 
## [1] "numeric"

Uprava dat podle modelu

GP=Oil$GasExp            # vybereme data 
PG=Oil$PG
Pnc=Oil$Pnc
Puc=Oil$Puc

GP=log(GP/Oil$Pop)
IP=log(Oil$Income/Oil$Pop)

LPg=log(PG)
LPnc=log(Pnc)
LPuc=log(Puc)
DelkaRady=length(Oil$Year)   # zjistíme delku casové rady
t=1:DelkaRady                # vygenerujeme data od 1 po délku casové rady

Provedeme regresi.

regrese1=lm(GP~IP+LPg+LPnc+LPuc+t)
summary(regrese1)
## 
## Call:
## lm(formula = GP ~ IP + LPg + LPnc + LPuc + t)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15524 -0.02749  0.01093  0.04054  0.06965 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.345224   0.641439 -11.451 4.60e-15 ***
## IP           1.274648   0.235739   5.407 2.22e-06 ***
## LPg          0.973218   0.051990  18.719  < 2e-16 ***
## LPnc        -0.448973   0.199111  -2.255  0.02894 *  
## LPuc         0.059901   0.121538   0.493  0.62446    
## t            0.010614   0.003165   3.353  0.00161 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0583 on 46 degrees of freedom
## Multiple R-squared:  0.9961, Adjusted R-squared:  0.9957 
## F-statistic:  2373 on 5 and 46 DF,  p-value: < 2.2e-16
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
Otestujeme na prítomnost heteroskedasticity.
Jedná se o casovou radu, takze tento postup není správný. Nejprve je nutné testovat náhodnou slozku na prítomnost autokorelace. 
Je to z toho duvodu, ze pouzívané testy na heteroskedasticitu nejsou robustní vuci autokorelaci. Zde si vsak chceme natrenovat testy. 
bptest(regrese1)
## 
##  studentized Breusch-Pagan test
## 
## data:  regrese1
## BP = 4.3255, df = 5, p-value = 0.5036
bptest(GP~IP+LPg+LPnc+LPuc+t)
## 
##  studentized Breusch-Pagan test
## 
## data:  GP ~ IP + LPg + LPnc + LPuc + t
## BP = 4.3255, df = 5, p-value = 0.5036

Pro daný test nemuzeme zamítnout nulovou hypotézu. Tedy nepodarilo se nám prokázat prítomnost heteroskedasticity.

White test pro heteroskedasticitu

residua=regrese1$residuals     # ulozíme si residua z regrese1 
residua_kvadrat=residua^2
y_hat=regrese1$fitted.values
y_hat2=(regrese1$fitted.values)^2

white=lm(residua_kvadrat~y_hat+y_hat2)
summary(white)
## 
## Call:
## lm(formula = residua_kvadrat ~ y_hat + y_hat2)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0060495 -0.0020436 -0.0006689  0.0010171  0.0161182 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.375984   0.076635   4.906 1.07e-05 ***
## y_hat       0.088626   0.018001   4.924 1.01e-05 ***
## y_hat2      0.005208   0.001047   4.974 8.45e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003912 on 49 degrees of freedom
## Multiple R-squared:  0.3526, Adjusted R-squared:  0.3261 
## F-statistic: 13.34 on 2 and 49 DF,  p-value: 2.368e-05

Mohli bychom si sestavit celý White test

PG2=PG^2
Pnc2=Pnc^2
PG_Pnc=PG*Pnc
white2=lm(residua_kvadrat~PG+Pnc+PG_Pnc+PG2+Pnc2)

Autokorelace

Otestujeme residua na prítomnost autokorelace. 
Nejprve si zkusíme sami "naprogramovat" test, abychom pochopili princip.
Budeme testovat sdruzenou signifikantnost parametru následujícího modelu.

\(\epsilon_t=\rho_1\epsilon_{t-1}+\rho_2\epsilon_{t-2}+\rho_3\epsilon_{t-3}+v_t\)

kdy \(v_t \sim iid(0,\sigma^2_vI)\)

Pomocná funkce pro tvorbu zpozdených promenných.

lagpad <- function(x, k) {            # x predstavuje casovou radu a k je delka zpozdení
    c(rep(NA, k), x)[1 : length(x)] 
}
residua=regrese1$residuals     # ulozíme si residua z regrese1 
plot(residua,type = "l")     

autokorelace=lm(residua~lagpad(residua, 1)+lagpad(residua, 2)+lagpad(residua, 3))
summary(autokorelace)
## 
## Call:
## lm(formula = residua ~ lagpad(residua, 1) + lagpad(residua, 2) + 
##     lagpad(residua, 3))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.049461 -0.013124  0.005255  0.011004  0.042396 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.0007753  0.0029430   0.263    0.793    
## lagpad(residua, 1)  1.3033641  0.1510804   8.627 4.29e-11 ***
## lagpad(residua, 2) -0.3575265  0.2381365  -1.501    0.140    
## lagpad(residua, 3) -0.1499035  0.1422896  -1.054    0.298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02024 on 45 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.8202, Adjusted R-squared:  0.8082 
## F-statistic: 68.42 on 3 and 45 DF,  p-value: < 2.2e-16

Odhady parametru \(\rho\) jsou lagpad. Z odhadu muzeme usuzovat, ze náhodná slozka se pravdepdobne rídí AR(1) procesem. Signifikantní podle t-testu je pouze první odhad. Nás vsak u techto testu zajímá sdruzená signifikantnost. p-value=0.00044 coz je mensí nez \(\alpha=0,01\) Muzeme tak prijmout hypotézu, ze je prítomna autokorelace.

Nyní jiz víme princip a muzeme vyuzít nekterý z testu pro autokorelaci (TOTO NENÍ NÁPLNÍ KURZU)

test_autokorelace=Box.test(residua,lag = 10)
Jelikoz je prítomna autokorelace, pokusíme se model upravit tak, abychom se jí zbavili.
Vzdy kdyz narazíta na heteroskedasticitu, nebo autokorelaci, zkuste se zamyslet, cím by mohla být zpusobena. Muze se jednat o spatne zvolený funkcní vztah (napr. mzda je vzdy v logaritmu) a nebo doslo k vynechání nekteré dulezité promenné.

Pokud se Vám nepodarí zjistit prícinu, az poté pouzíte robustní odhad chyb prípadne v kombinaci s FGLS.

My si zkusíme do modelu pridat zpozdenou závisle promennou.

Oil$GP1<-(GP1=lagpad(GP, 1))  # pridáme promennou do jiz existujícího data.frame
regrese2=lm(GP~IP+LPg+LPnc+LPuc+t+GP1)
summary(regrese2)
## 
## Call:
## lm(formula = GP ~ IP + LPg + LPnc + LPuc + t + GP1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.08669 -0.02649  0.00924  0.02707  0.09310 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.162831   0.802385  -5.188 5.17e-06 ***
## IP           0.952165   0.192605   4.944 1.16e-05 ***
## LPg          0.682171   0.069704   9.787 1.29e-12 ***
## LPnc        -0.253874   0.156993  -1.617   0.1130    
## LPuc        -0.076707   0.102753  -0.747   0.4593    
## t            0.006590   0.002479   2.658   0.0109 *  
## GP1          0.365510   0.076543   4.775 2.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04359 on 44 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.9978, Adjusted R-squared:  0.9975 
## F-statistic:  3353 on 6 and 44 DF,  p-value: < 2.2e-16
bptest(regrese2)
## 
##  studentized Breusch-Pagan test
## 
## data:  regrese2
## BP = 8.868, df = 6, p-value = 0.1811
residua2=regrese2$residuals
plot(residua2,type = "l")

test_autokorelace=Box.test(residua2,lag = 10)
test_autokorelace
## 
##  Box-Pierce test
## 
## data:  residua2
## X-squared = 27.635, df = 10, p-value = 0.002065

Muzeme videt, ze se nám nepodarilo odstranit autokorelaci. Nyní muzeme pracovat s predpoklladem, ze jsou splneny vsechny GM (pro casové rady máme slabsí predpoklad o exogenite)

Z tohoto duvodu odhadneme smerodatné chyby pomocí robustní metody HAC. (Newey-West estimator)

library(sandwich)
## Warning: package 'sandwich' was built under R version 3.2.5
coeftest(regrese1, vcov = vcovHAC)
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) -7.3452240  0.6259277 -11.7349 1.978e-15 ***
## IP           1.2746479  0.3207395   3.9741 0.0002471 ***
## LPg          0.9732176  0.0657148  14.8097 < 2.2e-16 ***
## LPnc        -0.4489732  0.1070268  -4.1950 0.0001231 ***
## LPuc         0.0599014  0.1155404   0.5184 0.6066321    
## t            0.0106140  0.0080422   1.3198 0.1934354    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Budeme chtít otestovat zdruzenou signifikantnost LPuc a t.

\(H_0: \beta_4=\beta_5=0\)

\(H_1: non \ H_0\)

Prestoze nejsou splneny GM, ukázeme si jak pouzít robustní F-test.

regrese1=lm(GP~IP+LPg+LPnc+LPuc+t)
regrese_omezena=lm(GP~IP+LPg+LPnc)
waldtest(regrese1,regrese_omezena,vcov = vcovHAC)  # robustní verze F-testu
## Wald test
## 
## Model 1: GP ~ IP + LPg + LPnc + LPuc + t
## Model 2: GP ~ IP + LPg + LPnc
##   Res.Df Df      F Pr(>F)
## 1     46                 
## 2     48 -2 1.7401 0.1868

Z p-value=0.1868 vidíme, ze nemuzeme zamítnosut nulovou hypotézu, takze se nám nepodarilo prokázat, ze by \(LPuc\) a \(t\) mely sdruzene signifikantní vliv na GP. \(\textbf{Jedná se o hypotetický test, casová promenná t, muze být klícovou. Více v casových radách ohledne stacionarity}\)