The data

Quandl.api_key('CEFP3eWxEJwr_uUP9a2D')
COP <- Quandl("FRED/MCOILWTICO", type="ts") 
RCGP <- Quandl("FRED/GASREGCOVM", type="ts")

Time series of the original data sets over sample period. The red trend line represents oil prices and blue trend line represents gas prices.

lCOP=100*log(COP)
lRCGP=100*log(RCGP)

Time series of the log of the data sets over sample period.The red trend line represents oil prices and blue trend line represents gas prices.

Unlike original data set trend lines the logged of the time series for oil and gas prices move together.

Unit root test

OIL PRICES

dlCOP <- diff(lCOP,lag=1)
dlRCGP <- diff(lRCGP,lag=1)
par(mfrow=c(2,1))
plot(dlCOP)

plot(dlRCGP)

adf.test(lCOP)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lCOP
## Dickey-Fuller = -1.7758, Lag order = 6, p-value = 0.6704
## alternative hypothesis: stationary
kpss.test(lCOP, null="Trend")
## Warning in kpss.test(lCOP, null = "Trend"): p-value smaller than printed p-
## value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  lCOP
## KPSS Trend = 0.82349, Truncation lag parameter = 3, p-value = 0.01
COP.urkpss <- ur.kpss(lCOP, type="tau", lags="short")
summary(COP.urkpss)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.5677 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
COP.ers1 <- ur.ers(lCOP, type="P-test", model="trend")
summary(COP.ers1)
## 
## ############################################### 
## # Elliot, Rothenberg and Stock Unit Root Test # 
## ############################################### 
## 
## Test of type P-test 
## detrending of series with intercept and trend 
## 
## Value of test-statistic is: 10.0096 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89

adf, kpss and ers test shows that we can reject the null hypothesis that the log of oilf pricrs has a unit root and is not stationary.

adf.test(dlCOP)
## Warning in adf.test(dlCOP): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlCOP
## Dickey-Fuller = -6.6161, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(dlCOP, null="Trend")
## Warning in kpss.test(dlCOP, null = "Trend"): p-value greater than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  dlCOP
## KPSS Trend = 0.05063, Truncation lag parameter = 3, p-value = 0.1
dlCOP.urkpss <- ur.kpss(dlCOP, type="tau", lags="short")
summary(dlCOP.urkpss)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.0498 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
dlCOP.ers1 <- ur.ers(dlCOP, type="P-test", model="trend")
summary(dlCOP.ers1)
## 
## ############################################### 
## # 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.5705 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89

These tests confirm that the first difference of log oil prices is stationary.

GAS PRICES

adf.test(lRCGP)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lRCGP
## Dickey-Fuller = -1.6181, Lag order = 6, p-value = 0.7368
## alternative hypothesis: stationary
kpss.test(lRCGP, null="Trend")
## Warning in kpss.test(lRCGP, null = "Trend"): p-value smaller than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  lRCGP
## KPSS Trend = 0.7202, Truncation lag parameter = 3, p-value = 0.01
lRCGP.urkpss <- ur.kpss(lRCGP, type="tau", lags="short")
summary(lRCGP.urkpss)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.5035 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
lRCGP.ers1 <- ur.ers(lRCGP, type="P-test", model="trend")
summary(lRCGP.ers1)
## 
## ############################################### 
## # 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.1303 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89

log of gas prices is found to be not stationary.

adf.test(dlRCGP)
## Warning in adf.test(dlRCGP): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlRCGP
## Dickey-Fuller = -7.9022, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(dlRCGP, null="Trend")
## Warning in kpss.test(dlRCGP, null = "Trend"): p-value greater than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  dlRCGP
## KPSS Trend = 0.037877, Truncation lag parameter = 3, p-value = 0.1
dlRCGP.urkpss <- ur.kpss(dlRCGP, type="tau", lags="short")
summary(dlRCGP.urkpss)
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 5 lags. 
## 
## Value of test-statistic is: 0.0403 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216
dlRCGP.ers1 <- ur.ers(dlRCGP, type="P-test", model="trend")
summary(dlRCGP.ers1)
## 
## ############################################### 
## # 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.0533 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89

