Create user functions
Introduction
Method
Data description
Data exploration and visualisation
Decomposition
Model fitting
Data transformation and differencing
- Model fitting on transformed and differenced data
Model selection
Forecasting
Conclusion

library(TSA)
library(forecast)
library(x12)
library(tseries)
library(dLagM)
library(readr)
library(stats)
library(car)
library(dynlm)
library(lmtest)
library(ggplot2)

Create user functions

sum.res = function(x){
  summary(x)
  checkresiduals(x$model)
  shapiro.test(x$model$residuals)
}

sum.res2 = function(x){
  summary(x)
  checkresiduals(x)
  shapiro.test(x$residuals)
}

explore=function(x){
  plot(x, xlab='Year', main="Time series plot of energy production series", ylab = "Industrial production (IP) index")
  points(y=x,x=time(x), pch=as.vector(season(x)))
  acf(x, lag.max = 48, main="Sample ACF plot")
  adf.test(x, k = ar(x)$order)
}

customise <- function(dataset, type){
  if (type == "hw"){
    seasonal = c("additive","multiplicative")
    dampeds = c(TRUE,FALSE)
    expand = expand.grid(seasonal,dampeds)
    fit.AIC = array(NA,4)
    fit.BIC = array(NA,4)
    fit.MASE=array(NA,4)
    levels = array(NA, dim=c(4,2))
    
    for(i in 1:4){fit.AIC[i]= hw(dataset, 
                                 seasonal = toString(expand[i ,1]),                    
                                 damped = expand[i ,2])$model$aic
    fit.BIC[i]= hw(dataset, 
                   seasonal = toString(expand[i ,1]),                    
                   damped = expand[i ,2])$model$bic
    fit.MASE[i]= accuracy(hw(dataset, 
                             seasonal = toString(expand[i ,1]),                    
                             damped = expand[i ,2]))[,6]
    levels[i,1]= toString(expand[i ,1])  
    levels[i,2]= expand[i ,2]
    }
    results = data.frame(levels, fit.MASE, fit.AIC, fit.BIC)
    colnames(results)= c("seasonal","damped","MASE", "AIC", "BIC")
    results
    } else if (type == "ets") {
      models = c("AAA","MAM","MMM", "MAA")
      dampeds = c(TRUE,FALSE)
      bounds = c("both", "admissible")
      expand = expand.grid(models,dampeds, bounds)
      fit.AIC = array(NA,16)
      fit.BIC = array(NA,16)
      fit.MASE=array(NA,16)
      levels = array(NA, dim=c(16,3))
      
      for(i in 1:16){fit.AIC[i]= ets(dataset, 
                                     model = toString(expand[i ,1]),                    
                                     damped = expand[i ,2],
                                     bounds = toString(expand[i, 3]))$aic 
      fit.BIC[i]= ets(dataset, 
                      model = toString(expand[i ,1]),                    
                      damped = expand[i ,2],
                      bounds = toString(expand[i, 3]))$bic
      fit.MASE[i]= accuracy(ets(dataset, 
                                model = toString(expand[i ,1]),                    
                                damped = expand[i ,2],
                                bounds = toString(expand[i, 3])))[,6]
      levels[i,1]= toString(expand[i ,1])  
      levels[i,2]= expand[i ,2]
      levels[i,3]= toString(expand[i ,3])
      }
      results = data.frame(levels, fit.MASE, fit.AIC, fit.BIC)
      colnames(results)= c("Model","Damped", "Bounds", "MASE","AIC", "BIC") 
      results
      } else {
      print("Specify correct model type `hw` or `ets`")
      }
  }

Introduction

The dataset used in this report energy is obtained from https://www.kaggle.com/sadeght/industrial-production-electric-and-gas-utilities. It contains information on monthly Industrial production of electric and gas utilities in the United States between January 1939 and May 2019. The aim of this investigation is to analyse and forecast the Industrial production (IP) index for the next 10 months.

Method

Firstly, a time series will be created for variable IPG2211A2N in dataset energy. Then we will examine them by looking at their stationarity and decompositions. After exploring the properties of this time series we will proceed to modelling. Three types of models will be used including dynamic models, exponential smoothing models and state-space models. These models will be fitted to both the raw data and data after transformation and differencing. The best model will be decided based on significance of the model and terms, MASE, information criteria and residual analysis. Finally, the best model chosen will be used to make 10 months ahead forecasts.

Data description

The following code chunk is used to load dataset energy and convert it to time series:

energy <- read_csv("energy.csv")
original<- ts(energy$IPG2211A2N, start=c(1939,1), frequency = 12)
plot(original, xlab='Year', main="Time series plot of energy production series", ylab = "Industrial production (IP) index")

Figure 1: Time series plot and sample ACF plot of complete monthly energy production

We decided to focus on series starting from 1970 Jan because prior to 1970 there was very little energy production activities nor change in patterns and the chance of it affecting intended forecast period is slim. Data from 1970 onward is also more relevant for prediction purpose.

energy <- energy$IPG2211A2N[373:965]
energy <- ts(energy, start=c(1970,1), frequency = 12)

Data exploration and visualisation

The following code chunk is used to explore and visualise the time series plot and sample ACF plot of monthly energy production:

explore(energy)

Figure 2: Time series plot and sample ACF plot of monthly energy production

## 
##  Augmented Dickey-Fuller Test
## 
## data:  x
## Dickey-Fuller = -0.5817, Lag order = 25, p-value = 0.9781
## alternative hypothesis: stationary

The time series plots for monthly energy production from Figure 2 shows that:

The series has an obvious upward trend and seasonality.
This series has no intervention point over the given time period.
It is hard to tell the behaviour or the existence of changing variance without removing the seasonal component first.

The sample ACF plot and the Augmented Dickey-Fuller test results for monthly energy production from Figure 2 shows that:

A decaying pattern in the ACF plot of monthly energy production suggests the existence of trend.
Obvious seasonality was observed in the ACF plot, this supports previous observation based on time series plot.
Augmented Dickey-Fuller test gives a p-value of 0.9781 (> 0.05), so there is not enough evidence to reject the null hypothesis of non-stationarity.

Based on the above evidence and observations made on the time series plot, energy is a non-stationary time series with seasonality.

Decomposition

In order to better investigate seasonal, trend and other components in the time series, X12 and STL decompositions are used. The following code chunk is used to perform and visualise X12 decomposition of energy series:

energy.x12 = x12(energy)
plot(energy.x12, sa = T, trend = T, main = "X12 decomposition of monthly energy production")

Figure 3: X12 decomposition of monthly energy production

Figure 3 shows the X12 decomposition of energy series, we can observe that:

The seasonally adjusted line is very different from the original, which means that seasonality is a main cause for variation in the original series.
The trend line shows a general upward trend.

The following code chunk is used to perform and visualise STL decomposition of energy series:

energy.decom <- stl(energy, t.window=15, s.window="periodic", robust=TRUE) 
plot(energy.decom, main="STL decomposition of monthly energy production time series")

