Question 1.Obtain monthly time series for Crude Oil Prices: West Texas Intermediate (WTI), FRED/MCOILWTICO, and the monthly time series for US Regular Conventional Gas Price,FRED/GASREGCOVM. 2.Create a single time series plot with two log prices log for the sample 1995M1-2017M4. 3.Perform unit root tests for both logs. 4.Determine the number of lags to include in cointegration analysis using Schwarz information criterion. 5. estimate a bivariate VEC model for the variables log. 6.Are the adjustment parameters α1 and α2 in the estimated VEC model statistically significant? Are their signs consistent with error correction mechanism that moves the system back to the long run equilibrium, whenever there is a disruption ? 7.eestimate the VEC model with a restriction α2= 0. 8.What is the intuition for imposing the restriction in (7)? 9.Create and plot sequence of one month ahead forecasts (using either fixed, rolling or recursive scheme) for the period 2011M1-2017M4. Obtain the RMSE for the forecast of the gas price,p GAS t. 10.Plot the ACF and PACF . 11.Create a sequence of one month ahead forecast. 12.Compare the RMSE of the forecast. 13.Based on your results in (12), is there any advantage to using the VEC model to forecast the price of the gas, instead of using a simple univariate AR/MA/ARMA model?

library(Quandl)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(zoo)
library(xts)
library(dygraphs)
library(knitr)
library(forecast)
library(urca)
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: lmtest

(1)

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

(2)

ly1=100*log(COP)
ly2=100*log(RCGP)
ly1 <- window(ly1, start=c(1995,1), end=c(2017,4))
## Warning in window.default(x, ...): 'end' value not changed
ly2 <- window(ly2, start=c(1995,1), end=c(2017,4))
plot(cbind(ly1,ly2), col=c(2,4), plot.type = "single", xlab="Years 1995-2017")
legend("right", c(expression(Log_Crude_Oil_Prices),expression(Log_Regular_Conventional_Gas_Price)), col=c(2,4), lty=1, bty = "n")

(3)

dl.y1 <- diff(ly1,lag=1)
dl.y2 <- diff(ly2,lag=1)
par(mfrow=c(2,1))
plot(dl.y1, xlab=" Years 1995-2017", ylab="", main="Diff Log Crude Oil Prices")
plot(dl.y2, xlab=" Years 1995-2017", ylab="", main="Diff Log Regular Conventional Gas Price")

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

All tests (ADF, KPSS, ERS) show that y1 is not stationary.

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

All tests Y1 is stationary.

