Textbook: Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. Link. Accessed on September 7, 2020.

# Required R packages
library(fpp2)
library(tidyverse)
library(gridExtra)

Exercise 3.1

For the following series, find an appropriate Box-Cox transformation in order to stabilize the variance.

# Function that finds an appropriate Box-Cox transformation in order to stabilize the variance
bct = function(ts, title){
  a = autoplot(ts) + labs(title = sprintf("Before: %s", title))
  lambda = BoxCox.lambda(ts)
  at = autoplot(BoxCox(ts, lambda)) + labs(title = sprintf("Transformed: %s", title), 
                                           subtitle = sprintf("lambda = %0.2f", lambda))
  grid.arrange(a, at)
}

bct(usnetelec, title = "Annual US Net Electricity Generation")

The plot above fits with a transformation of \(\lambda = 0.52\), however, there doesn’t appear to be any noticeable changes that can help with better interpretation. This may be due to the lack of any trend of increasing or decreasing variation in the time series data.

bct(usgdp, title = "Quarterly US GDP")

The plot above fits with a transformation of \(\lambda = 0.37\), and it is apparent that the curvature of the graph was reduced.

bct(mcopper, title = "Monthly Copper Prices")

The plot above fits with a transformation of \(\lambda = 0.19\), and the transformation allowed seasonal variation to be more observably stable.

bct(enplanements, title = "Monthly US Domestic Enplanements")

The plot above fits with a transformation of \(\lambda = -0.23\), and this transformation was also able to stabilize the variation in the data.

Exercise 3.2

Why is a Box-Cox transformation unhelpful for the cangas data?

bct(cangas, title = "Monthly Canadian Gas Production")

After examining the original and transformed plots for the Monthly Canadian Gas Production time series, there is no apparent stabilization of the variation. This may be due to the increasing and decreasing variations of the time series. Note that, from 1960 to 1975 there is an increasing trend, which tapers until 1987 before it increases again. Moreover, there is evidence of seasonality variance in which the middle portion of the data exhibits widely than the lower and upper portions. Therefore, the transformation does not show to be helpful for the plot and it would be helpful if we had a larger transforming \(\lambda\) for the lower and upper portion, and a smaller transforming \(\lambda\) for the middle portion.

Exercise 3.3

What Box-Cox transformation would you select for your retail data (from Exercise 3 in Section 2.10)?

retaildata = readxl::read_excel("retail.xlsx", skip = 1)
myts = ts(retaildata[, 13], frequency = 12, start = c(1982,4))
bct(myts, title = "New South Wales - Department Stores")

For this time series subset, the optimal Box-Cox transformation is with \(\lambda = 0.21\). The resulting transformation highlights the stabilization of the seasonal variability throughout the years. The plot shows how the data before 1990, with smaller variability, were stretched, while after the year 2000, with larger variability, were diminished.

Exercise 3.8

For your retail time series (from Exercise 3 in Section 2.10):

  1. Split the data into two parts using
myts.train = window(myts, end=c(2010,12))
myts.test = window(myts, start=2011)
  1. Check that your data have been split appropriately by producing the following plot.
autoplot(myts) +
autolayer(myts.train, series="Training") +
autolayer(myts.test, series="Test")

# Ensures that data was split appropriately
tail(myts.train)
##        Jul   Aug   Sep   Oct   Nov   Dec
## 2010 499.4 419.6 465.3 470.4 542.8 887.5
head(myts.test)
##        Jan   Feb   Mar   Apr   May   Jun
## 2011 461.2 345.4 415.7 479.9 453.9 471.3
  1. Calculate forecasts using snaive applied to myts.train.
# Seasonal Naive Method
fc = snaive(myts.train)
fc
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2011          491.3 459.6747 522.9253 442.9333 539.6667
## Feb 2011          354.1 322.4747 385.7253 305.7333 402.4667
## Mar 2011          443.7 412.0747 475.3253 395.3333 492.0667
## Apr 2011          447.4 415.7747 479.0253 399.0333 495.7667
## May 2011          471.7 440.0747 503.3253 423.3333 520.0667
## Jun 2011          512.6 480.9747 544.2253 464.2333 560.9667
## Jul 2011          499.4 467.7747 531.0253 451.0333 547.7667
## Aug 2011          419.6 387.9747 451.2253 371.2333 467.9667
## Sep 2011          465.3 433.6747 496.9253 416.9333 513.6667
## Oct 2011          470.4 438.7747 502.0253 422.0333 518.7667
## Nov 2011          542.8 511.1747 574.4253 494.4333 591.1667
## Dec 2011          887.5 855.8747 919.1253 839.1333 935.8667
## Jan 2012          491.3 446.5751 536.0249 422.8991 559.7009
## Feb 2012          354.1 309.3751 398.8249 285.6991 422.5009
## Mar 2012          443.7 398.9751 488.4249 375.2991 512.1009
## Apr 2012          447.4 402.6751 492.1249 378.9991 515.8009
## May 2012          471.7 426.9751 516.4249 403.2991 540.1009
## Jun 2012          512.6 467.8751 557.3249 444.1991 581.0009
## Jul 2012          499.4 454.6751 544.1249 430.9991 567.8009
## Aug 2012          419.6 374.8751 464.3249 351.1991 488.0009
## Sep 2012          465.3 420.5751 510.0249 396.8991 533.7009
## Oct 2012          470.4 425.6751 515.1249 401.9991 538.8009
## Nov 2012          542.8 498.0751 587.5249 474.3991 611.2009
## Dec 2012          887.5 842.7751 932.2249 819.0991 955.9009
  1. Compare the accuracy of your forecasts against the actual values stored in myts.test.