Figure 4: STL decomposition of monthly energy production

Figure 4 shows the STL decomposition of energy series, we can observe that:

There a general upward trend in the trend component.
There is changing variance in the remainder.

Model fitting

Dynamic linear regression method

We will use dynamic models without specifying the intervention function because there was no intervention point found in the series. Dynamic models fitted include lags of the series itself as well as trend and seasonality.

The following code implements a dynamic model with first lag of energy series and performs residual analysis.

Y.t = energy
model1 = dynlm(Y.t ~ L(Y.t , 1 ) + season(Y.t) + trend(Y.t))
sum.res2(model1)

Figure 5: Residual analysis plots for model1

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.98909, p-value = 0.0002179

summary(model1)

## 
## Time series regression with "ts" data:
## Start = 1970(2), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + season(Y.t) + trend(Y.t))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.4732  -1.7672   0.0871   2.0916  10.4218 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     18.79238    1.42373  13.199  < 2e-16 ***
## L(Y.t, 1)        0.75830    0.02710  27.986  < 2e-16 ***
## season(Y.t)Feb -11.51199    0.65866 -17.478  < 2e-16 ***
## season(Y.t)Mar -12.12588    0.63531 -19.087  < 2e-16 ***
## season(Y.t)Apr -15.86832    0.65211 -24.334  < 2e-16 ***
## season(Y.t)May -10.42282    0.73053 -14.268  < 2e-16 ***
## season(Y.t)Jun  -3.68305    0.73705  -4.997 7.73e-07 ***
## season(Y.t)Jul  -2.24870    0.66852  -3.364  0.00082 ***
## season(Y.t)Aug  -6.43639    0.63923 -10.069  < 2e-16 ***
## season(Y.t)Sep -13.87665    0.63872 -21.726  < 2e-16 ***
## season(Y.t)Oct -13.79282    0.67140 -20.543  < 2e-16 ***
## season(Y.t)Nov  -7.27006    0.72826  -9.983  < 2e-16 ***
## season(Y.t)Dec   1.21025    0.69790   1.734  0.08343 .  
## trend(Y.t)       0.33453    0.03917   8.541  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.16 on 578 degrees of freedom
## Multiple R-squared:  0.9788, Adjusted R-squared:  0.9783 
## F-statistic:  2053 on 13 and 578 DF,  p-value: < 2.2e-16

MASE(lm(model1))

##                 MASE
## lm(model1) 0.4230642

From the model summary, we can conclude that all terms of the dynamic model with the first lag are significant at 5% level except the seasonal term for December (p-value = 0.08343 > 0.05). The model is reported to be overall statistically significant at 5% level (p-value < 0.05) and its adjusted \(R^2=0.9783\) which means the model explains 97.83% of the variability in energy series.

According to the residual analysis in Figure 5, we can conclude the following:

The errors are not spread randomly, there is evidence of changing variance with the residuals being smaller in the center and larger towards the ends of the plot.
All the lags in ACF plot are significant and have a wave-like pattern, which suggests serial correlation and seasonality still remaining in the residuals.
The errors are not normal. The histogram and the Shapiro-Wilk normality test with p-value = 0.0002179 suggest not normal residuals.

Overall, we can conclude that this model is not successful at capturing the autocorrelation and seasonality in the series.

To overcome the previously found issues in the residual diagnostic, we will include the second lag of energy series as another explanatory variable.

The following code implements the dynamic model with both lags of the series plus trend and seasonality components.

model2 = dynlm(Y.t ~ L(Y.t ,1) + L(Y.t ,2) + season(Y.t) + trend(Y.t))
sum.res2(model2)

Figure 6: Residual analysis plots for model2

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.98134, p-value = 7.389e-07

summary(model2)

## 
## Time series regression with "ts" data:
## Start = 1970(3), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + season(Y.t) + trend(Y.t))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4417  -1.8703   0.1242   1.7371   9.9693 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     19.07366    1.38261  13.795  < 2e-16 ***
## L(Y.t, 1)        0.96504    0.04009  24.073  < 2e-16 ***
## L(Y.t, 2)       -0.27166    0.04027  -6.747 3.69e-11 ***
## season(Y.t)Feb -10.11926    0.67554 -14.979  < 2e-16 ***
## season(Y.t)Mar  -7.56372    0.91223  -8.291 7.92e-16 ***
## season(Y.t)Apr -11.91309    0.85826 -13.880  < 2e-16 ***
## season(Y.t)May  -6.32096    0.92788  -6.812 2.42e-11 ***
## season(Y.t)Jun  -1.57203    0.77523  -2.028  0.04304 *  
## season(Y.t)Jul  -1.60296    0.65135  -2.461  0.01415 *  
## season(Y.t)Aug  -5.36601    0.63608  -8.436 2.66e-16 ***
## season(Y.t)Sep -11.22891    0.72992 -15.384  < 2e-16 ***
## season(Y.t)Oct  -9.54850    0.90081 -10.600  < 2e-16 ***
## season(Y.t)Nov  -3.85370    0.86355  -4.463 9.74e-06 ***
## season(Y.t)Dec   2.67298    0.70619   3.785  0.00017 ***
## trend(Y.t)       0.42783    0.04050  10.564  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 576 degrees of freedom
## Multiple R-squared:  0.9803, Adjusted R-squared:  0.9798 
## F-statistic:  2046 on 14 and 576 DF,  p-value: < 2.2e-16

MASE(lm(model2))

##                 MASE
## lm(model2) 0.4013518

From the model summary, we can conclude that all terms of this dynamic model are statistically significant at 5% level. The model is reported to be overall statistically significant at 5% level (p-value < 0.05) and its adjusted \(R^2=0.9798\) which means the model explains 97.98% of the variability in energy series.

According to the residual analysis in Figure 6, we can conclude the following:

The errors are not spread randomly, in fact there seems to be even stronger changing variance than in the previous model; the residuals are much smaller in the center and larger towards the ends of the plot.
There is still seasonality left in the residuals from this model with significant seasonal lags observed in the ACF plot.
The errors are not normal. The histogram is left-skewed, and the Shapiro-Wilk normality test with p-value = 7.389e-07 suggests not normal residuals.

Overall, we can conclude that this model also fails to capture the autocorrelation and seasonality in the series.

We will include the third lag of the series to see if it makes any improvement in the residuals.

model3 = dynlm(Y.t ~ L(Y.t ,1) + L(Y.t ,2) + L(Y.t ,3) + season(Y.t) + trend(Y.t))
sum.res2(model3)

Figure 7: Residual analysis plots for model3

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.98113, p-value = 6.601e-07

summary(model3)

