Forecasting Principles and Practice

Load required packages

packages <- c("tidyverse", "fpp2", "forecast", "kableExtra", "broom", "ggplot2", "caret", "e1071", "knitr", "GGally", "VIM", "mlbench", "car", "corrplot", "mice", "seasonal", "fma", "latex2exp","gridExtra")
pacman::p_load(char = packages)

3.1

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

  • usnetelec
  • usgdp
  • mcopper
  • enplanements

Answer

A good value of λ is one which makes the size of the seasonal variation about the same across the whole series, as that makes the forecasting model simpler.

Let’s use BoxCox.lambda function to pick up the best lambda value for appropriate box cox transformation.

BoxCox has a built in function to find a good lambda value to stabilize the variance.

Step 1) Choose the best value of lambda via BoxCox.lambda

Step 2) Apply Box Cox on the time series data,

usnetelec : Annual US net electricity generation (billion kwh) for 1949-2003

## wt - without transformation

usnetelec_wt <- autoplot(usnetelec)
usnetelec_transformed <- BoxCox.lambda(usnetelec)
usn_BoxCoxgraph <-autoplot(BoxCox(usnetelec, usnetelec_transformed))
grid.arrange( grobs = list(usnetelec_wt,usn_BoxCoxgraph))

Inference: In the above graph, there is not much seasonal variation in the data to be smoothed over, the values have been normalized by the Box-Cox transformation so that may help in prediction if other values are as small.

usgdp: Quarterly US GDP. 1947:1 - 2006.1.

## wt - without transformation

usg_wt<- autoplot(usgdp)
usgdp_transformed <- BoxCox.lambda(usgdp)
usg_BoxCoxgraph<- autoplot(BoxCox(usgdp, usgdp_transformed))
grid.arrange( grobs = list(usg_wt,usg_BoxCoxgraph))

Inference: In the above graph, which is similar to usnetelec, there is not much seasonal variation in the data to be smoothed over, also the values have been normalized by the Box-Cox transformation so that may help in prediction if other values are as small.

