Homework 7
Abeer Althawad
May 7, 2017
library("Quandl")
library("urca")
library("vars")
library("forecast")
library("zoo")
library(readr)
Quandl.api_key("sfE3eBMHtyGwpQTDwkWQ")Import data
o <-Quandl("FRED/MCOILWTICO", type="zoo")
g <-Quandl("FRED/GASREGCOVM",type="zoo")a
lo <- log(o)
lg <- log(g)
lo <- window(lo, start=1995+0, end=2017+3/12)
lg <- window(lg, start=1995+0, end=2017+3/12)par(mfrow=c(1,1))
plot(lo, type='l', main="Crude Oil Prices vs. Regular Conventional Gas Price", col="blue", ylim=c(-1, 5))
lines(lg, col="red")
b
lo.urers1 <- ur.ers(lo, type="P-test", model="trend")
summary(lo.urers1)##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 9.6438
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 3.96 5.62 6.89lg.urers2 <- ur.ers(lg, type="P-test", model="trend")
summary(lg.urers2)##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 8.0956
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 3.96 5.62 6.89We reject the null of no unit root.
first difference
do <- diff(lo)
dg <- diff(lg)do.urers1 <- ur.ers(do, type="P-test", model="trend")
summary(do.urers1)##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 0.8329
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 3.96 5.62 6.89dg.urers2 <- ur.ers(dg, type="P-test", model="trend")
summary(dg.urers2)##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type P-test
## detrending of series with intercept and trend
##
## Value of test-statistic is: 0.4831
##
## Critical values of P-test are:
## 1pct 5pct 10pct
## critical values 3.96 5.62 6.89We fail to reject the null of no unit root. Thus, Crude Oil Prices and Regular Conventional Gas Price are time series at I(1).
par(mfrow=c(1,1))
plot(do, type='l', main="log difference of Crude Oil Prices vs. log difference Regular Conventional Gas Price", col="blue", ylim=c(-0.4, 0.4))
lines(dg, col="red")
c
y <- cbind(lo, lg)
colnames(y) <- c("log.wti","log.gas")
y <- na.trim(y)y.VAR.IC <- VARselect(y, type="const")
nlags <- y.VAR.IC$selection["SC(n)"]
nlags## SC(n)
## 3y <- window(y, start=1995+0, end=2010+11/12)y.CA <- ca.jo(y, ecdet="const", type="trace", K=nlags, spec="transitory")
summary(y.CA)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 1.624149e-01 1.223479e-02 2.818926e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 2.33 7.52 9.24 12.97
## r = 0 | 35.82 17.85 19.96 24.60
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## log.wti.l1 log.gas.l1 constant
## log.wti.l1 1.000000 1.00000000 1.0000000
## log.gas.l1 -1.584166 -0.06721038 0.5915855
## constant -2.754478 -3.99184158 -2.8496192
##
## Weights W:
## (This is the loading matrix)
##
## log.wti.l1 log.gas.l1 constant
## log.wti.d 0.0183103 -0.012570961 9.267236e-18
## log.gas.d 0.2109392 -0.005473235 -1.263507e-16y.CA <- ca.jo(y, ecdet="const", type="eigen", K=nlags, spec="transitory")
summary(y.CA)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 1.624149e-01 1.223479e-02 2.818926e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 2.33 7.52 9.24 12.97
## r = 0 | 33.50 13.75 15.67 20.20
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## log.wti.l1 log.gas.l1 constant
## log.wti.l1 1.000000 1.00000000 1.0000000
## log.gas.l1 -1.584166 -0.06721038 0.5915855
## constant -2.754478 -3.99184158 -2.8496192
##
## Weights W:
## (This is the loading matrix)
##
## log.wti.l1 log.gas.l1 constant
## log.wti.d 0.0183103 -0.012570961 9.267236e-18
## log.gas.d 0.2109392 -0.005473235 -1.263507e-16We reject the null hypothesis. Thus, oil and gas prices are cointegrated.
lttest(y.CA, r=1)## LR-test for no linear trend
##
## H0: H*2(r<=1)
## H1: H2(r<=1)
##
## Test statistic is distributed as chi-square
## with 1 degress of freedom
## test statistic p-value
## LR test 0.91 0.34The test suggests case 2(restricted constant) while the plots from (a) suggested Case 4(restricted trend).
d
y.VEC <- cajorls(y.CA, r=1)
y.VEC## $rlm
##
## Call:
## lm(formula = substitute(form1), data = data.mat)
##
## Coefficients:
## log.wti.d log.gas.d
## ect1 0.01831 0.21094
## log.wti.dl1 0.27532 0.13558
## log.gas.dl1 -0.13892 0.35368
## log.wti.dl2 0.17025 0.01702
## log.gas.dl2 -0.05737 -0.14318
##
##
## $beta
## ect1
## log.wti.l1 1.000000
## log.gas.l1 -1.584166
## constant -2.754478β = 1.000000 and α = 0.01831 -1.584166 0.21094 -2.754478
e
summary(y.VEC$rlm)## Response log.wti.d :
##
## Call:
## lm(formula = log.wti.d ~ ect1 + log.wti.dl1 + log.gas.dl1 + log.wti.dl2 +
## log.gas.dl2 - 1, data = data.mat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26410 -0.04813 0.01579 0.06504 0.23386
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ect1 0.01831 0.07813 0.234 0.8150
## log.wti.dl1 0.27532 0.10870 2.533 0.0121 *
## log.gas.dl1 -0.13892 0.15197 -0.914 0.3619
## log.wti.dl2 0.17025 0.10848 1.569 0.1183
## log.gas.dl2 -0.05737 0.14514 -0.395 0.6931
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08399 on 184 degrees of freedom
## Multiple R-squared: 0.07845, Adjusted R-squared: 0.0534
## F-statistic: 3.133 on 5 and 184 DF, p-value: 0.009728
##
##
## Response log.gas.d :
##
## Call:
## lm(formula = log.gas.d ~ ect1 + log.wti.dl1 + log.gas.dl1 + log.wti.dl2 +
## log.gas.dl2 - 1, data = data.mat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.210243 -0.024651 0.003275 0.032126 0.129029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ect1 0.21094 0.04902 4.304 2.73e-05 ***
## log.wti.dl1 0.13558 0.06819 1.988 0.048261 *
## log.gas.dl1 0.35368 0.09534 3.710 0.000275 ***
## log.wti.dl2 0.01702 0.06806 0.250 0.802782
## log.gas.dl2 -0.14318 0.09105 -1.572 0.117572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05269 on 184 degrees of freedom
## Multiple R-squared: 0.3709, Adjusted R-squared: 0.3538
## F-statistic: 21.69 on 5 and 184 DF, p-value: < 2.2e-16α1= 0.01831 is not significant, but α2= 0.07813 is significant. Also,as α1 >0 and α2 >0, that does not satisfy the condition for long run stable relationship.
f
rest.alpha <- matrix(c(1,0), c(2,1))
y.CA.ralpha <- alrtest(y.CA, A=rest.alpha, r=1)
summary(y.CA.ralpha)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Estimation and testing under linear restrictions on beta
##
## The VECM has been estimated subject to:
## beta=H*phi and/or alpha=A*psi
##
## [,1]
## [1,] 1
## [2,] 0
##
## Eigenvalues of restricted VAR (lambda):
## [1] 0.084 0.000 0.000
##
## The value of the likelihood ratio test statistic:
## 16.91 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0
##
## Eigenvectors, normalised to first column
## of the restricted VAR:
##
## [,1]
## RK.log.wti.l1 1.0000
## RK.log.gas.l1 -1.5423
## RK.constant -2.7886
##
## Weights W of the restricted VAR:
##
## [,1]
## [1,] -0.2291
## [2,] 0.0000g
If α1 >0 and α2 >0, y1 and y2 will be growing. However, if α2=0 that means the adjustment occurs by y1(oil price)
h
ly1=100*log(o)
ly2=100*log(g)
y.fcst <- cbind(ly1, ly2)
y.fcst <- na.trim(y.fcst)
y.fcst <- window(y.fcst, start=1995+0, end=2017+3/12)fstM <- 1995+0
lstM <- 2010+11/12
y.1 <- window(y.fcst, end=lstM)
y.2 <- window(y.fcst, start=lstM+1/12)
y.VAR.f1 <-data.frame()
y.VAR.f2 <-data.frame()
for(i in 1:length(y.2[,2]))
{
y <- window(y.fcst, end = lstM + (i-1)/12)
y.CA <- ca.jo(y, ecdet="const", type="eigen", K=3, spec="transitory")
y.VAR <- vec2var(y.CA, r=1)
y.VAR.updt <- predict(y.VAR, n.ahead=1)
y.VAR.f1 <-rbind(y.VAR.f1, as.ts(y.VAR.updt$fcst$ly1))
y.VAR.f2 <-rbind(y.VAR.f2, as.ts(y.VAR.updt$fcst$ly2))
}
y.VAR.f1 <- ts(y.VAR.f1, start=2011, frequency=12)
y.VAR.f2 <- ts(y.VAR.f2, start=2011, frequency=12)
y.pgas <-as.ts(y.fcst[,2])
y.pgas1 <-as.ts(y.1[,2])
y.pgas2 <-as.ts(y.2[,2])par(mfrow=c(1,1))
plot(y.pgas, type="l", pch=16, lty="solid", main="One month ahead forecast for gas price", xlim=c(2010,2017), ylim=c(45,140))
lines(y.pgas1, type="o", pch=16, lty="solid")
lines(y.VAR.f2[,1], type="o", pch=16, lty="dashed", col="red")
yall <- cbind(lo, lg)
colnames(yall) <- c("log.wti","log.gas")
yall <- na.trim(yall)
yall <- window(yall, start=1995+0, end=2017+3/12)
first.m <- 1995+0
last.m <-2010+11/12
y.p1 <- window(yall, end=last.m)
y.p2 <- window(yall, start=last.m+ 1/12)
y.VAR.f1 <-data.frame()
y.VAR.f2 <-data.frame()y.g <-as.ts(yall[,2])
y.g.p1 <-as.ts(y.p1[,2])
y.g.p2 <-as.ts(y.p2[,2])i
first.m <- 1995+0
last.m <-2010+11/12
dg.p1 <- window(dg, end=last.m)
dg.p2 <- window(dg, start=last.m+1/12)par(mfrow=c(1,2))
Acf(dg.p1, type="correlation", lag=48, main="ACF for Gas")
Acf(dg.p1, type="partial", lag=48, main="PACF for Gas")
arma26 <- Arima(dg.p1, order=c(2, 0, 6))
tsdiag(arma26, gof.lag=36)
j
rol.f <-zoo()
for(i in 1: length(dg.p2))
{
y <- window(dg, end = last.m+(i-1)/12)
rol.updt <- arima(y, order=c(2,0,6))
rol.f <- c(rol.f, forecast(rol.updt, 1)$mean)
}
rol.f <-as.ts(rol.f)
par(mfrow=c(1,1))
plot(rol.f, type="o", pch=16, xlim=c(1995, 2017), ylim=c(-0.4, 0.2), main="One month ahead Forecasts for log-diff for Regular Conventional Gas Price")
lines(rol.f, type="p", pch=16, lty="solid", col="red")
lines(dg, type="l", pch=16, lty="dashed")
lines(dg.p1, type="o", pch=16, lty="solid")
k & l
Since VEC model includes the error correction term and it has the ability to reduce the deviation from the equilibrium, it is better.