## 
## Time series regression with "ts" data:
## Start = 1970(4), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + L(Y.t, 3) + season(Y.t) + 
##     trend(Y.t))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9968  -1.6394   0.0416   1.8212  10.1707 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     16.78659    1.44920  11.583  < 2e-16 ***
## L(Y.t, 1)        1.01620    0.04100  24.786  < 2e-16 ***
## L(Y.t, 2)       -0.45479    0.05632  -8.075 3.98e-15 ***
## L(Y.t, 3)        0.19064    0.04129   4.617 4.81e-06 ***
## season(Y.t)Feb  -9.03866    0.70570 -12.808  < 2e-16 ***
## season(Y.t)Mar  -6.92341    0.91764  -7.545 1.78e-13 ***
## season(Y.t)Apr -13.41336    0.90504 -14.821  < 2e-16 ***
## season(Y.t)May  -7.19991    0.93270  -7.719 5.22e-14 ***
## season(Y.t)Jun  -2.80260    0.80768  -3.470  0.00056 ***
## season(Y.t)Jul  -1.80484    0.64218  -2.811  0.00512 ** 
## season(Y.t)Aug  -4.61774    0.64660  -7.142 2.81e-12 ***
## season(Y.t)Sep -10.56904    0.73332 -14.413  < 2e-16 ***
## season(Y.t)Oct  -9.61532    0.88745 -10.835  < 2e-16 ***
## season(Y.t)Nov  -5.04003    0.88747  -5.679 2.15e-08 ***
## season(Y.t)Dec   1.73780    0.72350   2.402  0.01663 *  
## trend(Y.t)       0.34493    0.04374   7.886 1.58e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.996 on 574 degrees of freedom
## Multiple R-squared:  0.9809, Adjusted R-squared:  0.9804 
## F-statistic:  1965 on 15 and 574 DF,  p-value: < 2.2e-16

MASE(lm(model3))

##                 MASE
## lm(model3) 0.3928438

From the model summary, we can conclude that all terms of this dynamic model are statistically significant at 5% level. The model is reported to be overall statistically significant at 5% level (p-value < 0.05) and its adjusted \(R^2=0.9804\) which means the model explains 98.04% of the variability in energy series.

According to the residual analysis in Figure 7, we can conclude the following:

The errors are not spread randomly, in fact there seems to be even stronger changing variance than in the previous model; the residuals are much smaller in the center and larger towards the ends of the plot.
There is still seasonality left in the residuals from this model with significant seasonal lags observed in the ACF plot.
The errors are not normal. The histogram is left-skewed, and the Shapiro-Wilk normality test with p-value = 2.2e-16 suggests not normal residuals.

Overall, we can conclude that there is no improvement in the residuals from this model, the model fails to capture the autocorrelation and seasonality in the series.

Exponential smoothing method

Another forecasting method we will try is exponential smoothing. Because we have found a strong seasonal component in energy series that we want to make predictions for, we will only consider Holt-Winters’ models that include either additive or multiplicative seasonality.

To find the best model from the set of candidate Holt-Winters’ models, we use a function customise.

customise(energy, type ="hw")

##         seasonal damped      MASE      AIC      BIC
## 1       additive   TRUE 0.7143593 4913.613 4992.547
## 2 multiplicative   TRUE 0.6438479 4702.382 4781.315
## 3       additive  FALSE 0.7154824 4910.398 4984.946
## 4 multiplicative  FALSE 0.6524627 4728.817 4803.366

According to the Information Criteria and MASE measure, we select to fit Holt-Winters’ model with damped trend and multiplicative seasonality. The following code implements this model and displays residual diagnostic plots.

hw.md <- hw(energy, seasonal = "multiplicative", damped = T, h=2*frequency(energy))
sum.res(hw.md)

## 
## Forecast method: Damped Holt-Winters' multiplicative method
## 
## Model Information:
## Damped Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = energy, h = 2 * frequency(energy), seasonal = "multiplicative",  
## 
##  Call:
##      damped = T) 
## 
##   Smoothing parameters:
##     alpha = 0.4924 
##     beta  = 1e-04 
##     gamma = 0.2427 
##     phi   = 0.9799 
## 
##   Initial states:
##     l = 39.8515 
##     b = 0.405 
##     s = 1.0804 1.0078 0.9655 0.9878 1.0114 0.9817
##            0.9396 0.8977 0.9224 1.0001 1.0745 1.1311
## 
##   sigma:  0.0283
## 
##      AIC     AICc      BIC 
## 4702.382 4703.574 4781.315 
## 
## Error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1526271 2.317136 1.740105 0.1209437 2.186815 0.6438479
##                   ACF1
## Training set 0.1721279
## 
## Forecasts:
##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## Jun 2019      104.50165 100.71190 108.29140  98.70573 110.29757
## Jul 2019      115.01424 110.36448 119.66399 107.90305 122.12543
## Aug 2019      114.03342 108.99290 119.07394 106.32461 121.74223
## Sep 2019      102.53372  97.64485 107.42260  95.05683 110.01061
## Oct 2019       94.49128  89.67941  99.30314  87.13216 101.85039
## Nov 2019       98.79101  93.45883 104.12318  90.63615 106.94586
## Dec 2019      112.43789 106.04436 118.83143 102.65982 122.21596
## Jan 2020      121.72293 114.46666 128.97920 110.62542 132.82044
## Feb 2020      108.42835 101.67942 115.17728  98.10675 118.74995
## Mar 2020      102.89353  96.22932 109.55774  92.70150 113.08556
## Apr 2020       89.75631  83.72475  95.78787  80.53184  98.98079
## May 2020       92.62609  86.18424  99.06794  82.77413 102.47805
## Jun 2020      104.51566  96.72836 112.30295  92.60602 116.42530
## Jul 2020      115.02957 106.21428 123.84486 101.54774 128.51139
## Aug 2020      114.04853 105.07233 123.02474 100.32061 127.77646
## Sep 2020      102.54724  94.26929 110.82518  89.88721 115.20726
## Oct 2020       94.50366  86.68887 102.31846  82.55196 106.45536
## Nov 2020       98.80389  90.44331 107.16446  86.01748 111.59029
## Dec 2020      112.45247 102.72524 122.17971  97.57595 127.32899
## Jan 2021      121.73863 110.98369 132.49358 105.29037 138.18690
## Feb 2021      108.44226  98.66605 118.21848  93.49083 123.39369
## Mar 2021      102.90667  93.44713 112.36621  88.43955 117.37378
## Apr 2021       89.76771  81.35982  98.17561  76.90894 102.62649
## May 2021       92.63780  83.80282 101.47279  79.12586 106.14974

Figure 8: Residual analysis plots for Holt-Winters’ multiplicative damped method

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt-Winters' multiplicative method
## Q* = 148.36, df = 7, p-value < 2.2e-16
## 
## Model df: 17.   Total lags used: 24

## 
##  Shapiro-Wilk normality test
## 
## data:  x$model$residuals
## W = 0.99717, p-value = 0.3982

Figure 8 shows residual analysis plots for the fitted model. We can observe that there is a slight improvement in the residuals compared to dynamic models. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly. They do not have as strong changing variance as before.
There is still seasonality left in the residuals from this model with significant lag at lag 12 observed in the ACF plot. In addition, the errors are autocorrelated, there are a number of significant lags in the ACF plot.
The errors are normal. The histogram and the Shapiro-Wilk normality test with p-value = 0.3982 suggest normal residuals.