accuracy(fc, myts.test)
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set  11.22432 24.67736 19.07958  3.240636 5.581014 1.000000
## Test set     -17.55417 25.70837 20.65417 -3.862206 4.506843 1.082527
##                     ACF1 Theil's U
## Training set -0.05265994        NA
## Test set     -0.31868147  0.212556

The mean error for the training set is 11.2 while the mean error for the test set is -17.6. Comparing its absolute value, both values are a bit far from zero which usually indicates its not as good of a fit. Next, the root mean square error is very similar for the training and test sets, 24.7 and 25.7 respectively. Likewise, the mean absolute error is very similar, 19.0 and 20.7 respectively. Because the RMSE and MAE are quite similar, this suggests that there is little variance within the errors and that the model does not over-fit.

Moreover, the mean percentage error for the training set is 3.24% and for the testing set, it is 3.86%. The mean absolute percentage error is 5.58% for the training set and 4.50% for the testing test. Lastly, the mean absolute scaled error for the training set is 1.00% and for the testing set, it is 1.08%.

  1. Check the residuals. Do the residuals appear to be uncorrelated and normally distributed?
checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 98.61, df = 24, p-value = 5.18e-11
## 
## Model df: 0.   Total lags used: 24

Here, there are a time plot, ACF plot, and histogram of the residuals (with an overlaid normal distribution for comparison). Firstly, the residual time plot highlights a few large positive and negative residuals after the year 2000 and 2005. Next, the histogram suggests that the residuals may be right-skewed and the mean of the residuals is far from zero. Moreover, there are a few significant correlations in the residuals series, suggesting the forecasts may not be as good. And lastly, based on the Ljung-Box test, the residuals are distinguishable from a white noise series, \(Q^* = 98.61, df = 24, p-value < 0.05\). Overall, it seems that forecasts from this method will probably be not as good.

  1. How sensitive are the accuracy measures to the training/test split?

To understand the sensitivity of the accuracy measure to the training/test split, the following performs multiple accuracy tests on different splits for comparison.

ts_accuracy = function(ts, year){
  train = window(ts, end = c(year, 12))
  test = window(ts, start = year + 1)
  accuracy = data.frame(accuracy(snaive(train), test))
  row.names(accuracy) = c(sprintf("Training Set, end %d", year), sprintf("Test Set, start %d", year + 1))
  return(accuracy)
}

sensitivity_test = data.frame()
for (year in c(2002:2007)){
  test = ts_accuracy(myts, year)
  sensitivity_test = rbind(sensitivity_test, test)
}

sensitivity_test
##                               ME     RMSE      MAE        MPE     MAPE
## Training Set, end 2002 11.871308 24.04889 18.00970  3.7286137 5.921119
## Test Set, start 2003   34.029167 41.34762 34.50417  7.6403375 7.752129
## Training Set, end 2003 12.350602 24.23704 18.23896  3.7904127 5.888030
## Test Set, start 2004   23.991667 34.39360 27.97500  5.4203181 6.426601
## Training Set, end 2004 12.905747 24.42501 18.52337  3.8713932 5.872569
## Test Set, start 2005   -2.125000 22.52676 16.90833 -0.5240122 3.911903
## Training Set, end 2005 12.300366 24.51798 18.62051  3.6913109 5.824473
## Test Set, start 2006    7.854167 28.79957 22.25417  1.2306722 4.814107
## Training Set, end 2006 11.676491 24.71612 18.71298  3.4963261 5.783780
## Test Set, start 2007   24.533333 32.82566 27.94167  4.8395878 5.612420
## Training Set, end 2007 12.042761 24.79339 18.91279  3.5219066 5.744110
## Test Set, start 2008   12.500000 27.15249 20.66667  2.4951476 4.289719
##                             MASE        ACF1 Theil.s.U
## Training Set, end 2002 1.0000000 -0.13169894        NA
## Test Set, start 2003   1.9158652 -0.07862768 0.4156029
## Training Set, end 2003 1.0000000 -0.13513269        NA
## Test Set, start 2004   1.5338049  0.04631379 0.3333066
## Training Set, end 2004 1.0000000 -0.10596807        NA
## Test Set, start 2005   0.9128108 -0.09885889 0.2315329
## Training Set, end 2005 1.0000000 -0.09270024        NA
## Test Set, start 2006   1.1951425 -0.05140473 0.2825023
## Training Set, end 2006 1.0000000 -0.09111525        NA
## Test Set, start 2007   1.4931701 -0.29007084 0.2763369
## Training Set, end 2007 1.0000000 -0.08988359        NA
## Test Set, start 2008   1.0927347  0.03754183 0.2032548

The accuracy test on different splits shows that there are variations in how well the model fits. For instance, the better forecast seems to be if the split was done ending in December 2004 to make a forecast for the following year. The RSME and MAE are 22.5 and 16.9, respectively, and are quite low compared to the values of the other splits. Moreover, the percentage errors are the lowest than the other splits. In conclusion, accuracy measures are very sensitive to the split.