Final Assignment - Forecasting MAth1307

Forecasting Crude Oil Prices based on its time-series data from 1986 - 2016

Student Details

• Dilip Chandra (S3574580)

Introduction

Non renewable sources of energy, albeit depleting, are still an integral part of our lives. Crude oil, a form of fossil fuel is used widely around the world. The advantages and use of petroleum vary from being a source of fuel to run planes, in manufacturing of clothes, for insulation and more.

In this project, we aim to understand how crude oil prices have been varying over the years (since 1986 till 2016), and attempt to forecast the next 10 months ahead prices using relevant model fitting and functions on R studio.

The dataset has been sourced from US Energy Information Adminstration website (https://www.eia.gov/dnav/pet/pet_pri_spt_s1_d.htm). The dataset was tidied up on excel and had 2 variables date and crude oil price and had 368 observations ranging from January, 1986 till August, 2016. In order to forecast this we also compared it with a few series to ensure a predictor series for it. These are as follows:

1. Crude Oil Production (From Jan, 1986 till Aug, 2016, sourced from https://www.eia.gov/dnav/pet/pet_crd_crpdn_adc_mbbl_m.htm)
2. Baltic Dirty Tanker Index (From Aug, 2002 till Aug, 2016, sourced from https://au.investing.com/indices/baltic-dirty-tanker-historical-data)
3. S&P 500 Historical Data (From Jan, 1986 till Aug, 2016, sourced from https://fred.stlouisfed.org/series/SP500)

Methodology

To undertake this research, forecasting methods on R Studio are being used to infer from the dataset.

Research and Inferences

Determining the seasonal, trend and remainder effects on the time series and commenting on stationarity and choosing an appropriate predictor series.

We Read in the dataset and convert it into a time series. We plot the series to check for visual inferences in the dataset.

cop <- read_csv("~/Desktop/Forecasting Final PRoject/Crude Oil Prices Dataset.csv")

## Parsed with column specification:
## cols(
##   Date = col_character(),
##   Cush = col_double()
## )

cop = ts(cop$Cush, start= c(1986,1), end = c(2016,8), frequency = 12)
plot(cop, type = "o", ylab = "Price", main = "Fig. 1 - Time series plot for Crude oil price")

From the time series above in Fig 1. we can see that uptil 2002, there has been a lot of variance in the crude oil price, and post 2002 the crude oil price has shot up. There is clear indication of a few interventions between 2007 and 2009 and again in 2014. To get a better understanding of the series, we plot it with the months as pointers.

plot(cop, type = "o", ylab = "Price", main = "Fig. 2 - Time series plot for Crude oil price")
points(y = cop, x = time(cop), pch = as.vector(season(cop)))

From Fig. 2, it is not very clear to determine a seasonal or trend factor in the time series. To get a better understanding, we plot the ACF and PACF charts and run a dickey fuller test to check for stationarity.

par(mfrow=c(1,2))
 acf(cop,  main="ACF")
 pacf(cop,  main="PACF")

 adf.test(cop)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  cop
## Dickey-Fuller = -2.1205, Lag order = 7, p-value = 0.5261
## alternative hypothesis: stationary

With a p value of .5261, clearly the time series is not stationary. The ACF plot doesnt have curves in the lags, so seasonality is not evident in the series. The high 1st lag in PACF also shows evidence of non stationarity in the series.

In order to make the series stationary, we first try a transformation using the BoxCox method to check for the lambda.

lambda = BoxCox.lambda(cop, method = 'loglik')
 lambda

## [1] -0.45

With a lambda value of -0.45, we use inverse power transformation on the series and run the augmented dickey fuller test to check for stationarity again.

cop.tr = (cop^lambda)
adf.test(cop.tr)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  cop.tr
## Dickey-Fuller = -2.2389, Lag order = 7, p-value = 0.4761
## alternative hypothesis: stationary

The p value has reduced slightly but is still at .4761 which is well above the significant value of 0.05.

Since transformation did not work in making the series stationary, we difference the series using ordinary differencing and plot and test for stationarity using the dickey fuller test.

cop.diff = diff(cop)
 plot(cop.diff,type = "o", ylab='Change in price index',xlab='Time', main = "Fig. 3 - Differenced time series plot for Crude oil price")

adf.test(cop.diff)

## Warning in adf.test(cop.diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  cop.diff
## Dickey-Fuller = -7.6526, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

The differenced series in Fig. 3 has a better variance around the mean level, and the peaks are evidence of the interventions in the orignal series. The p value of the dickey fuller test is significant and is now 0.01 and hence the series is now stationary.

We move on to decomposing the time series to remark on the seasonal, trend and remainder effects on the series.

We first run x12 decomposition.

cop.decom.x12 = x12(cop)
 plot(cop.decom.x12 , sa=TRUE , trend=TRUE, main = "Fig. 4 - X12 decomposed time series plot for Crude oil price ")

From Fig. 4 we can see that the seasonally adjusted series is closely following the original series and shows no evidence of a seasonal effect. The trend deviates slightly from the orignal, showing that it has an effect on the series.

To understand the seasonal effect of the series we plot the seasonal factors monthly.

plotSeasFac(cop.decom.x12, main = "Fig. 5 - Seasonal Factors by period and SI Rations for Crude Oil Prices")

We can see that in the months January, June, July and October, the expected pattern deviates from the mean values.

It is also observed from the SI Ratios, there exists influantial observations for all months.

We next undertake STL decomposition to isolate seasonal and trend effects to infer on the trend effects.

cop.decom <- stl(cop, t.window=15, s.window="periodic", robust=TRUE)
plot(cop.decom, main = "Fig. 6 - STL Decomposition of Crude Oil Prices")

From Fig. 6, we can see that the trend follows the overall pattern of the original series with an upward trend from 2003-2004. The remainder part shows the peaks which reflect the interventions in the series.

We break up the trend effect with the monthplot to understand the trends over the month.

monthplot(cop.decom, choice = "trend", main="Fig. 7 - Trend component of Crude Oil Prices", ylab="Trend")

Fig. 7 shows that the mean over the months remains pretty much the same, with values of prices vary a lot between 10 and 100. This clearly shows the trend effect on the original series.

Before we go into model fitting, we found three series that could potentially act as predictor series for forecasting future crude oil prices. The first one we used was crude oil production dataset. We load the dataset and join it with the original series and plot it simulataneously to check for visual inferences to correlation.

crudeprod <- read_csv("~/Desktop/Forecasting Final PRoject/crudeprod.csv")

## Parsed with column specification:
## cols(
##   Date = col_character(),
##   Total = col_double()
## )

crudeprod = ts(crudeprod$Total, start = c(1986,1), end = c(2016,8), frequency = 12)

data.joint = ts.intersect(cop,crudeprod)
colnames(data.joint) =  c("Crude Oil Price","Crude Oil Production")
plot(data.joint , yax.flip=T, main = "Fig. 8 - Timeseries plot of Crude Oil Price and Production")

Fig. 8, shows that both series have an overall upward trend, minus the intervention factors affecting the price.

cor(data.joint)

##                      Crude Oil Price Crude Oil Production
## Crude Oil Price             1.000000             0.731201
## Crude Oil Production        0.731201             1.000000

The series’ have a 0.731 correlation value which is quite high. To check whether or not the series are spuriously correlated we undertake CCF tests.

ccf(as.vector(data.joint[,1]), as.vector(data.joint[,2]),ylab='CCF', main = "Fig. 9 - Sample CCF between Crude Oil Price and Production")

There appears to be a high correlation structure in the CCF plot, and a lot of cross -correlations are significantly different from zero.

We will now take the CCF of the differenced series to see if it has any variation to the previous CCF Plot.

ccf(as.vector(diff(data.joint[,1])), as.vector(diff(data.joint[,2])),ylab='CCF', main = "Fig. 10 - Sample CCF between Differenced Crude Oil Price and Production")

There appear to be some significant correlation in the CCF between the differenced time-series. The number of significant lags have reduced significantly.

This is not enough to say for certain, that there is no spurious correlation between the two series.

We will move on to prewhitening the series.

data.dif=ts.intersect(diff(diff(cop, lag = 12)),diff(diff(crudeprod, lag = 12)))

prewhiten(as.vector(data.dif[,1]), as.vector(data.dif[,2]), ylab='CCF', main="Fig. 11 - Sample CFF after prewhitening")

From Fig. 11, it seems that there is no correlation between the two series after prewhitening.

Therefore, it seems that the two series are definitely correlated, and the strong correlation pattern found between them in the dataset is not spurious.

We next undertake the same steps for S&P 500 Historical data. We load the dataset and join it with the original series and plot it simulataneously to check for visual inferences to correlation.

sp<- read_csv("~/Desktop/Forecasting Final PRoject/S&P Historical Data.csv")

## Parsed with column specification:
## cols(
##   Date = col_character(),
##   `S&P 500 Price` = col_number()
## )

sp = ts(sp$`S&P 500 Price`, start = c(1986,1), end = c(2016,8), frequency = 12)

data2.joint = ts.intersect(cop,sp)
colnames(data2.joint) =  c("Crude Oil Price","S&P 500 Price Index")
plot(data2.joint , yax.flip=T, main = "Fig. 12 - Timeseries plot of Crude Oil Price and S&P 500 Price Index")

Fig. 12, shows that both series have an overall upward trend, minus the intervention factors affecting the price.

cor(data2.joint)

##                     Crude Oil Price S&P 500 Price Index
## Crude Oil Price           1.0000000           0.6344546
## S&P 500 Price Index       0.6344546           1.0000000

The series’ have a 0.634 correlation value which is quite high, but not as high as the correlation with the crude oil production dataset. To check whether or not the series are spuriously correlated we undertake CCF tests.

ccf(as.vector(data2.joint[,1]), as.vector(data2.joint[,2]),ylab='CCF', main = "Fig. 13 - Sample CCF between Crude Oil Price and S&P 500 Price Index")

There appears to be a high correlation structure in the CCF plot, and a lot of cross -correlations are significantly different from zero.

We will now take the CCF of the differenced series to see if it has any variation to the previous CCF Plot.

ccf(as.vector(diff(data2.joint[,1])), as.vector(diff(data2.joint[,2])),ylab='CCF', main = "Fig. 14 - Sample CCF of Differenced Crude Oil Price and S&P 500 Price Index")

There appear to be some significant correlation in the CCF between the differenced time-series. The number of significant lags have reduced significantly.

This is not enough to say for certain, that there is no spurious correlation between the two series.

We will move on to prewhitening the series.

data.dif2=ts.intersect(diff(diff(cop, lag = 12)),diff(diff(sp, lag = 12)))

prewhiten(as.vector(data.dif2[,1]), as.vector(data.dif2[,2]), ylab='CCF', main="Fig. 15 - Sample CFF after prewhitening")

From Fig15. it seems that there is no correlation between crude oil prices and S&P 500 historical price data..

The significant correlations in Fig 13 and 14 can be said to be related to false alarm rate of CCF.

Therefore, it seems that the two series are uncorrelated, and the strong correlation pattern found between them in the dataset is indeed spurious.

We next undertake the same steps with a third potential predictor series, the Baltic Dirty Tanker Price Index. We load the dataset and join it with the original series and plot it simulataneously to check for visual inferences to correlation.

baltic <- read_csv("~/Desktop/Forecasting Final PRoject/forecastingShippingPrices.csv")

## Parsed with column specification:
## cols(
##   date = col_character(),
##   price = col_number()
## )

baltic.ts <- ts(baltic[182:1,2], start = c(2002,8),end = c(2016,8), frequency = 12)

data3.joint = ts.intersect(cop, baltic.ts)
plot(data3.joint, yax.flip=T, main = "Fig. 16 - Timeseries plot of Crude Oil Price and Baltic Dirty Tanker Price")

The above figure shows that there might be an inverse correlation between the two series. We run a correlation test.

oil.ts.baltic <- ts(cop[200:368])
cor(oil.ts.baltic, baltic.ts[1:169])

## [1] -0.1861289

The negative value of -0.186, confirms our inference from the visualisation and shows a very low inverse correlation.

We move on to plot the CCF chart.

ccf(as.vector(data3.joint[,1]), as.vector(data3.joint[,2]),ylab='CCF', main = "Fig. 17 - Sample CCF between Crude Oil Price and Baltic Dirty Shipping Index")

There appears to be a high correlation structure in the CCF plot, and a lot of cross -correlations are significantly different from zero.

We will now take the CCF of the differenced series to see if it has any variation to the previous CCF Plot.

ccf(as.vector(diff(data3.joint[,1])), as.vector(diff(data3.joint[,2])),ylab='CCF', main = "Fig. 18 - Sample CCF of Differenced Crude Oil Price and Baltic Shipping Index")

There appear to be some significant correlation in the CCF between the differenced time-series. The number of significant lags have reduced significantly.

This is not enough to say for certain, that there is no spurious correlation between the two series.

We will move on to prewhitening the series.

shipping.dif=ts.intersect(diff(diff(baltic.ts,12)),diff(diff(cop,12)))

prewhiten(as.vector(shipping.dif[,1]),as.vector(shipping.dif[,2]),ylab='CCF', main = "Fig. 19 - Sample CCF post prewhitening")

From Fig 19. it seems that there is no correlation between residential property price index (PPI) and quarterly poupulation change.

The significant correlations in Fig 17 and 18 can be said to be related to false alarm rate of CCF.

Therefore, it seems that the two series are uncorrelated, and the strong correlation pattern found between them in the dataset is indeed spurious.

We therefore create a predictor series using the 10 months ahead rates of the crude oil production series to model fit and forecast our crude oil price series.

pred <- read_csv("~/Desktop/Forecasting Final PRoject/predictor.csv")

## Parsed with column specification:
## cols(
##   Date = col_character(),
##   Total = col_integer()
## )

pred = ts(pred$Total, start = c(2016,9), frequency=)

Model Fitting, Specification and Forecasting

To find the best model that would fit the crude oil price time series, we undertake a systematic approach of finding models from DLM, Dynlm and ETS models.

We first try to find a good DLM model that fits the data appropriately.

model1p = dlm(x = as.vector(crudeprod) , y = as.vector(cop), q = 1)
vif(model1p$model)

##      x.t      x.1 
## 75.13132 75.13132

We see that as we increase lags from 1 through 12, at q = 1 the R squared value is coming to .5352 and hence its not a very good fit. We also find that the lags are affected by multicollinearity, as values of VIF test are well above 10. although the p value is significant the model doesnt seem very good.

We continue to check residuals.

checkresiduals(model1p$model$residuals)

The residuals seem to be slightly random, but there are many extreme lags in the ACF test, suggesting evidence of serial correlation.

bgtest(model1p$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model1p$model
## LM test = 349.32, df = 1, p-value < 2.2e-16

The Breusch-Godfrey test having a significant p value (< 2.2e-16) shows that the model is significant and fits the data well, although the coefficients are insignificant.

We next try the polydlm fitting.

model2p = polyDlm(x = as.vector(crudeprod) , y = as.vector(cop) , q = 1 , k = 1, show.beta = TRUE)

## Estimates and t-tests for beta coefficients:
##        Estimate Std. Error t value P(>|t|)
## beta.0  0.00328    0.00193   1.700  0.0903
## beta.1  0.00131    0.00192   0.679  0.4980

vif(model2p$model)

##     z.t0     z.t1 
## 300.8932 300.8932

Model and intercept are significant, however all parameters are insignificant, and R squared value of fit is low ie. 0.535. We also find that the lags are affected by multicollinearity, as values of VIF test are well above 10. We continue to check for residuals.

checkresiduals(model2p$model$residuals)

The residuals randomness is still not good, and there are still many extreme lags in the ACF test suggesting serial correlation.

bgtest(model2p$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model2p$model
## LM test = 349.32, df = 1, p-value < 2.2e-16

The test has a significant p value shows that the model is significant and fits the data well, despite the low R squared value.

We next use the ARDLM fitting.

model3p = ardlDlm(x = as.vector(crudeprod) , y = as.vector(cop) , p = 3 , q = 3)

The R squared value for the ardlm is coming to .983 which shows a very good fit for the model. Intercept is insignificant and a few significant y,t lags. We move on to diagnostic testing.

checkresiduals(model3p$model$residuals)

We see that the randsomness is much better than the previous models and even less extreme lags showing low correlation and symmetric histogram.

bgtest(model3p$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model3p$model
## LM test = 0.10834, df = 1, p-value = 0.742

With a p value < 4.136e-15 the model is statistically significant and the best DLM model fit so far.

We also undertake the koyck dlm fit test.

model4p = koyckDlm(x = as.vector(crudeprod) , y = as.vector(cop))
vif(model4p$model)

##      Y.1      X.t 
## 2.208511 2.208511

We see that the Koyck model shows a high R squared value of .9797 and one lag and model are significant, hence it is a good fit. We run diagnostic check.

checkresiduals(model4p$model$residuals)

The randomness is much better than the dlm and polydlm fits, but there are still some extreme values in the ACF plot showing serial correlation.

bgtest(model4p$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model4p$model
## LM test = 61.541, df = 1, p-value = 4.336e-15

With a p value < 4.336e-15 the model is statistically significant.

For dynlm fitting we first start off with SES fitting.

fit1.ses = ses(cop, initial="simple", h=10) 
summary(fit1.ses)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = cop, h = 10, initial = "simple") 
## 
##   Smoothing parameters:
##     alpha = 1 
## 
##   Initial states:
##     l = 22.93 
## 
##   sigma:  4.3152
## Error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.05921196 4.315164 2.725299 -0.2112265 6.579527 0.2447538
##                   ACF1
## Training set 0.3955798
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 2016          44.72 39.18990 50.25010 36.26243 53.17757
## Oct 2016          44.72 36.89925 52.54075 32.75920 56.68080
## Nov 2016          44.72 35.14158 54.29842 30.07107 59.36893
## Dec 2016          44.72 33.65979 55.78021 27.80487 61.63513
## Jan 2017          44.72 32.35431 57.08569 25.80831 63.63169
## Feb 2017          44.72 31.17407 58.26593 24.00328 65.43672
## Mar 2017          44.72 30.08872 59.35128 22.34339 67.09661
## Apr 2017          44.72 29.07850 60.36150 20.79839 68.64161
## May 2017          44.72 28.12969 61.31031 19.34730 70.09270
## Jun 2017          44.72 27.23227 62.20773 17.97483 71.46517

With a simple ses fitting we see that the MASE value is .2447, which is not bad. We do diagnostic check.

checkresiduals(fit1.ses)

## 
##  Ljung-Box test
## 
## data:  Residuals from Simple exponential smoothing
## Q* = 129.29, df = 22, p-value < 2.2e-16
## 
## Model df: 2.   Total lags used: 24

The randomness in the series is evident, and extreme lags show existence of serial correlation.

We next try the holt model.

fit2.holt = holt(cop, initial="simple", h=10)
summary(fit2.holt)

## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = cop, h = 10, initial = "simple") 
## 
##   Smoothing parameters:
##     alpha = 1 
##     beta  = 0.4626 
## 
##   Initial states:
##     l = 22.93 
##     b = -7.47 
## 
##   sigma:  4.3678
## Error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.04336389 4.367767 2.901516 0.5123177 7.19339 0.2605795
##                    ACF1
## Training set 0.09404874
## 
## Forecasts:
##          Point Forecast       Lo 80     Hi 80      Lo 95     Hi 95
## Sep 2016       44.63290  39.0353835  50.23042  36.072236  53.19357
## Oct 2016       44.54580  34.6280179  54.46359  29.377859  59.71375
## Nov 2016       44.45871  29.8127606  59.10465  22.059666  66.85775
## Dec 2016       44.37161  24.5431224  64.20009  14.046557  74.69666
## Jan 2017       44.28451  18.8331689  69.73585   5.360043  83.20897
## Feb 2017       44.19741  12.7057878  75.68903  -3.964870  92.35969
## Mar 2017       44.11031   6.1835864  82.03704 -13.893610 102.11423
## Apr 2017       44.02321  -0.7129715  88.75940 -24.394878 112.44131
## May 2017       43.93612  -7.9657581  95.83799 -35.440951 123.31318
## Jun 2017       43.84902 -15.5587605 103.25680 -47.007339 134.70537

The Holt model gives us a MASE value of .2605, which is higher than the SES model. We move on to diagnostic check.

checkresiduals(fit2.holt)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt's method
## Q* = 66.669, df = 20, p-value = 6.31e-07
## 
## Model df: 4.   Total lags used: 24

The randomness in the data is still evident, and there is still existence of serial correlation.

fit3.hw = hw(cop,seasonal="additive", h=10)
summary(fit3.hw) 

## 
## Forecast method: Holt-Winters' additive method
## 
## Model Information:
## Holt-Winters' additive method 
## 
## Call:
##  hw(y = cop, h = 10, seasonal = "additive") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.0086 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 13.734 
##     b = 0.5659 
##     s = -2.4547 -0.7013 0.3433 1.2946 1.7774 1.4374
##            1.1682 1.2228 0.937 -0.435 -2.126 -2.4636
## 
##   sigma:  4.385
## 
##      AIC     AICc      BIC 
## 3279.764 3281.512 3346.201 
## 
## Error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.2182562 4.288613 2.774115 -1.168417 7.031722 0.2491379
##                   ACF1
## Training set 0.3608986
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 2016       44.11133 38.49173 49.73094 35.51689 52.70577
## Oct 2016       43.03134 35.05001 51.01267 30.82495 55.23773
## Nov 2016       41.85787 32.04074 51.67501 26.84386 56.87189
## Dec 2016       39.97757 28.59301 51.36213 22.56639 57.38876
## Jan 2017       39.84164 27.05874 52.62455 20.29188 59.39141
## Feb 2017       40.04810 25.98524 54.11096 18.54081 61.55539
## Mar 2017       41.61199 26.35758 56.86641 18.28238 64.94161
## Apr 2017       42.85538 26.47836 59.23240 17.80889 67.90187
## May 2017       43.01419 25.57007 60.45831 16.33570 69.69267
## Jun 2017       42.83234 24.36686 61.29782 14.59182 71.07285

The Holt-Winter’s additive method gives us a MASE value of .2526, which is also good, and we move on to diagnostic check

checkresiduals(fit3.hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' additive method
## Q* = 111.35, df = 8, p-value < 2.2e-16
## 
## Model df: 16.   Total lags used: 24

The randomness for the series is very good. Histogram is symmetric. Few extreme lags.

We next check the Holt Winters with an additive seasonal and damped criteria.

fit4.hw = hw(cop,seasonal="additive", damped = TRUE , h=10)
summary(fit4.hw) 

## 
## Forecast method: Damped Holt-Winters' additive method
## 
## Model Information:
## Damped Holt-Winters' additive method 
## 
## Call:
##  hw(y = cop, h = 10, seasonal = "additive", damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.1412 
##     gamma = 1e-04 
##     phi   = 0.8001 
## 
##   Initial states:
##     l = 14.1892 
##     b = -0.064 
##     s = -2.5215 -1.2214 0.583 1.7441 1.5899 2.2624
##            1.0995 0.5277 0.8193 -0.1246 -2.0798 -2.6786
## 
##   sigma:  4.2768
## 
##      AIC     AICc      BIC 
## 3262.334 3264.294 3332.679 
## 
## Error measures:
##                     ME     RMSE      MAE         MPE    MAPE      MASE
## Training set 0.0509217 4.176879 2.824332 -0.06622927 7.32543 0.2536477
##                   ACF1
## Training set 0.2511475
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
## Sep 2016       44.89649 39.41551 50.37747 36.514053 53.27893
## Oct 2016       43.75817 35.55791 51.95843 31.216951 56.29938
## Nov 2016       41.97099 31.44793 52.49405 25.877364 58.06462
## Dec 2016       40.68111 28.04748 53.31475 21.359638 60.00259
## Jan 2017       40.53801 25.94294 55.13307 18.216782 62.85923
## Feb 2017       41.14022 24.70262 57.57782 16.001085 66.27936
## Mar 2017       43.10368 24.92448 61.28289 15.300988 70.90638
## Apr 2017       44.05492 24.22263 63.88721 13.724043 74.38579
## May 2017       43.76924 22.36288 65.17560 11.031032 76.50745
## Jun 2017       44.34109 21.43192 67.25026  9.304532 79.37764

The MASE of the model is 0.2511.

checkresiduals(fit4.hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt-Winters' additive method
## Q* = 81.758, df = 7, p-value = 5.995e-15
## 
## Model df: 17.   Total lags used: 24

The randomness for the series is very good. Histogram is symmetric. Few extreme lags.

We next go for a Holt-Winters model with a multiplicative seasonal component.

fit5.hw = hw(cop,seasonal="multiplicative" , h=10)
summary(fit5.hw)

## 
## Forecast method: Holt-Winters' multiplicative method
## 
## Model Information:
## Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = cop, h = 10, seasonal = "multiplicative") 
## 
##   Smoothing parameters:
##     alpha = 0.9095 
##     beta  = 1e-04 
##     gamma = 0.0902 
## 
##   Initial states:
##     l = 15.4879 
##     b = 0.036 
##     s = 0.92 0.9919 1.0571 0.9215 0.9491 0.8725
##            1.0978 1.0572 1.0455 1.0097 1.058 1.0197
## 
##   sigma:  0.0996
## 
##      AIC     AICc      BIC 
## 3087.834 3089.583 3154.271 
## 
## Error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.02832798 4.485616 2.937047 -0.3523797 7.278162 0.2637704
##                   ACF1
## Training set 0.3735573
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 2016       44.11577 38.48600 49.74555 35.50577 52.72577
## Oct 2016       42.74986 35.36151 50.13820 31.45036 54.04935
## Nov 2016       41.58520 32.90752 50.26289 28.31383 54.85658
## Dec 2016       40.17091 30.55415 49.78768 25.46334 54.87848
## Jan 2017       40.69371 29.83763 51.54980 24.09077 57.29666
## Feb 2017       42.95818 30.42495 55.49141 23.79026 62.12610
## Mar 2017       45.89270 31.44085 60.34456 23.79050 67.99491
## Apr 2017       46.80380 31.04931 62.55830 22.70938 70.89822
## May 2017       47.01407 30.22417 63.80396 21.33614 72.69199
## Jun 2017       46.96427 29.27503 64.65350 19.91092 74.01762

The MASE is good, with a value of 0.261. We go for diagnostic test.

checkresiduals(fit5.hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 58.94, df = 8, p-value = 7.521e-10
## 
## Model df: 16.   Total lags used: 24

The randomness for the series is very good. Histogram is symmetric. Fewer extreme lags than before.

We next undertake the multiplicative seasonal and exponential component in the Holt- Winters’s model.

fit6.hw = hw(cop,seasonal="multiplicative", exponential = TRUE , h=10)
summary(fit6.hw)

## 
## Forecast method: Holt-Winters' multiplicative method with exponential trend
## 
## Model Information:
## Holt-Winters' multiplicative method with exponential trend 
## 
## Call:
##  hw(y = cop, h = 10, seasonal = "multiplicative", exponential = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9846 
##     beta  = 0.0483 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 13.9698 
##     b = 0.608 
##     s = 0.9423 0.971 0.9887 1.0179 1.0271 1.0548
##            1.0402 1.0296 1.027 1 0.9536 0.9478
## 
##   sigma:  0.1575
## 
##      AIC     AICc      BIC 
## 3411.788 3413.537 3478.226 
## 
## Error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.02883452 4.367788 2.897789 1.297859 7.629545 0.2602448
##                  ACF1
## Training set 0.395576
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95     Hi 95
## Sep 2016       43.86149 35.02370 52.86096 30.37046  57.24947
## Oct 2016       42.19781 30.13422 54.98070 24.62330  62.60310
## Nov 2016       41.03853 26.70327 56.98111 21.44719  67.03017
## Dec 2016       39.44646 23.65515 57.93659 18.13187  70.19553
## Jan 2017       39.30831 21.98901 59.94947 16.25777  75.79111
## Feb 2017       39.15776 20.37133 62.89304 14.44239  83.44962
## Mar 2017       40.67086 19.57099 69.79708 13.66712  94.53450
## Apr 2017       41.37649 18.42556 73.63236 12.41565 102.11371
## May 2017       41.08450 17.25634 77.39567 10.98935 110.53490
## Jun 2017       41.10594 16.07972 79.92068 10.22552 119.56235

checkresiduals(fit6.hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method with exponential trend
## Q* = 323.76, df = 8, p-value < 2.2e-16
## 
## Model df: 16.   Total lags used: 24

The MASE is coming at 0.28 which is higher than before and the randomness and extreme lags are still evident.

Based on the MASE value (0.252), extreme lags and residual randomness we choose the Holt-Winter’s seasonal additive model to forecast.

plot(cop, type="o", xlim = c(1986, 2018),  ylab = "Crude Oil Price", xlab = "Year", 
     main="Fig. 21 - Crude Oil Price Forecast")                       
lines(fit3.hw$mean, type="o", col="purple")
lines(fitted(fit3.hw), col="purple")
legend("topleft",lty=1, pch = 1, text.width = 20, col=c("black","purple"), 
       c("Crude Oil Price","Holt-Winter's Additive"))

In order to conduct a systematic analysis, many state space models were applied to the series.

fit1.etsA = ets(cop, model = "ANN")
summary(fit1.etsA)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = cop, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 22.9331 
## 
##   sigma:  4.3271
## 
##      AIC     AICc      BIC 
## 3256.335 3256.401 3268.059 
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.05920942 4.315335 2.725375 -0.2113092 6.579705 0.2447606
##                   ACF1
## Training set 0.3956712

checkresiduals(fit1.etsA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,N)
## Q* = 129.33, df = 22, p-value < 2.2e-16
## 
## Model df: 2.   Total lags used: 24

The first state space model to be explored was the “ANN” model, with additive error and no trend or seasonality. It produced an AIC of 3256, BIC of 3268 and MASE of 0.2447. This turned out to be a rather good model. The residual check showed errors to be seemingly random, with stable variance until around 2005, after which it increased somewhat. The ACF of residuals still had several significant lags, but the histogram appears to indicate a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit2.etsA = ets(cop, model = "AAN")
summary(fit2.etsA)

## ETS(A,Ad,N) 
## 
## Call:
##  ets(y = cop, model = "AAN") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.476 
##     phi   = 0.8 
## 
##   Initial states:
##     l = 17.5462 
##     b = -1.1735 
## 
##   sigma:  4.1206
## 
##      AIC     AICc      BIC 
## 3223.324 3223.556 3246.772 
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.02859116 4.092558 2.753527 0.1190606 6.80118 0.2472889
##                    ACF1
## Training set 0.04340901

checkresiduals(fit2.etsA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,N)
## Q* = 52.997, df = 19, p-value = 4.686e-05
## 
## Model df: 5.   Total lags used: 24

The second state space model to be explored was the “AAN” model, with additive error, additive trend and no seasonality. It produced an AIC of 3222, BIC of 3245 and MASE of 0.2468. The residual check was similar to the previous model, showing errors to be random with roughly uniform variance until around 2005, after which it increased. The ACF of residuals had a few more significant lags at later lag periods, but the histogram also appeared to indicate normally distributed errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit3.etsA = ets(cop, model = "AAA")
summary(fit3.etsA)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.1412 
##     gamma = 1e-04 
##     phi   = 0.8001 
## 
##   Initial states:
##     l = 14.1892 
##     b = -0.064 
##     s = -2.5215 -1.2214 0.583 1.7441 1.5899 2.2624
##            1.0995 0.5277 0.8193 -0.1246 -2.0798 -2.6786
## 
##   sigma:  4.2768
## 
##      AIC     AICc      BIC 
## 3262.334 3264.294 3332.679 
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE    MAPE      MASE
## Training set 0.0509217 4.176879 2.824332 -0.06622927 7.32543 0.2536477
##                   ACF1
## Training set 0.2511475

checkresiduals(fit3.etsA)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 81.758, df = 7, p-value = 5.995e-15
## 
## Model df: 17.   Total lags used: 24

The third state space model to be explored was the “AAA” model, with additive error, trend and seasonality. It produced an AIC of 3262, BIC of 3333 and MASE of 0.2511. The residual check showed similarly random errors, with stable variance until around 2005. The ACF of residuals had one quite significant lag and a few other borderline lags, but the histogram indicates a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit3.etsM = ets(cop, model = "MAA")
summary(fit3.etsM)

## ETS(M,A,A) 
## 
## Call:
##  ets(y = cop, model = "MAA") 
## 
##   Smoothing parameters:
##     alpha = 0.9998 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 15.211 
##     b = 0.266 
##     s = -2.1088 -1.1073 0.2377 0.865 1.511 1.6105
##            1.8953 1.5146 0.7141 -0.7249 -2.7785 -1.6286
## 
##   sigma:  0.0973
## 
##      AIC     AICc      BIC 
## 3076.478 3078.227 3142.915 
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.1868352 4.276536 2.772854 -1.077244 7.020437 0.2490246
##                   ACF1
## Training set 0.3643574

checkresiduals(fit3.etsM)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,A,A)
## Q* = 72.22, df = 8, p-value = 1.774e-12
## 
## Model df: 16.   Total lags used: 24

The fourth state space model to be explored was the “MAA” model, with multiplicative error, additive trend and additive seasonality. It produced an AIC of 3075, BIC of 3142 and MASE of 0.2520. The residual check showed errors with possible seasonality but rather uniform variance across the entire period. The ACF of residuals again had one quite significant lag and two other borderline significant lags, but the histogram indicated a perfect normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit4.etsM = ets(cop, model="MAM")
summary(fit4.etsM)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = cop, model = "MAM") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9734 
## 
##   Initial states:
##     l = 14.2355 
##     b = 0.1358 
##     s = 0.9885 1.0081 1.0304 1.0315 1.0193 0.9973
##            0.9894 0.9974 0.9926 0.9856 0.9672 0.9928
## 
##   sigma:  0.0934
## 
##      AIC     AICc      BIC 
## 3041.105 3043.064 3111.450 
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE    MAPE      MASE
## Training set 0.05546565 4.358848 2.776126 -0.1745263 6.63985 0.2493184
##                   ACF1
## Training set 0.3797468

checkresiduals(fit4.etsM)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,M)
## Q* = 50.228, df = 7, p-value = 1.303e-08
## 
## Model df: 17.   Total lags used: 24

The fifth state space model to be explored was the “MAM” model, with multiplicative error, additive trend and multiplicative seasonality. It produced an AIC of 3037, BIC of 3107 and MASE of 0.2510. The residual check again showed errors with possible seasonality but rather uniform variance across the entire period with a couple of spikes. The ACF of residuals once again had one quite significant lag and two other borderline significant lags like the previous two, but the histogram indicated a decent normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit5 = ets(cop)
summary(fit5)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = cop) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 23.1734 
## 
##   sigma:  0.088
## 
##      AIC     AICc      BIC 
## 2982.795 2982.861 2994.519 
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.05855649 4.315354 2.726028 -0.214157 6.582553 0.2448193
##                   ACF1
## Training set 0.3959303

The sixth state space model left out the model specification, to let the function itself suggest an optimal model. The resulting model was of type “MNN”, with multiplicative error, no trend and no seasonality. It produced an AIC of 2982, BIC of 2994 and MASE of 0.2448. The residual check again showed errors with possible seasonality but rather uniform variance across the entire period with a couple of spikes. The ACF of residuals showed signs of seasonality and had several significant lags. The histogram indicated a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.AAA1 = ets(cop, model = "AAA", bounds="admiss" )
summary(fit.AAA1)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA", bounds = "admiss") 
## 
##   Smoothing parameters:
##     alpha = 1.3141 
##     beta  = -0.0133 
##     gamma = 0 
##     phi   = 0.9899 
## 
##   Initial states:
##     l = 14.7396 
##     b = 0.1363 
##     s = -2.5427 -0.8006 0.3966 1.4113 1.5303 1.6275
##            1.3721 1.0329 0.7333 -0.4013 -2.2255 -2.1341
## 
##   sigma:  4.099
## 
##      AIC     AICc      BIC 
## 3231.077 3233.037 3301.422 
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.1707125 4.003205 2.654505 -0.07171114 6.632331 0.2383959
##                    ACF1
## Training set 0.04577771

checkresiduals(fit.AAA1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 53.441, df = 7, p-value = 3.034e-09
## 
## Model df: 17.   Total lags used: 24

The seventh state space model once more was of type “AAA”, with additive error, additive trend and additive seasonality. The stability condition constraints were relaxed in this model, producing an AIC of 3231, BIC of 3301 and MASE of 0.2380. The residual check was similar to the original AAA model, with random errors of low uniform variance until 2005 and higher variance afterward. The ACF of residuals showed several significant lags. The histogram indicated a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.AAA3 = ets(cop, model = "AAA", upper=rep(1,4))
summary(fit.AAA3)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA", upper = rep(1, 4)) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.0406 
##     gamma = 1e-04 
##     phi   = 0.8097 
## 
##   Initial states:
##     l = 14.3128 
##     b = 0.3818 
##     s = -2.0552 -0.7183 0.4851 1.1429 1.5037 1.4076
##            1.7882 1.5333 0.708 -0.9923 -2.4227 -2.3804
## 
##   sigma:  4.3349
## 
##      AIC     AICc      BIC 
## 3272.257 3274.216 3342.602 
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.06331985 4.233573 2.754155 -0.1177275 6.860767 0.2473453
##                   ACF1
## Training set 0.3334076

checkresiduals(fit.AAA3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 100.86, df = 7, p-value < 2.2e-16
## 
## Model df: 17.   Total lags used: 24

The eighth state space model was another variation of type “AAA”, with additive error, additive trend and additive seasonality. This time the upper argument was used to ensure stability by constraining alpha and beta to be non-negative and that 2*alpha + beta was less than or equal to 4. It produced an AIC of 3271, BIC of 3341 and MASE of 0.2497. The residual check was pretty much identical to the original AAA model, with random errors of low uniform variance until 2005 and higher variance afterward. The ACF of residuals showed one significant lag and a couple borderline. The histogram indicated a roughly normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.AAA4 = ets(cop, model = "AAA", opt.crit = "mse")
summary(fit.AAA4)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA", opt.crit = "mse") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.1412 
##     gamma = 1e-04 
##     phi   = 0.8001 
## 
##   Initial states:
##     l = 14.1892 
##     b = -0.064 
##     s = -2.5215 -1.2214 0.583 1.7441 1.5899 2.2624
##            1.0995 0.5277 0.8193 -0.1246 -2.0798 -2.6786
## 
##   sigma:  4.2768
## 
##      AIC     AICc      BIC 
## 3262.334 3264.294 3332.679 
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE    MAPE      MASE
## Training set 0.0509217 4.176879 2.824332 -0.06622927 7.32543 0.2536477
##                   ACF1
## Training set 0.2511475

checkresiduals(fit.AAA4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 81.758, df = 7, p-value = 5.995e-15
## 
## Model df: 17.   Total lags used: 24

The ninth state space model was yet another variation of type “AAA”, with additive error, additive trend and additive seasonality. This time the optimality criterion was set to optimise for MSE. It produced an AIC of 3262, BIC of 3333 and MASE of 0.2511. The residual check was also more or less the same as the original AAA model, with random errors of low uniform variance until 2005 and higher variance afterward. The ACF of residuals showed one significant lag and a couple borderline. The histogram indicated a roughly normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.AAdA1 = ets(cop, model = "AAA", damped = TRUE, bounds="admiss" )
summary(fit.AAdA1)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA", damped = TRUE, bounds = "admiss") 
## 
##   Smoothing parameters:
##     alpha = 1.3141 
##     beta  = -0.0133 
##     gamma = 0 
##     phi   = 0.9899 
## 
##   Initial states:
##     l = 14.7396 
##     b = 0.1363 
##     s = -2.5427 -0.8006 0.3966 1.4113 1.5303 1.6275
##            1.3721 1.0329 0.7333 -0.4013 -2.2255 -2.1341
## 
##   sigma:  4.099
## 
##      AIC     AICc      BIC 
## 3231.077 3233.037 3301.422 
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.1707125 4.003205 2.654505 -0.07171114 6.632331 0.2383959
##                    ACF1
## Training set 0.04577771

checkresiduals(fit.AAdA1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 53.441, df = 7, p-value = 3.034e-09
## 
## Model df: 17.   Total lags used: 24

The tenth state space model once more was of type “AAA”, with additive error, additive trend and additive seasonality. The stability condition constraints were again relaxed in this model but also a damping factor was included, producing an AIC of 3231, BIC of 3301 and MASE of 0.2380. The residual check was similar to the original AAA model, with random errors of low uniform variance until 2005 and higher variance afterward. The ACF of residuals showed several significant lags across the entire lag spectrum. The histogram indicated a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.AAdA4 = ets(cop, model = "AAA", damped = TRUE, opt.crit = "mse")
summary(fit.AAdA4)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = cop, model = "AAA", damped = TRUE, opt.crit = "mse") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.1412 
##     gamma = 1e-04 
##     phi   = 0.8001 
## 
##   Initial states:
##     l = 14.1892 
##     b = -0.064 
##     s = -2.5215 -1.2214 0.583 1.7441 1.5899 2.2624
##            1.0995 0.5277 0.8193 -0.1246 -2.0798 -2.6786
## 
##   sigma:  4.2768
## 
##      AIC     AICc      BIC 
## 3262.334 3264.294 3332.679 
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE    MAPE      MASE
## Training set 0.0509217 4.176879 2.824332 -0.06622927 7.32543 0.2536477
##                   ACF1
## Training set 0.2511475

checkresiduals(fit.AAdA4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 81.758, df = 7, p-value = 5.995e-15
## 
## Model df: 17.   Total lags used: 24

The eleventh state space model once more was of type “AAA”, with additive error, additive trend and additive seasonality. The optimality criterion was set to optimise for MSE and a damping factor was included, producing an AIC of 3262, BIC of 3333 and MASE of 0.2511. The residual check was similar to the original AAA model, with random errors of low uniform variance until 2005 and higher variance afterward. The ACF of residuals showed two highly significant lags. The histogram indicated a normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

fit.MAdM = ets(cop, model = "MAM", damped = TRUE, opt.crit = "mse")
summary(fit.MAdM)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = cop, model = "MAM", damped = TRUE, opt.crit = "mse") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.0502 
##     gamma = 1e-04 
##     phi   = 0.8001 
## 
##   Initial states:
##     l = 14.3303 
##     b = 0.0533 
##     s = 0.9289 0.9482 0.966 1.0064 1.0321 1.0649
##            1.0551 1.037 1.0426 1.0159 0.9594 0.9435
## 
##   sigma:  0.0988
## 
##      AIC     AICc      BIC 
## 3081.691 3083.651 3152.037 
## 
## Training set error measures:
##                      ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.07716276 4.112226 2.742166 -0.09763832 6.872172 0.2462686
##                   ACF1
## Training set 0.3159422

checkresiduals(fit.MAdM)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,M)
## Q* = 51.663, df = 7, p-value = 6.804e-09
## 
## Model df: 17.   Total lags used: 24

The eleventh state space model once more was of type “MAM”, with multiplicative error, additive trend and multiplicative seasonality. The optimality criterion was set to optimise for MSE and a damping factor was included, producing an AIC of 3075, BIC of 3145 and MASE of 0.2448. The residual check showed random errors of roughly similar variance across the whole period. The ACF of residuals was rather different than the other models, with only one significant lag albeit some signs of seasonality. The histogram indicated a somewhat bimodal but normally distributed set of errors. The Ljung-Box test result indicated the model was significant with a p-value of <0.01

aic.models = AIC(fit1.etsA,fit2.etsA, fit3.etsA, fit3.etsM, fit4.etsM, fit.AAA1, fit.AAA3, fit.AAA4, fit.MAdM,fit.AAdA1, fit.AAdA4)
bic.models = BIC(fit1.etsA,fit2.etsA, fit3.etsA, fit3.etsM, fit4.etsM, fit.AAA1, fit.AAA3, fit.AAA4, fit.MAdM, fit.AAdA1, fit.AAdA4)

The above and below code snippets collect all models’ AIC and BIC values, to return a sorted list of the best models by these metrics.

sort.score(aic.models, score = "aic")

##           df      AIC
## fit4.etsM 18 3041.105
## fit3.etsM 17 3076.478
## fit.MAdM  18 3081.691
## fit2.etsA  6 3223.324
## fit.AAA1  18 3231.077
## fit.AAdA1 18 3231.077
## fit1.etsA  3 3256.335
## fit3.etsA 18 3262.334
## fit.AAA4  18 3262.334
## fit.AAdA4 18 3262.334
## fit.AAA3  18 3272.257

sort.score(bic.models, score = "bic")

##           df      BIC
## fit4.etsM 18 3111.450
## fit3.etsM 17 3142.915
## fit.MAdM  18 3152.037
## fit2.etsA  6 3246.772
## fit1.etsA  3 3268.059
## fit.AAA1  18 3301.422
## fit.AAdA1 18 3301.422
## fit3.etsA 18 3332.679
## fit.AAA4  18 3332.679
## fit.AAdA4 18 3332.679
## fit.AAA3  18 3342.602

forecast.fit4.etsM = forecast(fit4.etsM, h = 10)
forecast.fit.AAdA1 = forecast(fit.AAdA1, h = 10)
forecast.fit4.etsM

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 2016       45.25335 39.83476 50.67194 36.96633 53.54037
## Oct 2016       45.20101 37.52990 52.87213 33.46906 56.93297
## Nov 2016       44.21947 35.00768 53.43125 30.13125 58.30768
## Dec 2016       43.35610 32.90339 53.80881 27.37005 59.34214
## Jan 2017       43.55179 31.78601 55.31757 25.55758 61.54600
## Feb 2017       42.42077 29.83828 55.00326 23.17752 61.66402
## Mar 2017       43.22523 29.34546 57.10499 21.99796 64.45250
## Apr 2017       43.53229 28.55488 58.50970 20.62631 66.43826
## May 2017       43.74277 27.74376 59.74179 19.27439 68.21116
## Jun 2017       43.38789 26.62224 60.15355 17.74704 69.02875

forecast.fit.AAdA1

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 2016       44.96109 39.70801 50.21418 36.92720 52.99499
## Oct 2016       43.79857 35.17900 52.41814 30.61608 56.98106
## Nov 2016       42.45495 31.49733 53.41256 25.69672 59.21317
## Dec 2016       40.56795 27.72542 53.41048 20.92700 60.20890
## Jan 2017       40.83319 26.37986 55.28653 18.72872 62.93766
## Feb 2017       40.59981 24.72554 56.47408 16.32221 64.87741
## Mar 2017       42.28350 25.13042 59.43659 16.05012 68.51688
## Apr 2017       43.27894 24.95890 61.59898 15.26085 71.29702
## May 2017       43.44089 24.04512 62.83666 13.77762 73.10416
## Jun 2017       43.64374 23.24862 64.03886 12.45210 74.83538

The mean forecasts shown in the above tables and plotted below clearly indicate very closely aligned predictions from these three models. The forecasts for the tenth period ahead, clearly susceptible to the wide variance associated with such long-range forecasts, were nevertheless quite similar. The only comment one may make on these are that the long-range forecasts appear to be decreasing for the MAdM model, staying roughly the same for the AAdA model and somewhat increasing for the Holt Winters Additive Seasonality model. The confidence intervals for all three models also track relatively closely with each other throughout the forecast period.

plot(cop, xlim = c(1986, 2018),  ylab = "Crude Oil Price", xlab = "Year", 
     main="Fig. 22 - 10 Month Ahead Crude Oil Price Forecast")                       
lines(fit3.hw$mean, col="purple", type = "o")
lines(fitted(fit3.hw), col="purple")
lines(fitted(fit.AAdA1), col="red")
lines(forecast.fit.AAdA1$mean, col = "red", type = "o")
lines(fitted(fit4.etsM), col="green3", lty=1)
lines(forecast.fit4.etsM$mean, col = "green3", type = "o")
legend("topleft",lty=1, pch = 1, text.width = 12, col=c("black","purple", "green3", "red"), 
       c("Crude Oil Price","Holt-Winter's Additive", "ETS MAdM", "ETS AAdA"))

Conclusion

The models that appear to produce the best results, according to their AIC, BIC and MASE metrics, are the ETS(MAdM), ETS(AAdA) and Holt Winters (Additive Seasonality) models. The above plot of each of these, with forecasts, overlaid on the original series shows that all three track the historical crude price rather closely. The forecasts are also tracking relatively closely to each other.

The best three models had the following metrics: 1. ETS(MAdM): AIC 3037, BIC 3108, MASE 0.251 2. ETS(AAdA): AIC 3231, BIC 3301, MASE 0.238 3. Holt Winters Additive Seasonality: AIC 3279, BIC 3346, MASE 0.252

It is interesting that the exponential smoothing models did not perform as well as the state space models, given the high correlation of the price series to the predictor series of OPEC crude oil production.