adf.test(ly2)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ly2
## Dickey-Fuller = -1.6181, Lag order = 6, p-value = 0.7368
## alternative hypothesis: stationary
kpss.test(ly2, null="Trend")
## Warning in kpss.test(ly2, null = "Trend"): p-value smaller than printed p-
## value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  ly2
## KPSS Trend = 0.7202, Truncation lag parameter = 3, p-value = 0.01
y2.urkpss <- ur.kpss(ly2, type="tau", lags="short")
summary(y2.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
y2.ers1 <- ur.ers(ly2, type="P-test", model="trend")
summary(y2.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

All three performed tests (ADF, KPSS, ERS) sugeest that y2 is not stationary.

adf.test(dl.y2)
## Warning in adf.test(dl.y2): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dl.y2
## Dickey-Fuller = -7.9022, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(dl.y2, null="Trend")
## Warning in kpss.test(dl.y2, null = "Trend"): p-value greater than printed
## p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  dl.y2
## KPSS Trend = 0.037877, Truncation lag parameter = 3, p-value = 0.1
dl.y2.urkpss <- ur.kpss(dl.y2, type="tau", lags="short")
summary(dl.y2.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
dl.y2.ers1 <- ur.ers(dl.y2, type="P-test", model="trend")
summary(dl.y2.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

first difference logarithmic of y2 is stationary.

(4)

y <- cbind(ly1, ly2)
y <- window(y, start=c(1995,1), end=c(2010,12))
VARselect(y, type="const", lag.max=12)
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      6      2      6 
## 
## $criteria
##                  1           2           3           4           5
## AIC(n)    7.353484    7.050363    7.000602    7.015203    6.968809
## HQ(n)     7.396638    7.122286    7.101294    7.144664    7.127039
## SC(n)     7.459916    7.227750    7.248944    7.334499    7.359059
## FPE(n) 1561.637976 1153.310834 1097.380193 1113.619041 1063.280429
##                  6           7           8           9          10
## AIC(n)    6.908009    6.931080    6.921235    6.926159    6.945912
## HQ(n)     7.095008    7.146848    7.165772    7.199464    7.247987
## SC(n)     7.369214    7.463239    7.524349    7.600227    7.690935
## FPE(n) 1000.758328 1024.392235 1014.716198 1020.183056 1041.114353
##                 11          12
## AIC(n)    6.956168    6.983254
## HQ(n)     7.287012    7.342867
## SC(n)     7.772146    7.870186
## FPE(n) 1052.556306 1082.326306

Schwarz information criterion suggests using 2 lags.

y.CA <- ca.jo(y, ecdet="const", type="trace", K=2, spec="transitory", season=12)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1]  1.980704e-01  9.957529e-03 -3.903128e-18
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  1.90  7.52  9.24 12.97
## r = 0  | 43.84 17.85 19.96 24.60
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               ly1.l1       ly2.l1     constant
## ly1.l1      1.000000    1.0000000    1.0000000
## ly2.l1     -1.587125   -0.1377993    0.2504144
## constant -275.337070 -429.0868103 -322.0681760
## 
## Weights W:
## (This is the loading matrix)
## 
##           ly1.l1       ly2.l1     constant
## ly1.d 0.04824429 -0.008302568 1.915537e-16
## ly2.d 0.24802732 -0.003053918 9.713652e-16
y.CA <- ca.jo(y, ecdet="const", type="eigen", K=2, spec="transitory", season=12)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1]  1.980704e-01  9.957529e-03 -3.903128e-18
## 
## Values of teststatistic and critical values of test:
## 
##           test 10pct  5pct  1pct
## r <= 1 |  1.90  7.52  9.24 12.97
## r = 0  | 41.94 13.75 15.67 20.20
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               ly1.l1       ly2.l1     constant
## ly1.l1      1.000000    1.0000000    1.0000000
## ly2.l1     -1.587125   -0.1377993    0.2504144
## constant -275.337070 -429.0868103 -322.0681760
## 
## Weights W:
## (This is the loading matrix)
## 
##           ly1.l1       ly2.l1     constant
## ly1.d 0.04824429 -0.008302568 1.915537e-16
## ly2.d 0.24802732 -0.003053918 9.713652e-16

(5)

y.VEC <- cajorls(y.CA, r=1)
y.VEC
## $rlm
## 
## Call:
## lm(formula = substitute(form1), data = data.mat)
## 
## Coefficients:
##          ly1.d     ly2.d   
## ect1      0.04824   0.24803
## sd1       6.08302   4.57391
## sd2       2.31774   2.22917
## sd3       8.03228   6.21726
## sd4       4.09514   5.92019
## sd5       5.27266   6.26258
## sd6       4.71141   4.85558
## sd7       4.94215   3.27466
## sd8       5.02513   4.23428
## sd9       3.40821   4.01838
## sd10      1.76899  -0.01183
## sd11     -0.20628  -0.33635
## ly1.dl1   0.27002   0.16276
## ly2.dl1  -0.11007   0.18754
## 
## 
## $beta
##                 ect1
## ly1.l1      1.000000
## ly2.l1     -1.587125
## constant -275.337070
summary(y.VEC$rlm)
## Response ly1.d :
## 
## Call:
## lm(formula = ly1.d ~ ect1 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + 
##     sd7 + sd8 + sd9 + sd10 + sd11 + ly1.dl1 + ly2.dl1 - 1, data = data.mat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.112  -4.685   1.140   6.214  19.117 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)   
## ect1     0.04824    0.08132   0.593  0.55375   
## sd1      6.08302    3.00078   2.027  0.04416 * 
## sd2      2.31774    3.05064   0.760  0.44842   
## sd3      8.03228    2.99738   2.680  0.00807 **
## sd4      4.09514    3.07874   1.330  0.18520   
## sd5      5.27266    3.14608   1.676  0.09552 . 
## sd6      4.71141    3.17557   1.484  0.13969   
## sd7      4.94215    3.14817   1.570  0.11825   
## sd8      5.02513    3.07894   1.632  0.10445   
## sd9      3.40821    3.04442   1.119  0.26446   
## sd10     1.76899    3.03716   0.582  0.56101   
## sd11    -0.20628    2.96477  -0.070  0.94461   
## ly1.dl1  0.27002    0.11374   2.374  0.01867 * 
## ly2.dl1 -0.11007    0.14766  -0.745  0.45702   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.325 on 176 degrees of freedom
## Multiple R-squared:  0.134,  Adjusted R-squared:  0.06508 
## F-statistic: 1.945 on 14 and 176 DF,  p-value: 0.02469
## 
## 
## Response ly2.d :
## 
## Call:
## lm(formula = ly2.d ~ ect1 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + 
##     sd7 + sd8 + sd9 + sd10 + sd11 + ly1.dl1 + ly2.dl1 - 1, data = data.mat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.3100  -2.3394   0.3539   3.5036  11.8312 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## ect1     0.24803    0.04799   5.169 6.35e-07 ***
## sd1      4.57391    1.77075   2.583 0.010606 *  
## sd2      2.22917    1.80017   1.238 0.217249    
## sd3      6.21726    1.76874   3.515 0.000559 ***
## sd4      5.92019    1.81675   3.259 0.001343 ** 
## sd5      6.26258    1.85649   3.373 0.000913 ***
## sd6      4.85558    1.87389   2.591 0.010367 *  
## sd7      3.27466    1.85772   1.763 0.079681 .  
## sd8      4.23428    1.81687   2.331 0.020912 *  
## sd9      4.01838    1.79650   2.237 0.026556 *  
## sd10    -0.01183    1.79221  -0.007 0.994740    
## sd11    -0.33635    1.74950  -0.192 0.847763    
## ly1.dl1  0.16276    0.06712   2.425 0.016318 *  
## ly2.dl1  0.18754    0.08713   2.152 0.032732 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.913 on 176 degrees of freedom
## Multiple R-squared:  0.4769, Adjusted R-squared:  0.4353 
## F-statistic: 11.46 on 14 and 176 DF,  p-value: < 2.2e-16

(6)

adjustment parameters are estimated as: α1 is not significant and α2 is significant at the level of 1%.and the signs are not consistent with error correction mechanism.

(7)

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.082 0.000 0.000
## 
## The value of the likelihood ratio test statistic:
## 25.68 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.ly1.l1      1.0000
## RK.ly2.l1     -1.5589
## RK.constant -278.3270
## 
## Weights W of the restricted VAR:
## 
##         [,1]
## [1,] -0.2446
## [2,]  0.0000

VEC with restriction on alpha_2

y.VEC.rest1 <- cajorls(y.CA.ralpha, r=1)
y.VEC.rest1
## $rlm
## 
## Call:
## lm(formula = substitute(form1), data = data.mat)
## 
## Coefficients:
##          ly1.d     ly2.d   
## ect1      0.01989   0.22771
## sd1       6.16105   4.65599
## sd2       2.35313   2.25799
## sd3       8.05228   6.27209
## sd4       4.01057   5.86157
## sd5       5.08712   6.11692
## sd6       4.39434   4.58973
## sd7       4.58845   2.97139
## sd8       4.73235   3.97373
## sd9       3.18491   3.80830
## sd10      1.57007  -0.20095
## sd11     -0.32171  -0.44939
## ly1.dl1   0.29063   0.18202
## ly2.dl1  -0.12807   0.17142
## 
## 
## $beta
##                ect1
## ly1.l1      1.00000
## ly2.l1     -1.55894
## constant -278.32704
summary(y.VEC.rest1$rlm)
## Response ly1.d :
## 
## Call:
## lm(formula = ly1.d ~ ect1 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + 
##     sd7 + sd8 + sd9 + sd10 + sd11 + ly1.dl1 + ly2.dl1 - 1, data = data.mat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.946  -4.719   1.290   6.271  18.953 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)   
## ect1     0.01989    0.08040   0.247  0.80487   
## sd1      6.16105    3.00238   2.052  0.04165 * 
## sd2      2.35313    3.05307   0.771  0.44189   
## sd3      8.05228    2.99965   2.684  0.00796 **
## sd4      4.01057    3.08111   1.302  0.19473   
## sd5      5.08712    3.14682   1.617  0.10776   
## sd6      4.39434    3.17109   1.386  0.16758   
## sd7      4.58845    3.14109   1.461  0.14586   
## sd8      4.73235    3.07381   1.540  0.12546   
## sd9      3.18491    3.04160   1.047  0.29648   
## sd10     1.57007    3.03530   0.517  0.60562   
## sd11    -0.32171    2.96558  -0.108  0.91374   
## ly1.dl1  0.29063    0.11263   2.580  0.01069 * 
## ly2.dl1 -0.12807    0.14716  -0.870  0.38534   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.332 on 176 degrees of freedom
## Multiple R-squared:  0.1325, Adjusted R-squared:  0.06354 
## F-statistic: 1.921 on 14 and 176 DF,  p-value: 0.02696
## 
## 
## Response ly2.d :
## 
## Call:
## lm(formula = ly2.d ~ ect1 + sd1 + sd2 + sd3 + sd4 + sd5 + sd6 + 
##     sd7 + sd8 + sd9 + sd10 + sd11 + ly1.dl1 + ly2.dl1 - 1, data = data.mat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.4165  -1.8241   0.8595   3.8376  11.9833 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## ect1     0.22771    0.04790   4.754 4.12e-06 ***
## sd1      4.65599    1.78847   2.603 0.010020 *  
## sd2      2.25799    1.81866   1.242 0.216049    
## sd3      6.27209    1.78684   3.510 0.000569 ***
## sd4      5.86157    1.83537   3.194 0.001665 ** 
## sd5      6.11692    1.87451   3.263 0.001323 ** 
## sd6      4.58973    1.88896   2.430 0.016114 *  
## sd7      2.97139    1.87110   1.588 0.114070    
## sd8      3.97373    1.83102   2.170 0.031329 *  
## sd9      3.80830    1.81183   2.102 0.036984 *  
## sd10    -0.20095    1.80807  -0.111 0.911631    
## sd11    -0.44939    1.76655  -0.254 0.799493    
## ly1.dl1  0.18202    0.06709   2.713 0.007331 ** 
## ly2.dl1  0.17142    0.08766   1.955 0.052119 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.963 on 176 degrees of freedom
## Multiple R-squared:  0.4661, Adjusted R-squared:  0.4236 
## F-statistic: 10.97 on 14 and 176 DF,  p-value: < 2.2e-16

(8)

It is not possible to achieve the long run equilibrium, if we do not restrict one of the alphas to be equal to 0.

(9)

Rolling Scheme one-month ahead forecast

ly1=100*log(COP)
ly2=100*log(RCGP)
y.fcst <- cbind(ly1, ly2)
y.fcst <- na.trim(y.fcst)
y.fcst <- window(y.fcst, start=1995+0, end=2017+3/12)
## Warning in window.default(x, ...): 'end' value not changed
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="blue")

VEC.f.rol <- as.ts(y.pgas2)
accuracy(VEC.f.rol, y.VAR.f2[,1])
##                ME     RMSE      MAE       MPE     MAPE       ACF1
## Test set -1.41565 4.930924 3.547289 -1.667412 3.922441 0.05211672
##          Theil's U
## Test set 0.5858428

RMSE is 4.930924.

(10)

par(mfrow=c(1,2))
acf(dl.y2, type="correlation", lag=48, xlab="Lag", ylab="ACF", main="Gas Price")
acf(dl.y2, type="partial", lag=48, xlab="Lag", ylab="PACF", main="Gas Price")

results of ACF and PACF, we can suggest different ARMA models.

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

#(11)

rol.f <-zoo()
for(i in 1: length(y.2[,2]))
{
  y <- window(dl.y2, end = lstM+(i-1)/12)
  rol.updt <- arima(y, order=c(2,0,2))
  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(2010, 2017), main="One month ahead forecast for ARMA (2,2)")
lines(rol.f, type="p", pch=16, lty="solid", col="blue")
lines(dl.y2, type="l", pch=16, lty="solid")
lines(ly1, type="o", pch=16, lty="dashed")

accuracy(rol.f, dl.y2)
##                  ME     RMSE      MAE MPE MAPE        ACF1 Theil's U
## Test set -0.5683734 5.221744 3.814195 Inf  Inf -0.06078359       NaN

(12,13)

VEC model preforms a better forcast.