Overall, we can conclude that Holt-Winters’ method has improved the residuals in terms of randomness and changing variance. However, this model is not successful at capturing the autocorrelation and seasonality in the series.

State-space models

For each exponential smoothing method, there are two corresponding state-space models (with additive or multiplicative errors). There are 8 state-space variations which include seasonality that we can implement in R (some combinations are forbidden due to their stability issues). We use the created function to fit a set of possible candidate ETS models and display their accuracy measures to select the best out of the set.

customise(energy, type ="ets")

##    Model Damped     Bounds      MASE      AIC      BIC
## 1    AAA   TRUE       both 0.7143593 4913.613 4992.547
## 2    MAM   TRUE       both 0.6434487 4696.925 4775.859
## 3    MMM   TRUE       both 0.6494204 4704.646 4783.580
## 4    MAA   TRUE       both 0.7041751 4881.960 4960.893
## 5    AAA  FALSE       both 0.7154824 4910.398 4984.946
## 6    MAM  FALSE       both 0.6633972 4733.071 4807.619
## 7    MMM  FALSE       both 0.6690150 4736.624 4811.172
## 8    MAA  FALSE       both 0.6777347 4778.423 4852.972
## 9    AAA   TRUE admissible 0.6652187 4847.302 4926.235
## 10   MAM   TRUE admissible 0.6426198 4694.253 4773.187
## 11   MMM   TRUE admissible 0.6417418 4698.158 4777.091
## 12   MAA   TRUE admissible 0.6829879 4809.238 4888.171
## 13   AAA  FALSE admissible 0.7216804 4919.003 4993.551
## 14   MAM  FALSE admissible 0.6633972 4733.071 4807.619
## 15   MMM  FALSE admissible 0.6498015 4701.887 4776.436
## 16   MAA  FALSE admissible 0.7237773 4931.756 5006.305

The model that minimises AIC and BIC is ETS(M,Ad,M) with the parameters lying in the admissible parameter region.

The following code implements ETS(M,Ad,M) and displays residual analysis results.

ets.MAdM = ets(energy, model="MAM", damped = T, bounds = "admissible")
sum.res2(ets.MAdM)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = energy, model = "MAM", damped = T, bounds = "admissible") 
## 
##   Smoothing parameters:
##     alpha = 0.5098 
##     beta  = -0.0015 
##     gamma = 0.2567 
##     phi   = 0.9973 
## 
##   Initial states:
##     l = 40.7824 
##     b = 0.212 
##     s = 1.0674 1.003 0.9606 1.0044 1.021 0.9898
##            0.9467 0.9024 0.932 1.0026 1.0638 1.1062
## 
##   sigma:  0.028
## 
##      AIC     AICc      BIC 
## 4694.253 4695.445 4773.187 
## 
## Training set error measures:
##                       ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.006655702 2.313386 1.736785 -0.07892702 2.181245 0.6426198
##                   ACF1
## Training set 0.1557571

Figure 9: Residual analysis plots for ETS(M,Ad,M)

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

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.9971, p-value = 0.3766

Figure 9 shows residual analysis plots for the fitted model. In fact, we can observe the same picture in the residuals as in the Holt-Winters’ method. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly. They do not show evidence of changing variance.
There is still seasonality left in the residuals from this model with significant lag at lag 12 observed in the ACF plot. In addition, the errors are autocorrelated, there are a number of significant lags in the ACF plot.
The errors are normal. The histogram and the Shapiro-Wilk normality test with p-value = 0.3766 suggest normal residuals.

Overall, we can conclude that this model is also not successful at capturing the autocorrelation and seasonality in the series.

The model that minimises MASE was ETS(M,Md,M) with the parameters in the admissible region. We will fit this model in attempt to overcome the issues in the residuals.

ets.MMdM = ets(energy, model="MMM", damped = T, bounds = "admissible")
sum.res2(ets.MMdM)

## ETS(M,Md,M) 
## 
## Call:
##  ets(y = energy, model = "MMM", damped = T, bounds = "admissible") 
## 
##   Smoothing parameters:
##     alpha = 0.5706 
##     beta  = -0.0081 
##     gamma = 0.252 
##     phi   = 0.9559 
## 
##   Initial states:
##     l = 40.5075 
##     b = 1.0112 
##     s = 1.0747 1.0025 0.9602 0.9971 1.0145 0.9841
##            0.9431 0.9066 0.9371 0.9966 1.0678 1.1158
## 
##   sigma:  0.0282
## 
##      AIC     AICc      BIC 
## 4698.158 4699.349 4777.091 
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.2307523 2.324038 1.734413 0.2365962 2.169279 0.6417418
##                   ACF1
## Training set 0.1153484

Figure 10: Residual analysis plots for ETS(M,Md,M)

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

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99782, p-value = 0.6454

According to the residual analysis plots in Figure 10, we can observe a very similar picture to the previous model. There are no improvements in terms of serial correlation and seasonality in the residuals.

The following ETS(M,M,M) model was fitted as well to see if there are any improvements in the residuals.

ets.MMM = ets(energy, model="MMM", damped = F)
sum.res2(ets.MMM)

## ETS(M,M,M) 
## 
## Call:
##  ets(y = energy, model = "MMM", damped = F) 
## 
##   Smoothing parameters:
##     alpha = 0.2838 
##     beta  = 1e-04 
##     gamma = 0.3338 
## 
##   Initial states:
##     l = 43.0285 
##     b = 1.0016 
##     s = 1.0354 0.9842 0.9751 0.9989 1.059 1.0399
##            0.9568 0.9135 0.9459 0.9982 1.0444 1.0489
## 
##   sigma:  0.0291
## 
##      AIC     AICc      BIC 
## 4736.624 4737.688 4811.172 
## 
## Training set error measures:
##                       ME     RMSE      MAE        MPE     MAPE     MASE
## Training set -0.07606488 2.397138 1.808123 -0.1168709 2.267875 0.669015
##                   ACF1
## Training set 0.3231045

Figure 11: Residual analysis plots for ETS(M,M,M)

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

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.98985, p-value = 0.0004088

Figure 11 shows residual analysis plots for the fitted model. In fact, we can observe that the seasonal lag 12 is no longer significant in the ACF plot which means this model is more successful than the previous ones in capturing seasonality in the data. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly. They do not show evidence of changing variance.
There is some autocorrelation left in the residuals according to significant lags in ACF.
The errors are not normal. Long tails in the histogram and the Shapiro-Wilk normality test with p-value = 0.0004088 suggest not normal residuals.

Overall, we can conclude that this model is performing better at capturing seasonality in the series.

The auto ETS model was fitted to see what the software automatically suggested model is.

