Intro

This week’s discussion topic is comparing ETS & ARIMA models to a GARCH model. ARCH (autoregressive conditionally heteroscedastic) model is a model for the variance of a time series. ARCH models are used to describe a changing, possibly volatile variance. Although an ARCH model could possibly be used to describe a gradually increasing variance over time, most often it is used in situations in which there may be short periods of increased variation. A GARCH (generalized autoregressive conditionally heteroscedastic) model uses values of the past squared observations and past variances to model the variance at time t. (11.1 ARCH/GARCH Models) The models I will build will be for Electricity net generation measured in billions of kilowatt hours (kWh). The data is from 1973-01-01, to 2013-06-01. This data comes from the fpp2 package.

Visualizations

This season has an increasing trend, and elements of seasonality. Additionally, it does not appear so much that the variance is increasing over time, but the increasing trend is causing the peaks to be higher, and the troughs to be lower. Let’s decompose this series, and see the elects separate from one another.

autoplot(decompose(usmelec, "a")) + 
  ggtitle("Additive Decomposition")

autoplot(decompose(usmelec, "m")) + 
  ggtitle("Multiplicative Decomposition")

library(seasonal)
usmelec %>% 
 stl(t.window=13, s.window="periodic", robust=TRUE) %>%
  autoplot() +
  ggtitle("STL Decomposition")

usmelec %>% 
  seas(x11 = "") %>%
  autoplot() +
  ggtitle("X11 Decomposition")

This data is very seasonal. The X11 decomposition seems to do the best job removing all of the components from the remainder series. It is very stationary, and follows no patter, unlike the other 3.

Models

The three models we will be building and comparing are ETS, ARIMA, and GARCH. First, let’s build out the ETS models, and select the best performing one. we’ll need to split the data into a training set and a test set. We will use until the end of 2004 to train, and 2005: 2013-06 to test.

ETS

The ETS models I will test are “ZZZ”, which is auto selected, “ANN” is simple exponential smoothing with additive errors, “MNN” is simple exponential smoothing with multiplicative errors, “MAM” is multiplicative Holt-Winters’ method with multiplicative errors, and “MAN” Holt’s linear method with multiplicative errors.

train <- window(usmelec, end = c(2004,12))

test<- window(usmelec, start = c(2005,1))

library(kableExtra)
ets_model_types <- c("ZZZ", "ANN", "MNN", "MAM", "MAN")

ts_model_list <- lapply(setNames(unique(ets_model_types), 
                                  unique(ets_model_types)),
                         
                         function(x){
                           
                           model <- ets(y = train, model = x)
                           
                           fcsts <- forecast(model, h = 32)
                           
                           eval_metrics <-
                             accuracy(object = fcsts,
                                      x = test)
                           return(eval_metrics)
                           
                         })

kable(ts_model_list[["ZZZ"]], caption = "ZZZ")
ZZZ
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.5738632 6.694844 4.975248 0.1781215 2.059266 0.5794835 0.0159539 NA
Test set 4.8199138 14.973278 11.588235 1.0866991 3.228432 1.3497200 0.4636921 0.4532237 
kable(ts_model_list[["ANN"]], caption = "ANN")
ANN
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.5157664 21.16762 17.35242 -0.1517985 7.159502 2.021094 0.1296537 NA
Test set 2.2053259 36.81947 29.98589 -0.4386358 8.565380 3.492556 0.5647535 1.145587 
kable(ts_model_list[["MNN"]], caption = "MNN")
MNN
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.4761309 21.02734 17.14512 -0.1458328 7.083829 1.996949 0.0791501 NA
Test set -0.4294072 36.75587 30.47991 -1.2185408 8.777484 3.550095 0.5647535 1.153688 
kable(ts_model_list[["MAM"]], caption = "MAM")
MAM
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.5738632 6.694844 4.975248 0.1781215 2.059266 0.5794835 0.0159539 NA
Test set 4.8199138 14.973278 11.588235 1.0866991 3.228432 1.3497200 0.4636921 0.4532237 
kable(ts_model_list[["MAN"]], caption = "MAN")
MAN
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 0.4599866 21.03587 17.16474 -0.1567924 7.096769 1.999234 0.0781813 NA
Test set -0.4390869 36.75533 30.48178 -1.2213853 8.778278 3.550314 0.5647436 1.153708

It appears the auto selection found “MAM”. It significantly out performs the others, so we shall use it to test. Let’s fit the model.

my_ets <- ets(y = train, model = "MAM")

ARIMA

For the ARIMA models, I will derive two \(ARIMA(p,d,q)\) orders automatically. The first will use the functionality of forecast::auto.arima to account for seasonality directly by estimating the \((PDQ)\) component of the model is as such

$$

ARIMA(p,d,q)(P,Q,D)[m]

$$

The second will use forecast::stlf(method = "arima"). This function will forecast STL objects are obtained by applying a non-seasonal forecasting method to the seasonally adjusted data and re-seasonalizing using the last year of the seasonal component.

First, let’s look at some features of this series that will give use a sense of its structure from an “ARIMA” point of view. We will determine if our series has a unit root, and needs to be differenced. I will use the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test to stay consistent w/ the method used in the text. In this test, \(Ho\) is that the data are stationary, and we look for evidence that the null hypothesis is false. Consequently, rejections of \(Ho\) suggest that differencing is required.

library(urca)

diffs <- usmelec %>% 
  ndiffs()

print(paste("The number of differences this series requires is", diffs))
## [1] "The number of differences this series requires is 1"
sdiffs<- usmelec %>% 
  diff(lag=12) %>% 
  ndiffs()