Log of gas prices is also stationary.

Hence, these tests overall shows that both log series follow I(1) process and these two time series data set are weakly stationary after taking their first differences.

Schwarz information criterion for number of lags

We need to use Schwars (BIC) information criteria.

x <- cbind(lCOP, lRCGP)
x <- window(x, start = 1995, end= 2010 + 11/12)
x.VAR.IC <- VARselect(x, type = "const")
nlags <- x.VAR.IC$selection["SC(n)"]

Now we need to run both the Johansen’s trace and maximum eigenvalue cointegration tests for log of prices og gas and oil.

x.CA.trace <- ca.jo(x, ecdet="trend", type="trace", K=nlags, spec="transitory")
summary(x.CA.trace)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1]  2.276875e-01  4.474969e-02 -1.762481e-19
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  8.70 10.49 12.25 16.26
## r = 0  | 57.79 22.76 25.32 30.45
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              lCOP.l1   lRCGP.l1   trend.l1
## lCOP.l1   1.00000000  1.0000000  1.0000000
## lRCGP.l1 -1.69305972 -0.4201790 -0.7171124
## trend.l1  0.07012778 -0.7603874  1.7286282
## 
## Weights W:
## (This is the loading matrix)
## 
##            lCOP.l1    lRCGP.l1      trend.l1
## lCOP.d  0.07355671 -0.08891914 -1.042340e-17
## lRCGP.d 0.24206508 -0.03215544 -4.001239e-18

This test suggest that r should be equal to 1.

x.CA.eigen <- ca.jo(x, ecdet="trend", type="eigen", K=nlags, spec="transitory"  )
summary(x.CA.eigen)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1]  2.276875e-01  4.474969e-02 -1.762481e-19
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  8.70 10.49 12.25 16.26
## r = 0  | 49.09 16.85 18.96 23.65
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##              lCOP.l1   lRCGP.l1   trend.l1
## lCOP.l1   1.00000000  1.0000000  1.0000000
## lRCGP.l1 -1.69305972 -0.4201790 -0.7171124
## trend.l1  0.07012778 -0.7603874  1.7286282
## 
## Weights W:
## (This is the loading matrix)
## 
##            lCOP.l1    lRCGP.l1      trend.l1
## lCOP.d  0.07355671 -0.08891914 -1.042340e-17
## lRCGP.d 0.24206508 -0.03215544 -4.001239e-18

We find same result from this test.

bivariate VEC model

x.VEC <- cajorls(x.CA.eigen, r=1)
summary(x.VEC$rlm)
## Response lCOP.d :
## 
## Call:
## lm(formula = lCOP.d ~ ect1 + constant + lCOP.dl1 + lRCGP.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.518  -5.350   0.407   5.837  21.570 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## ect1        0.073557   0.064048   1.148   0.2523  
## constant  -19.720485  17.736872  -1.112   0.2676  
## lCOP.dl1    0.212477   0.106215   2.000   0.0469 *
## lRCGP.dl1   0.005519   0.135883   0.041   0.9676  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.38 on 186 degrees of freedom
## Multiple R-squared:  0.07275,    Adjusted R-squared:  0.05281 
## F-statistic: 3.648 on 4 and 186 DF,  p-value: 0.00691
## 
## 
## Response lRCGP.d :
## 
## Call:
## lm(formula = lRCGP.d ~ ect1 + constant + lCOP.dl1 + lRCGP.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.2101  -2.5861   0.3027   3.0801  12.6767 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## ect1        0.24207    0.03990   6.067 7.16e-09 ***
## constant  -66.76909   11.04993  -6.042 8.11e-09 ***
## lCOP.dl1    0.12761    0.06617   1.928 0.055321 .  
## lRCGP.dl1   0.32742    0.08465   3.868 0.000152 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.22 on 186 degrees of freedom
## Multiple R-squared:  0.3757, Adjusted R-squared:  0.3623 
## F-statistic: 27.99 on 4 and 186 DF,  p-value: < 2.2e-16

e

We see that alpha1 is not significant however alpha2 is statistically significant at the 0.01% level and is positive.