fit.auto <- ets(energy)
summary(fit.auto)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = energy) 
## 
##   Smoothing parameters:
##     alpha = 0.5265 
##     beta  = 1e-04 
##     gamma = 0.2527 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 39.7118 
##     b = 0.415 
##     s = 1.0685 0.9986 0.963 0.9993 1.0189 0.9876
##            0.9458 0.9088 0.9315 0.999 1.0657 1.1135
## 
##   sigma:  0.0282
## 
##      AIC     AICc      BIC 
## 4696.925 4698.117 4775.859 
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1427332 2.321808 1.739026 0.1087827 2.177447 0.6434487
##                  ACF1
## Training set 0.147378

The model suggested automatically is ETS(M,Ad,M) which is the same model we have selected based on Information Criteria values.

Data transformation and differencing

Because we have found that there is some changing variance in the residuals as well as strong seasonal correlation, we will perform logarithmic transformation and differencing on the energy series in attempt to overcome these issues.

The following code chunk performs the log transformation on the data and displays it in the time series plot.

bc <- BoxCox.ar(energy, method = "yule-walker")

Figure 12: Time Series plot of log energy series

bc$ci

## [1] -0.1  0.3

energy.tr<- log(energy)
plot(energy.tr, xlab='Year', main="Time series plot of log energy production series", ylab = "log Industrial production (IP) index")

Figure 12: Time Series plot of log energy series

From the time series plot in Figure 12, we can see that the log transformation has smoothed our series a little bit.

We continue with the seasonal and ordinary differencing of the energy series.

energy.diff = diff(diff(energy.tr,lag = 12))
plot(energy.diff, xlab='Year', main = "Time series plot of transformed & differenced energy production series", ylab = "Industrial production (IP) index")

Figure 13: Time series plot of transformed & differenced series

Figure 13 shows the plot of transformed and differenced data. We can observe that the series is stationary. To confirm stationarity, we display the ACF plot and perform the Augmented Dickey-Fuller test:

acf(energy.diff, main = "Sample ACF plot of transformed & differenced energy production series")

Figure 14: ACF plot of transformed & differenced series