print(paste("The number of seasonal differences this series requires is", diffs))
## [1] "The number of seasonal differences this series requires is 1"

The results of the tests indicate we should do both a seasonal difference and a first difference in our ARIMA models.

Now let’s take a look at the ACF and PACF plots, and fit the models.

train %>% 
  diff %>% 
  ggAcf()

train %>%
  diff %>% 
  ggPacf()

my_auto_arima <- auto.arima(train, seasonal = TRUE)

my_stlf_arima <- stlf(train, method = "arima", h = length(test))

Here are the results using auto.arima()

my_auto_arima %>% 
  summary()
## Series: train 
## ARIMA(1,0,2)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1      ma2     sma1   drift
##       0.9501  -0.4299  -0.2637  -0.6893  0.4796
## s.e.  0.0261   0.0588   0.0582   0.0362  0.0584
## 
## sigma^2 estimated as 47.83:  log likelihood=-1248.55
## AIC=2509.11   AICc=2509.34   BIC=2532.62
## 
## Training set error measures:
##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.01886106 6.761253 5.066625 -0.1352361 2.067803 0.5901265
##                     ACF1
## Training set 0.006801918

And here is using STLF:

my_stlf_arima$model
## Series: x 
## ARIMA(0,1,2) with drift 
## 
## Coefficients:
##           ma1      ma2   drift
##       -0.4890  -0.2396  0.4791
## s.e.   0.0507   0.0517  0.0798
## 
## sigma^2 estimated as 32.76:  log likelihood=-1210.44
## AIC=2428.87   AICc=2428.98   BIC=2444.66
accuracy(my_stlf_arima)
##                      ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.01574581 5.693439 4.327251 -0.04961887 1.786792 0.5040092
##                      ACF1
## Training set -0.001212224

The results from the STLF make quite a bit of sense, since the ACF is sinusoidal, and the PACF has about an order of 2. This model performs best in terms of AICc, RMSE, MAPE, and MASE. I will use this model to forecast test data.

GARCH

Finally, we need to fit the GARCH model for our training set. I will fit the model to the order of \[GARCH(1,1)\]

library(fGarch)
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following objects are masked from 'package:seasonal':
## 
##     outlier, series
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
## 
##     volatility
garch_11=garchFit(data=train, formula = ~garch(1,1), cond.dist="QMLE", trace=FALSE)

summary(garch_11)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = train, cond.dist = "QMLE", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000001dad8f78>
##  [data = train]
## 
## Conditional Distribution:
##  QMLE 
## 
## Coefficient(s):
##        mu      omega     alpha1      beta1  
## 197.72706   10.15699    0.15775    0.84475  
## 
## Std. Errors:
##  robust 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     197.72706     4.48756   44.061  < 2e-16 ***
## omega   10.15699     7.92538    1.282      0.2    
## alpha1   0.15775     0.03275    4.817 1.46e-06 ***
## beta1    0.84475     0.03126   27.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -2021.1    normalized:  -5.263281 
## 
## Description:
##  Thu Jul 30 13:35:58 2020 by user: jawilliams 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  20.42772  3.665868e-05
##  Shapiro-Wilk Test  R    W      0.9399884 2.447251e-11
##  Ljung-Box Test     R    Q(10)  1516.179  0           
##  Ljung-Box Test     R    Q(15)  2479.345  0           
##  Ljung-Box Test     R    Q(20)  3211.581  0           
##  Ljung-Box Test     R^2  Q(10)  212.9904  0           
##  Ljung-Box Test     R^2  Q(15)  613.0244  0           
##  Ljung-Box Test     R^2  Q(20)  682.6537  0           
##  LM Arch Test       R    TR^2   262.6909  0           
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 10.54740 10.58855 10.54718 10.56372

Test

Let’s evaluate our models performance on the test data.

fcst_ets <- forecast(my_ets, h = length(test))

fcst_arima <- my_stlf_arima$mean

fcst_garch <- fGarch::predict(garch_11, n.ahead = length(test))

fcst_garch <- ts(fcst_garch$meanForecast, start = c(2005,1), frequency = 12)

# MAM
accuracy(fcst_ets$mean, test)
##               ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 4.22075 13.79647 10.50331 0.9615199 2.967047 0.5847733 0.4153442
#ARIMA
accuracy(fcst_arima, test)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -19.95926 27.12156 22.22037 -6.194206 6.775985 0.7973904 0.8887439
#GARCH 
accuracy(fcst_garch, test)
##                ME     RMSE      MAE      MPE     MAPE      ACF1 Theil's U
## Test set 141.3629 145.5129 141.3629 41.11859 41.11859 0.5583063  4.522462

Based on the test set evaluation metrics, the “MAM” model greatly outperforms the other two.

Let’s take a look at the forecasts versus the observed values in the test set.

autoplot(test)+ 
  autolayer(fcst_ets$mean, series = "MAM")+
  ggtitle("MAM Forecasts v Observed Values")+
  ylab("kwH")

autoplot(test)+ 
  autolayer(fcst_arima, series = "ARIMA")+
  ggtitle("ARIMA(0,1,2) with drift Forecasts v Observed Values")+
  ylab("kwH")

autoplot(test)+ 
  autolayer(fcst_garch, series = "GARCH")+
  ggtitle("GARCH(1,1)) Forecasts v Observed Values")+
  ylab("kwH")

So it appears that I completely missed the mark w/ my Garch model. The ’MAM" model performed the best, just slightly under forecasting the values in the test set.