mcopper: Monthly copper prices. Copper, grade A, electrolytic wire bars/cathodes,LME,cash (pounds/ton) Source: UNCTAD (http://stats.unctad.org/Handbook).

## wt - without transformation

mc_wt <- autoplot(mcopper)
mcopper_transformed <- BoxCox.lambda(mcopper)
mc_BoxCoxgraph <-autoplot(BoxCox(mcopper, mcopper_transformed))
grid.arrange(mc_wt, mc_BoxCoxgraph)

Inference: In the above graph, there seems to be a shift in the mcopper dataset, in the Box-Cox transformation graph due to the spike in copper prices.

enplanements: "Domestic Revenue Enplanements (millions): 1996-2000. SOURCE: Department of Transportation, Bureau of Transportation Statistics, Air Carrier Traffic Statistic Monthly.

## wt - without transformation

enp_wt<- autoplot(enplanements)
enplanements_transformed <- BoxCox.lambda(enplanements)
enp_BoxCoxgraph <- autoplot(BoxCox(enplanements, enplanements_transformed))
grid.arrange( grobs = list(enp_wt,enp_BoxCoxgraph))

Inference: In the above graph, there is not much seasonal variance in the data to be smoothed, needed values have been normalized by box cox transformation.

3.2

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

Answer

Let’s explore cangas data:

Cangas is a time series data. It is the monthly Canadian gas production, billions of cubic metres, January 1960 - February 2005

## wt - without transformation

cangas_transformed <- BoxCox.lambda(cangas)
cangas_wt<- autoplot(cangas) 
cangas_BoxCoxgraph<- autoplot(BoxCox(cangas, cangas_transformed))
grid.arrange( grobs = list(cangas_wt,cangas_BoxCoxgraph))

Inference:

  • As BoxCox.lambda picks the best value of lambda and after applying box cox on the time series data, we were not able to stabilize the variance.

  • As we can see the comparison of both autoplots, the box-cox transformation is generating similar results to the original, in other words the transformation does not provide stationarity to time series.

  • Probably as the variation in the cangas data set is not extreme enough for a B box-Cox transformation to have any real effect. We saw this behaviour with usnetelec and usgdp in the previous exercise, and even enplanements to a certain extent.

3.3

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

Answer

Let’s replicate our answer from Exercise 3 in Section 2.10:

  1. Read the data into R with the following script:
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
#View(retaildata)
  1. Select one of the time series as follows
new_series <- ts(retaildata[,"A3349873A"],
  frequency=12, start=c(1982,4))
  1. Explore your chosen retail time series using the following functions: autoplot(), ggseasonplot(), ggsubseriesplot(), gglagplot(), ggAcf(). Can you spot any seasonality, cyclicity and trend? What do you learn about the series?
autoplot(new_series)

## wt - without transformation

retaildat_wt <- BoxCox.lambda(new_series)
autoplot(BoxCox(new_series, retaildat_wt))

Inference:

  • From the above transformed plot, we see that transformation has smoothed out the data to a greater extent. The seasonal variations are stabilizied and time series is almost stationary in the transformed plot.

  • The seasonal peaks are more stationary than the original timeseries data plot.

3.8

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

a. Split the data into two parts using

new_series.train <- window(new_series, end=c(2010,12))
new_series.test <- window(new_series, start=2011)

b. Check that your data have been split appropriately by producing the following plot.</strong>

autoplot(new_series) +
  autolayer(new_series.train, series="Training") +
  autolayer(new_series.test, series="Test")

Yes, split is correct in the above plot.

Let’s confirm it via head() and tail() function

tail(new_series.train)
##        Jul   Aug   Sep   Oct   Nov   Dec
## 2010 268.5 277.0 278.7 279.0 319.3 400.2
head(new_series.test)
##        Jan   Feb   Mar   Apr   May   Jun
## 2011 296.2 302.5 310.8 274.8 267.0 263.8

c. Calculate forecasts using snaive applied to new_series.train.

fc <- snaive(new_series.train)

fc
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2011          266.2 240.2540 292.1460 226.5190 305.8810
## Feb 2011          240.0 214.0540 265.9460 200.3190 279.6810
## Mar 2011          267.5 241.5540 293.4460 227.8190 307.1810
## Apr 2011          260.7 234.7540 286.6460 221.0190 300.3810
## May 2011          272.8 246.8540 298.7460 233.1190 312.4810
## Jun 2011          260.5 234.5540 286.4460 220.8190 300.1810
## Jul 2011          268.5 242.5540 294.4460 228.8190 308.1810
## Aug 2011          277.0 251.0540 302.9460 237.3190 316.6810
## Sep 2011          278.7 252.7540 304.6460 239.0190 318.3810
## Oct 2011          279.0 253.0540 304.9460 239.3190 318.6810
## Nov 2011          319.3 293.3540 345.2460 279.6190 358.9810
## Dec 2011          400.2 374.2540 426.1460 360.5190 439.8810
## Jan 2012          266.2 229.5068 302.8932 210.0826 322.3174
## Feb 2012          240.0 203.3068 276.6932 183.8826 296.1174
## Mar 2012          267.5 230.8068 304.1932 211.3826 323.6174
## Apr 2012          260.7 224.0068 297.3932 204.5826 316.8174
## May 2012          272.8 236.1068 309.4932 216.6826 328.9174
## Jun 2012          260.5 223.8068 297.1932 204.3826 316.6174
## Jul 2012          268.5 231.8068 305.1932 212.3826 324.6174
## Aug 2012          277.0 240.3068 313.6932 220.8826 333.1174
## Sep 2012          278.7 242.0068 315.3932 222.5826 334.8174
## Oct 2012          279.0 242.3068 315.6932 222.8826 335.1174
## Nov 2012          319.3 282.6068 355.9932 263.1826 375.4174
## Dec 2012          400.2 363.5068 436.8932 344.0826 456.3174

d. Compare the accuracy of your forecasts against the actual values stored in new_series.test.

accuracy(fc,new_series.test)
##                     ME     RMSE      MAE       MPE      MAPE     MASE      ACF1
## Training set  7.772973 20.24576 15.95676  4.702754  8.109777 1.000000 0.7385090
## Test set     55.300000 71.44309 55.78333 14.900996 15.082019 3.495907 0.5315239
##              Theil's U
## Training set        NA
## Test set      1.297866

Inference:

  • The seasonal naive forecasts accuracy in the test set is significantly lower than the training set.

  • Because the root mean square error for test set is 71.443 and the mean absolute error for test set is about 55.783.

  • Now as we see that the time series data from Jan 2011 to Dec 2012 ranges between 269 and 400, and accuracy scores tells us that it is relatively large error.

  • Also, the mean absolute percentage error for test set is 15.082%, by which we can infer that on average, forecast is 15.082% off.

e. Check the residuals 

checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 624.45, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

Do the residuals appear to be uncorrelated and normally distributed?

Inference on Ljung Box Test

Hypotheses

  • The null hypothesis of the Box Ljung Test, H0, is that our model does not show lack of fit (or in simple terms—the model is just fine). The alternate hypothesis, Ha, is just that the model does show a lack of fit.

  • In this case there is a significant p-value in this test which rejects the null hypothesis that the time series isn’t autocorrelated

Information on Residuals:

Residuals are useful in checking whether a model has adequately captured the information in the data. A good forecasting method will yield residuals with the following properties:

  • The residuals are uncorrelated. If there are correlations between residuals, then there is information left in the residuals which should be used in computing forecasts.

  • The residuals have zero mean. If the residuals have a mean other than zero, then the forecasts are biased.

Also, In addition to these essential properties, it is useful (but not necessary) for the residuals to also have the following two properties.

  • The residuals have constant variance.
  • The residuals are normally distributed.

Inference on Plots:

  • From the plots, it is apparent that the residuals are highly correlated. This means that there is information left in the residuals that can be used to improve the forecast.

  • The mean of the residuals is not centered around zero, which means the forecast is biased. Although the residuals appear to be normally distributed, the variance is not constant. Therefore the prediction interval is not reliable.

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

Answer

To test how sensitive are the accuracy measures to the training/test split, let’s perform forecasts multiple times, each time using different year to split the data, and compate the accuracies calculated below.

sen <- function(split_year){
  trainset <- window(new_series, end=c(split_year, 12))
  testset <- window(new_series, start=split_year+1)
  acc <- accuracy(snaive(trainset), testset)
  return(acc)
}

splits <- c(2000:2011)

accs <- data.frame()
for (year in splits){
  acc <- sen(year)
  temp <- data.frame(t(acc[2,c(1:6)]))
  accs <- rbind(accs, temp)
}
row.names(accs) <- splits
accs
##               ME      RMSE       MAE        MPE      MAPE      MASE
## 2000  46.9958333 51.253085 46.995833 16.0619584 16.061958 3.1627528
## 2001  20.2958333 23.203403 20.454167  6.9090812  6.969629 1.2630187
## 2002   0.9541667 15.196395 12.495833  0.3421824  4.229447 0.7756915
## 2003 -23.9208333 28.733930 25.020833 -8.4660641  8.865573 1.5670274
## 2004  -0.6875000  7.066087  5.820833 -0.2190730  2.103361 0.3571399
## 2005  10.3375000 19.096673 13.870833  3.1513972  4.490635 0.8747372
## 2006  11.5500000 18.573794 13.708333  3.4641276  4.285225 0.8820923
## 2007   4.9208333 24.992541 20.020833  1.9475253  6.575480 1.2948451
## 2008  -0.2375000 25.096256 18.470833  0.1008644  6.074517 1.1985484
## 2009  -5.6958333 34.321622 27.137500 -3.3305299  8.948578 1.7364971
## 2010  55.3000000 71.443089 55.783333 14.9009957 15.082019 3.4959067
## 2011  65.4291667 78.693755 67.345833 15.8536151 16.520064 4.0179004

  • From the above table, tt is apparent that the accuracy measures are very sensitive to the split.

  • For example, if we see lowest MAPE value i.e. 2.1%, it came through data from Apr 1982 to Dec 2004 taken as train set and data from Jan 2005 to Dec 2006 taken as test set. the accuracy is be very good which also means that on average the forecast is just 2.1% off.

  • On the other hand, if we pick Dec 2011 to do the split, the accuracy will be very bad, with a MAPE of 15.1% which also means that on average the forecast is 15.1% off.