adf.test(energy.diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  energy.diff
## Dickey-Fuller = -11.642, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

According to Figure 14, there is no trend in the series but significant seasonal lags suggest the existence of seasonality. The results of ADF test report the data is stationary at 5% level of significance (p-value = 0.01). We will use this transformed and differenced series energy.diff for further model fitting.

Model fitting on transformed and differenced data

Dynamic linear regression method

The following code implements a dynamic model with first lag of energy.diff series and performs residual analysis.

Y.t = energy.diff
model1.1 = dynlm(Y.t ~ L(Y.t , 1 ) + season(Y.t) + trend(Y.t))
sum.res2(model1.1)

Figure 15: Residual analysis plots for model1.1

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99371, p-value = 0.01651

summary(model1.1)

## 
## Time series regression with "ts" data:
## Start = 1971(3), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + season(Y.t) + trend(Y.t))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.094797 -0.020389  0.000897  0.019744  0.122945 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.222e-04  5.368e-03   0.153    0.878    
## L(Y.t, 1)      -1.706e-01  4.152e-02  -4.108 4.58e-05 ***
## season(Y.t)Feb -2.271e-03  6.772e-03  -0.335    0.737    
## season(Y.t)Mar -1.765e-03  6.739e-03  -0.262    0.793    
## season(Y.t)Apr -3.944e-03  6.738e-03  -0.585    0.559    
## season(Y.t)May -2.995e-04  6.740e-03  -0.044    0.965    
## season(Y.t)Jun  7.532e-04  6.773e-03   0.111    0.911    
## season(Y.t)Jul  4.200e-04  6.773e-03   0.062    0.951    
## season(Y.t)Aug -1.358e-03  6.773e-03  -0.201    0.841    
## season(Y.t)Sep -2.586e-03  6.773e-03  -0.382    0.703    
## season(Y.t)Oct -8.536e-04  6.773e-03  -0.126    0.900    
## season(Y.t)Nov  8.085e-04  6.773e-03   0.119    0.905    
## season(Y.t)Dec  7.237e-04  6.772e-03   0.107    0.915    
## trend(Y.t)     -2.806e-06  9.900e-05  -0.028    0.977    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03318 on 565 degrees of freedom
## Multiple R-squared:  0.03077,    Adjusted R-squared:  0.008473 
## F-statistic:  1.38 on 13 and 565 DF,  p-value: 0.164

MASE(lm(model1.1))

##                   MASE
## lm(model1.1) 0.6487095

From the model summary, we can conclude that the first lag of the series is significant at 5% level but neither seasonal nor trend components are significant at 5% level which is expected since we use differenced series. The model is reported to be overall statistically insignificant at 5% level (p-value = 0.164) and its adjusted \(R^2=0.008473\) which means the model explains 0.85% of the variability in energy series. Obviously, this model has very low explainability and is not suitable.

Figure 15 shows residual analysis plots for the fitted model. In fact, we can observe slight improvement in the residuals in terms of serial correlation. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly. They do not show evidence of changing variance.
There is still seasonality left in the residuals from this model with significant lag at lag 12 and 24 observed in the ACF plot. There is slight autocorrelation in the residuals according to ACF plot.
The errors are normal. The histogram and the Shapiro-Wilk normality test with p-value = 0.3766 suggest normal residuals.

Overall, we can conclude that this model is not suitable for forecasting.

We will then include the second lag of the energy series and fit this model.

model2.1 = dynlm(Y.t ~ L(Y.t ,1) + L(Y.t ,2) + season(Y.t) + trend(Y.t))
sum.res2(model2.1)

Figure 16: Residual analysis plots for model2.1

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99673, p-value = 0.2927

summary(model2.1)

## 
## Time series regression with "ts" data:
## Start = 1971(4), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + season(Y.t) + trend(Y.t))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.095007 -0.019032  0.000184  0.020861  0.106291 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.379e-03  5.151e-03   0.268    0.789    
## L(Y.t, 1)      -2.187e-01  4.037e-02  -5.417 8.98e-08 ***
## L(Y.t, 2)      -2.904e-01  4.053e-02  -7.165 2.46e-12 ***
## season(Y.t)Feb -2.395e-03  6.494e-03  -0.369    0.712    
## season(Y.t)Mar -1.895e-03  6.495e-03  -0.292    0.771    
## season(Y.t)Apr -4.946e-03  6.463e-03  -0.765    0.444    
## season(Y.t)May -1.165e-03  6.464e-03  -0.180    0.857    
## season(Y.t)Jun -2.162e-04  6.496e-03  -0.033    0.973    
## season(Y.t)Jul  1.458e-05  6.494e-03   0.002    0.998    
## season(Y.t)Aug -1.384e-03  6.494e-03  -0.213    0.831    
## season(Y.t)Sep -2.859e-03  6.494e-03  -0.440    0.660    
## season(Y.t)Oct -1.660e-03  6.496e-03  -0.256    0.798    
## season(Y.t)Nov -1.795e-04  6.495e-03  -0.028    0.978    
## season(Y.t)Dec  3.498e-04  6.494e-03   0.054    0.957    
## trend(Y.t)     -5.015e-06  9.518e-05  -0.053    0.958    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03181 on 563 degrees of freedom
## Multiple R-squared:  0.1117, Adjusted R-squared:  0.08966 
## F-statistic: 5.059 on 14 and 563 DF,  p-value: 5.648e-09

MASE(lm(model2.1))

##                   MASE
## lm(model2.1) 0.6290276

From the model summary, we can conclude that the first and the second lags of the series is significant at 5% level but neither seasonal nor trend components are again significant at 5% level. The model is reported to be overall statistically significant at 5% level (p-value = 5.648e-09) and its adjusted \(R^2=0.08966\) which means the model explains about 9% of the variability in energy series. This is also very low explainability.

Figure 16 shows residual analysis plots for the fitted model. The picture that we can observe in the residuals from this model is very similar to the previous one except there is no serial correlation according to the ACF plot.

Overall, we can conclude that this model is also not suitable for forecasting due to low explainability and seasonality in the residuals.

Finally, we include the third lag of the energy series as another explanatory variable in attempt to overcome the issues in previous models.

model3.1 = dynlm(Y.t ~ L(Y.t ,1) + L(Y.t ,2) +L(Y.t ,3)+ season(Y.t) + trend(Y.t))
sum.res2(model3.1)

Figure 17: Residual analysis plots for model3.1

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99677, p-value = 0.303

summary(model3.1)

## 
## Time series regression with "ts" data:
## Start = 1971(5), End = 2019(5)
## 
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + L(Y.t, 3) + season(Y.t) + 
##     trend(Y.t))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.094989 -0.018952 -0.000352  0.021009  0.107281 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.390e-03  5.158e-03   0.269    0.788    
## L(Y.t, 1)      -2.308e-01  4.220e-02  -5.470 6.79e-08 ***
## L(Y.t, 2)      -2.996e-01  4.160e-02  -7.202 1.91e-12 ***
## L(Y.t, 3)      -4.259e-02  4.241e-02  -1.004    0.316    
## season(Y.t)Feb -2.349e-03  6.500e-03  -0.361    0.718    
## season(Y.t)Mar -1.895e-03  6.501e-03  -0.291    0.771    
## season(Y.t)Apr -5.048e-03  6.501e-03  -0.777    0.438    
## season(Y.t)May -1.313e-03  6.471e-03  -0.203    0.839    
## season(Y.t)Jun -2.727e-04  6.501e-03  -0.042    0.967    
## season(Y.t)Jul -7.192e-05  6.500e-03  -0.011    0.991    
## season(Y.t)Aug -1.392e-03  6.500e-03  -0.214    0.830    
## season(Y.t)Sep -2.832e-03  6.500e-03  -0.436    0.663    
## season(Y.t)Oct -1.685e-03  6.501e-03  -0.259    0.796    
## season(Y.t)Nov -2.623e-04  6.502e-03  -0.040    0.968    
## season(Y.t)Dec  2.607e-04  6.500e-03   0.040    0.968    
## trend(Y.t)     -3.940e-06  9.551e-05  -0.041    0.967    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03184 on 561 degrees of freedom
## Multiple R-squared:  0.1133, Adjusted R-squared:  0.08963 
## F-statistic: 4.781 on 15 and 561 DF,  p-value: 9.108e-09

MASE(lm(model3.1))

##                   MASE
## lm(model3.1) 0.6278907

From the model summary, we can conclude that the third lag of the series is not significant at 5% level so there is no point to further increase the number of lags. The model is reported to be overall statistically significant at 5% level (p-value = 9.108e-09) and its adjusted \(R^2=0.08963\) which means the model explains about 9% of the variability in energy series. This is also very low explainability.

Figure 17 shows residual analysis plots for the fitted model. The picture that we can observe in the residuals from this model is very similar to the previous one.

Overall, we can conclude that this model is also not suitable for forecasting due to low explainability and seasonality in the residuals.

Exponential smoothing method

We will try exponential smoothing on the transformed and differenced data to see whether there are any improvements in the residuals. To find the best model from the set of candidate Holt-Winters’ models, we use a function customise.

customise(energy.diff + 0.2, type = "hw")

##         seasonal damped      MASE       AIC       BIC
## 1       additive   TRUE 0.5894597 -218.5649 -140.0304
## 2 multiplicative   TRUE 0.5867337 -220.7570 -142.2225
## 3       additive  FALSE 0.5955193 -204.0216 -129.8501
## 4 multiplicative  FALSE 0.5897816 -217.3290 -143.1575

hw.md.diff <- hw(energy.diff+0.2, seasonal = "multiplicative", damped = T, h=2*frequency(energy.diff))
sum.res(hw.md.diff)

## 
## Forecast method: Damped Holt-Winters' multiplicative method
## 
## Model Information:
## Damped Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = energy.diff + 0.2, h = 2 * frequency(energy.diff), seasonal = "multiplicative",  
## 
##  Call:
##      damped = T) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9688 
## 
##   Initial states:
##     l = 0.1982 
##     b = 0 
##     s = 0.9969 1.014 0.9982 0.9951 1.0087 1.0036
##            0.996 0.9904 1.0074 1.0004 0.9947 0.9947
## 
##   sigma:  0.169
## 
##       AIC      AICc       BIC 
## -220.7570 -219.5378 -142.2225 
## 
## Error measures:
##                        ME       RMSE        MAE      MPE     MAPE
## Training set 8.266058e-05 0.03330002 0.02545316 -3.00969 13.61359
##                   MASE       ACF1
## Training set 0.5867337 -0.1694023
## 
## Forecasts:
##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## Jun 2019      0.1980570 0.1551530 0.2409610 0.1324410 0.2636730
## Jul 2019      0.1991651 0.1560211 0.2423091 0.1331820 0.2651482
## Aug 2019      0.2006640 0.1571952 0.2441327 0.1341843 0.2671436
## Sep 2019      0.2016823 0.1579930 0.2453716 0.1348653 0.2684994
## Oct 2019      0.1989802 0.1558762 0.2420842 0.1330584 0.2649021
## Nov 2019      0.1995989 0.1563609 0.2428369 0.1334721 0.2657258
## Dec 2019      0.2027405 0.1588219 0.2466591 0.1355728 0.2699081
## Jan 2020      0.1993391 0.1561573 0.2425208 0.1332983 0.2653799
## Feb 2020      0.1988948 0.1558093 0.2419804 0.1330012 0.2647885
## Mar 2020      0.1988870 0.1558032 0.2419709 0.1329959 0.2647781
## Apr 2020      0.2000065 0.1566801 0.2433328 0.1337445 0.2662684
## May 2020      0.2014201 0.1577875 0.2450527 0.1346898 0.2681505
## Jun 2020      0.1980406 0.1551400 0.2409412 0.1324298 0.2636514
## Jul 2020      0.1991491 0.1560084 0.2422899 0.1331710 0.2651272
## Aug 2020      0.2006484 0.1571828 0.2441139 0.1341735 0.2671232
## Sep 2020      0.2016671 0.1579808 0.2453534 0.1348547 0.2684795
## Oct 2020      0.1989657 0.1558646 0.2420668 0.1330482 0.2648832
## Nov 2020      0.1995848 0.1563495 0.2428201 0.1334621 0.2657075
## Dec 2020      0.2027266 0.1588107 0.2466425 0.1355630 0.2698902
## Jan 2021      0.1993258 0.1561466 0.2425051 0.1332888 0.2653629
## Feb 2021      0.1988820 0.1557988 0.2419652 0.1329920 0.2647721
## Mar 2021      0.1988746 0.1557930 0.2419563 0.1329869 0.2647624
## Apr 2021      0.1999944 0.1566701 0.2433187 0.1337356 0.2662532
## May 2021      0.2014083 0.1577777 0.2450390 0.1346810 0.2681357

Figure 18: Residual analysis plots for damped Holt-Winters’ multiplicative method

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt-Winters' multiplicative method
## Q* = 193.75, df = 7, p-value < 2.2e-16
## 
## Model df: 17.   Total lags used: 24

## 
##  Shapiro-Wilk normality test
## 
## data:  x$model$residuals
## W = 0.99289, p-value = 0.007487

Figure 18 shows residual analysis plots for the fitted model. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly.
There is still seasonality left in the residuals from this model with significant lag at lag 12 and 14 observed in the ACF plot. In addition from the ACF plot, there is still autocorrelation in the errors.
The residuals are not normal. The histogram and the Shapiro-Wilk normality test with p-value = 0.007487 suggest not normal residuals.

Overall, we can conclude that Holt-Winters’ method has not improved the residuals in terms of autocorrelation and seasonality.

State-space models

We use the created function to fit a set of possible candidate ETS models on the transformed and differenced data and display their accuracy measures to select the best out of the set.

customise(energy.diff+0.2, type="ets")

##    Model Damped     Bounds      MASE       AIC       BIC
## 1    AAA   TRUE       both 0.5894597 -218.5649 -140.0304
## 2    MAM   TRUE       both 0.5915705 -220.6736 -142.1391
## 3    MMM   TRUE       both 0.5966812 -215.7046 -137.1700
## 4    MAA   TRUE       both 0.5908081 -218.0424 -139.5079
## 5    AAA  FALSE       both 0.5955193 -204.0216 -129.8501
## 6    MAM  FALSE       both 0.5909966 -225.2689 -151.0974
## 7    MMM  FALSE       both 0.5933987 -219.8977 -145.7263
## 8    MAA  FALSE       both 0.5925894 -222.2298 -148.0583
## 9    AAA   TRUE admissible 0.5887306 -218.2297 -139.6952
## 10   MAM   TRUE admissible 0.5805912 -235.8559 -157.3214
## 11   MMM   TRUE admissible 0.5764635 -256.4136 -177.8791
## 12   MAA   TRUE admissible 0.5879057 -229.6944 -151.1599
## 13   AAA  FALSE admissible 0.5967392 -203.6162 -129.4447
## 14   MAM  FALSE admissible 0.5892372 -223.4786 -149.3071
## 15   MMM  FALSE admissible 0.5961053 -206.5394 -132.3679
## 16   MAA  FALSE admissible 0.6017450 -196.7697 -122.5982

The model that minimises AIC and BIC is ETS(M,Md,M) with the parameters lying in the admissible parameter region.

The following code implements ETS(M,Md,M) and displays residual analysis results.

ets.MMdM.diff = ets(energy.diff+0.2, model="MMM", damped = T, bounds = "admissible")
sum.res2(ets.MMdM.diff)

## ETS(M,Md,M) 
## 
## Call:
##  ets(y = energy.diff + 0.2, model = "MMM", damped = T, bounds = "admissible") 
## 
##   Smoothing parameters:
##     alpha = -0.1143 
##     beta  = 0.0157 
##     gamma = 0.008 
##     phi   = 0.8907 
## 
##   Initial states:
##     l = 0.2071 
##     b = 0.9956 
##     s = 1.066 1.0619 1.0408 0.9718 0.9917 0.9746
##            0.9767 0.9964 0.9672 0.9522 0.9943 1.0065
## 
##   sigma:  0.1661
## 
##       AIC      AICc       BIC 
## -256.4136 -255.1943 -177.8791 
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE     MAPE
## Training set 0.002557368 0.03260537 0.02500763 -1.616267 13.23798
##                   MASE       ACF1
## Training set 0.5764635 -0.1251102

Figure 19: Residual analysis plots for ETS(M,Md,M) model

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

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99464, p-value = 0.03982

Figure 19 shows residual analysis plots for the fitted model. In fact, we can observe a very similar picture to Holt-Winters’ method. We can make the following comments about diagnostic checking:

The errors are spread more or less randomly. They do not show evidence of changing variance.
There is some autocorrelation left in the residuals according to significant lags in ACF.
The errors are not normal. Long tails in the histogram and the Shapiro-Wilk normality test with p-value = 0.03982 suggest not normal residuals.

Overall, we can conclude that this model is not performing better at capturing seasonality in the series.

The auto ETS model was fitted to see what the software automatically suggested model is.

auto.fit <- ets(energy.diff)
sum.res2(auto.fit)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = energy.diff) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = -1e-04 
## 
##   sigma:  0.0333
## 
##       AIC      AICc       BIC 
## -251.2721 -251.2304 -238.1830 
## 
## Training set error measures:
##                         ME       RMSE        MAE      MPE     MAPE
## Training set -1.395225e-06 0.03326331 0.02547093 101.4309 101.4309
##                   MASE       ACF1
## Training set 0.5871433 -0.1694556