reestimation of the VEC model with restriction

restricted.alpha <- matrix(c(1,0), c(2,1))
x.CA.restricted.alpha <- alrtest(x.CA.eigen, A=restricted.alpha, r=1)
summary(x.CA.restricted.alpha)
## 
## ###################### 
## # 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.1031 0.0000 0.0000
## 
## The value of the likelihood ratio test statistic:
## 28.42 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.lCOP.l1   1.0000
## RK.lRCGP.l1 -1.4117
## RK.trend.l1 -0.1135
## 
## Weights W of the restricted VAR:
## 
##         [,1]
## [1,] -0.2439
## [2,]  0.0000
x.VEC.restricted <- cajorls(x.CA.restricted.alpha, r=1)
x.VEC.restricted
## $rlm
## 
## Call:
## lm(formula = substitute(form1), data = data.mat)
## 
## Coefficients:
##            lCOP.d     lRCGP.d  
## ect1        -0.03144    0.19514
## constant     9.23707  -53.13998
## lCOP.dl1     0.28482    0.17653
## lRCGP.dl1   -0.07055    0.25005
## 
## 
## $beta
##                ect1
## lCOP.l1   1.0000000
## lRCGP.l1 -1.4116528
## trend.l1 -0.1134815
summary(x.VEC.restricted$rlm)
## Response lCOP.d :
## 
## Call:
## lm(formula = lCOP.d ~ ect1 + constant + lCOP.dl1 + lRCGP.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.523  -5.630   0.818   5.999  21.572 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)   
## ect1      -0.03144    0.07126  -0.441  0.65953   
## constant   9.23707   19.49799   0.474  0.63624   
## lCOP.dl1   0.28482    0.10626   2.680  0.00801 **
## lRCGP.dl1 -0.07055    0.13314  -0.530  0.59683   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.405 on 186 degrees of freedom
## Multiple R-squared:  0.06715,    Adjusted R-squared:  0.04709 
## F-statistic: 3.347 on 4 and 186 DF,  p-value: 0.01128
## 
## 
## Response lRCGP.d :
## 
## Call:
## lm(formula = lRCGP.d ~ ect1 + constant + lCOP.dl1 + lRCGP.dl1 - 
##     1, data = data.mat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.9057  -2.6117   0.0566   3.2729  16.0288 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## ect1        0.19514    0.04628   4.216 3.87e-05 ***
## constant  -53.13998   12.66332  -4.196 4.20e-05 ***
## lCOP.dl1    0.17653    0.06901   2.558  0.01132 *  
## lRCGP.dl1   0.25005    0.08647   2.892  0.00429 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.459 on 186 degrees of freedom
## Multiple R-squared:  0.3174, Adjusted R-squared:  0.3028 
## F-statistic: 21.63 on 4 and 186 DF,  p-value: 1.148e-14

g

We should restrict the model i.e. alpha equal to zero otherwise we cannot achieve long run equilibrium.

h

forcasting using by using rolling scheme

VEC.f.rol <- as.ts(x.priceofgas2)
accuracy(VEC.f.rol, x.VAR.f2[,1])
##                 ME     RMSE      MAE       MPE     MAPE       ACF1
## Test set -1.478578 4.949706 3.582169 -1.748946 3.974304 0.05192568
##          Theil's U
## Test set 0.5926838

RMSE of the forcast is equal to 4.94%.

i

The ARMA model

par(mfrow=c(1,2))
acf(dlRCGP, type="correlation")
acf(dlRCGP, type="partial")

Considering the ACF and PACF, different ARMA model is suggested. By chossing ARMA (2,2)

arma22 <- Arima(x.priceofgas1, order=c(2, 0, 2))
tsdiag(arma22, gof.lag=36)

j

FOrcasting by using the Estimated VEC model

accuracy(rol.f, dlRCGP)
##                  ME     RMSE      MAE MPE MAPE        ACF1 Theil's U
## Test set -0.5239191 5.197283 3.800984 Inf  Inf -0.05525264       NaN

RMSE is 5.19%.

Conclusion

VEC model is good fit for forcasting as compare to ARMA model.