Figure 20: Residual analysis plots for ETS(A,N,N) model

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

## 
##  Shapiro-Wilk normality test
## 
## data:  x$residuals
## W = 0.99292, p-value = 0.007651

The model suggested automatically is ETS(A,N,N) which is a simple exponential smoothing method with additive errors.

The residual analysis in Figure 20 shows the following:

Residuals are spread more or less randomly.
There are some significant lags in the ACF plot. Ljung-Box test reports p-value < 0.05, therefore there is evidence of autocorrelation in the residuals. In addition, lag 12 and 24 are significant which suggests seasonality effect is still present in the residuals.
The errors are not normal. Long tails in the histogram and the Shapiro-Wilk normality test with p-value = 0.007651 suggest not normal residuals.

Overall, the residual analysis suggests the auto model is not successful at capturing autocorrelation and seasonality in energy series.

Model selection

The following code chunk is used to compare the different models estimated based on Information Criteria and MASE:

dataset <- c("Original","Original","Original","Original","Original","Original","Original","Original","Trans & Diff","Trans & Diff","Trans & Diff","Trans & Diff","Trans & Diff","Trans & Diff")
methods <- c("DLM","DLM","DLM","Holt Winters'","State-space","State-space","State-space","State-space","DLM","DLM","DLM","Holt Winters'","State-space","State-space")
models <- c("model1","model2","model3","hw.md", "ets.MAdM", "ets.MMdM","ets.MMM","fit.auto","model1.1","model2.1","model3.1","hw.md.diff", "ets.MMdM.diff", "auto.fit")

AIC <- c(AIC(model1), AIC(model2), AIC(model3),hw.md$model$aic, ets.MAdM$aic, ets.MMdM$aic, ets.MMM$aic, fit.auto$aic
         ,AIC(model1.1), AIC(model2.1), AIC(model3.1),hw.md.diff$model$aic, ets.MMdM.diff$aic, auto.fit$aic)

BIC <- c(BIC(model1), BIC(model2), BIC(model3),hw.md$model$bic, ets.MAdM$bic, ets.MMdM$bic, ets.MMM$bic, fit.auto$bic,
         BIC(model1.1), BIC(model2.1), BIC(model3.1),hw.md.diff$model$bic, ets.MMdM.diff$bic, auto.fit$bic)

MASE <- c(MASE(lm(model1))[1,1], MASE(lm(model2))[1,1], MASE(lm(model3))[1,1],accuracy(hw.md)[,6], 
          accuracy(ets.MAdM)[,6], accuracy(ets.MMdM)[,6],accuracy(ets.MMM)[,6],accuracy(fit.auto)[,6],
          MASE(lm(model1.1))[1,1], MASE(lm(model2.1))[1,1], MASE(lm(model3.1))[1,1],accuracy(hw.md.diff)[,6], 
          accuracy(ets.MMdM.diff)[,6], accuracy(auto.fit)[,6])

ALL <- data.frame(dataset,methods,models, AIC, BIC, MASE)
colnames(ALL) <- c("Dataset","Methods","Model",  "AIC", "BIC","MASE")

ALL

##         Dataset       Methods         Model        AIC        BIC
## 1      Original           DLM        model1  3058.2922  3124.0448
## 2      Original           DLM        model2  3010.5370  3080.6460
## 3      Original           DLM        model3  2986.9377  3061.3998
## 4      Original Holt Winters'         hw.md  4702.3820  4781.3155
## 5      Original   State-space      ets.MAdM  4694.2532  4773.1867
## 6      Original   State-space      ets.MMdM  4698.1578  4777.0913
## 7      Original   State-space       ets.MMM  4736.6237  4811.1720
## 8      Original   State-space      fit.auto  4696.9250  4775.8585
## 9  Trans & Diff           DLM      model1.1 -2285.0477 -2219.6282
## 10 Trans & Diff           DLM      model2.1 -2328.6710 -2258.9178
## 11 Trans & Diff           DLM      model3.1 -2322.6353 -2248.5520
## 12 Trans & Diff Holt Winters'    hw.md.diff  -220.7570  -142.2225
## 13 Trans & Diff   State-space ets.MMdM.diff  -256.4136  -177.8791
## 14 Trans & Diff   State-space      auto.fit  -251.2721  -238.1830
##         MASE
## 1  0.4230642
## 2  0.4013518
## 3  0.3928438
## 4  0.6438479
## 5  0.6426198
## 6  0.6417418
## 7  0.6690150
## 8  0.6434487
## 9  0.6487095
## 10 0.6290276
## 11 0.6278907
## 12 0.5867337
## 13 0.5764635
## 14 0.5871433

The Information Criteria and the MASE in Figure 21 shows the following:

The model with lowest MASE is the DLM with 3 lags from the original dataset, MASE = 0.3928438.
The model with the highest MASE is the ETS model with multiplicative error, trend and seasonality from the original dataset, MASE = 0.6690150.

Based on the residual analysis the state-space model with multiplicative error, trend and seasonality is the best model. Even though the IC and MASE is high for this model, we select this model for forecasting, based on the fact that this model got rid of the seasonal effect in the residuals. We think that this is more important for reliable forecasts than minimising ICs and MASE.

Forecasting

The following code chunk is used to visualize the simulated future sample paths.

A = ts(matrix(NA,10,5000), start = c(2019, 6), frequency = 12)
M = 5000
for(i in 1:M){
  A[,i]= simulate(ets.MMM,initstate = ets.MMM$states[25,], nsim = 10)
  }

plot(energy, xlim=c(1970, 2021), ylim = range(energy, A), ylab = "Industrial production (IP) index", xlab = "Year",
     main="Simulated future sample paths")

for(i in 1:4000){
  lines(A[,i], col="cornflowerblue")
}

text(1970,100,"Historical data", adj=0)
text(2012,70,"4000 simulated\n future sample\n paths",adj=0)

Figure 22: Simulated future sample paths

Figure 22 shows the distribution of predictions on our time series plot based on 4000 simulated sample paths. Based on these simulations the characteristics, such as confidence interval and mean, of the prediction distribution is estimated.

The following code chunk is used to visualize the 10 months ahead predictions of monthly energy production.

N = 10
data = energy
xlim=c(1970,2021)
Pi = array(NA, dim=c(N,2))
avrg = array(NA, N)

for (i in 1:N){
  Pi[i,] = quantile(A[i,],type=8,prob=c(.05,.95))
  avrg[i] = mean(A[i,])
}

Pi.lb = ts(Pi[,1],start=end(data),f=12)
Pi.ub = ts(Pi[,2],start=end(data),f=12)
avrg.pred = ts(avrg,start=end(data),f=12)

plot(data,xlim=xlim , ylim=range(data,A),ylab="Industrial production (IP) index",xlab="Year", 
     main="10 months ahead predictions")
lines(Pi.lb,col="blue", type="l")
lines(Pi.ub,col="blue", type="l")
lines(avrg.pred,col="red", type="l")
legend("topleft", lty=1, pch=1, col=c("black","blue","red"), text.width = 17,
       c("Data","90% confidence interval","Mean prediction"))

Figure 23: 10 months ahead predictions of monthly energy production

Figure 23 shows the 90% confidence interval and the mean predictions 10 months ahead.

Conclusion

The best model to forecast industrial production of electric and gas utilities is the state-space model with multiplicative error, trend and seasonality on the original dataset. This model is best in terms of the residual analysis. The forecast based on the prediction distribution looks good as it has captured the seasonality and the trend which the energy production is affected by.

MATH1307 Final Project

Monthly Industrial production (IP) index data foreacsting analysis

Xiyue Shu (s3705474), Anna Krinochkina (s3712761), Anne Mette Hesselholt (s3797131)

Create user functions

Introduction

Method

Data description

Data exploration and visualisation

Decomposition

Model fitting

Dynamic linear regression method

Exponential smoothing method

State-space models

Data transformation and differencing

Model fitting on transformed and differenced data

Dynamic linear regression method

Exponential smoothing method

State-space models

Model selection

Forecasting

Conclusion