Exercise 3.1

a. usnetelec

Box-Cox transformation: \(\lambda = 0.52\)

myts <- usnetelec
(lambda <- BoxCox.lambda(myts))
## [1] 0.5167714
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
    labs(title = "Annual US net electricity generation, 1949-2003", 
         subtitle = "Billion kWh", 
         x = "Year", y = "")

b. usgdp

Box-Cox transformation: \(\lambda = 0.37\)

myts <- usgdp
(lambda <- BoxCox.lambda(myts))
## [1] 0.366352
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
    labs(title = "Quarterly US GDP", 
         subtitle = "1947/Q1 - 2006/Q1", 
         x = "Year", y = "")

c. mcopper

Box-Cox transformation: \(\lambda = 0.19\)

myts <- mcopper
(lambda <- BoxCox.lambda(myts))
## [1] 0.1919047
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
    labs(title = "Monthly copper prices, 1960/1 - 2006/12", 
         subtitle = "Pounds/ton",
         x = "Year", y = "")

d. enplanements

Box-Cox transformation: \(\lambda = -0.23\)

myts <- enplanements
(lambda <- BoxCox.lambda(myts))
## [1] -0.2269461
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
    labs(title = "Monthly US domestic enplanements, 1996-2000", 
         subtitle = "Millions",
         x = "Year", y = "")

Exercise 3.2

Box-Cox transformation: \(\lambda = 0.58\)

The Box-Cox transformation isn’t as useful here, because the seasonal variation doesn’t increase or decrease monotonically and instead follows a complex pattern. For instance, the seasonal variation at first increases slowly (from 1960 to early 1970s), then increases rapidly (from early 1970s to late 1980s), and finally decreases slowly (from late 1980s onward). As a result the transformed time series continues to exhibit a seasonal pattern with a non-constant variance.

myts <- cangas
(lambda <- BoxCox.lambda(myts))
## [1] 0.5767759
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
    labs(title = "Monthly Canadian gas production, 1960/1 - 2005/2", 
         subtitle = "Billions of cubic meters",
         x = "Year", y = "")

Exercise 3.3

For the Australian retail data, I choose the 7th variable column in the retail dataset, which relates to “Turnover - New South Wales - Hardware, building and garden supplies retailing”.

The Box-Cox function computes a value of \(\lambda = 0.92\); however this is close enough to 1, so we choose 1 (i.e., no transformation) for ease of interpretation.

retaildata <- readxl::read_excel("../Week2/retail.xlsx", skip=1)
colID <- colnames(retaildata)[8]
myts <- ts(retaildata[ , colID], frequency=12, start=c(1982,4))

(lambda <- BoxCox.lambda(myts))
## [1] 0.9165544
df <- cbind(Raw = myts, Transformed = BoxCox(myts, lambda)) 
autoplot(df, facet=TRUE) + 
  labs(title = "Turnover - NSW - Hardware, building & garden supplies retailing", 
       subtitle = "1982/04 - 2013/12", 
       x = "Month", y = "")

Exercise 3.8

a. Split data into training and test sets

The training set consists of observations through 2010/12, while the test set consists of all observations thereafter.

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

b. Plot the data

The dashed line at 2011/01 divides the time series into the training set (prior to 2011) and the test set (from and after 2011).

autoplot(myts) +
  autolayer(myts.train, series="Training") +
  autolayer(myts.test, series="Test") + 
  geom_vline(xintercept = 2011, lty = 2) + 
  labs(title = "Turnover - NSW - Hardware, building & garden supplies retailing", 
       subtitle = "1982/04 - 2013/12", 
       x = "Month", y = "") + 
  guides(colour=guide_legend(title="Series"))

c. Calculate forecasts on training data

We apply the snaive function to the myts.train series to generate seasonal naive forecasts. Note that by default, the forecast horizon of snaive is 2 years (from 2011/1 to 2012/12). The plots below show the forecasts against the test data, with the forecast window highlighted.

fc <- snaive(myts.train)
summary(fc)
## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = myts.train) 
## 
## Residual sd: 24.5845 
## 
## Error measures:
##                    ME     RMSE      MAE     MPE     MAPE MASE      ACF1
## Training set 9.460661 26.30758 21.23363 4.65569 12.76289    1 0.8070166
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2011          309.2 275.4855 342.9145 257.6381 360.7619
## Feb 2011          272.6 238.8855 306.3145 221.0381 324.1619
## Mar 2011          311.1 277.3855 344.8145 259.5381 362.6619
## Apr 2011          298.2 264.4855 331.9145 246.6381 349.7619
## May 2011          313.1 279.3855 346.8145 261.5381 364.6619
## Jun 2011          305.8 272.0855 339.5145 254.2381 357.3619
## Jul 2011          307.3 273.5855 341.0145 255.7381 358.8619
## Aug 2011          330.9 297.1855 364.6145 279.3381 382.4619
## Sep 2011          362.8 329.0855 396.5145 311.2381 414.3619
## Oct 2011          361.7 327.9855 395.4145 310.1381 413.2619
## Nov 2011          364.2 330.4855 397.9145 312.6381 415.7619
## Dec 2011          395.4 361.6855 429.1145 343.8381 446.9619
## Jan 2012          309.2 261.5205 356.8795 236.2804 382.1196
## Feb 2012          272.6 224.9205 320.2795 199.6804 345.5196
## Mar 2012          311.1 263.4205 358.7795 238.1804 384.0196
## Apr 2012          298.2 250.5205 345.8795 225.2804 371.1196
## May 2012          313.1 265.4205 360.7795 240.1804 386.0196
## Jun 2012          305.8 258.1205 353.4795 232.8804 378.7196
## Jul 2012          307.3 259.6205 354.9795 234.3804 380.2196
## Aug 2012          330.9 283.2205 378.5795 257.9804 403.8196
## Sep 2012          362.8 315.1205 410.4795 289.8804 435.7196
## Oct 2012          361.7 314.0205 409.3795 288.7804 434.6196
## Nov 2012          364.2 316.5205 411.8795 291.2804 437.1196
## Dec 2012          395.4 347.7205 443.0795 322.4804 468.3196
autoplot(myts) + 
  autolayer(myts.train, series = "Training") + 
  autolayer(myts.test, series = "Test") + 
  autolayer(fc, series = "Forecast", PI = FALSE) + 
  geom_vline(xintercept = c(2011, 2013), lty = 2) + 
  labs(title = "Turnover - NSW - Hardware, building & garden supplies retailing", 
       subtitle = "1982/04 - 2013/12", 
       x = "Month", y = "") + 
  guides(colour=guide_legend(title="Series"))

Let’s take a closer look at the forecast window.

# close-up of 24-month forecast period
autoplot(myts) + 
  autolayer(myts.train, series = "Training") + 
  autolayer(myts.test, series = "Test") + 
  autolayer(fc, series = "Forecast", PI = FALSE) + 
  geom_vline(xintercept = c(2011, 2013), lty = 2) + 
  labs(title = "Turnover - NSW - Hardware, building & garden supplies retailing", 
       subtitle = "Close-up: 2009/01 - 2013/12", 
       x = "Month", y = "") + 
  guides(colour=guide_legend(title="Series")) + 
  coord_cartesian(xlim=c(2009, 2014), ylim=c(250, 425))

d. Evaluate accuracy on test data

We compute accuracy metrics on the test data, by comparing the forecasts against actual values in myts.test. It appears that the forecast accuracy measured on the test data is better than that on the training data, on the basis of several accuracy measures. For instance, the RMSE is 21.3 and MAPE is 4.8 on the test set, compared to 26.3 and 12.8 on the training set, respectively.

accuracy(fc, myts.test) %>% kable(digits = 2)
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Training set 9.46 26.31 21.23 4.66 12.76 1.00 0.81 NA
Test set 17.21 21.26 17.40 4.75 4.81 0.82 0.48 0.69

e. Check the residuals

For the seasonal naive forecast model, it is evident from the plots below that:

  • The residuals are strongly correlated as can be seen in the time series chart and the Ljung-Box test result. In fact, from the ACF chart, we see that autocorrelation is statistically significant for up to 12-month lags.

  • The residuals have non-constant variance over time, as we can see in the time series chart.

  • The residuals appear approximately normally distributed, although the distribution shown in the histogram appears slightly skewed and the left tail exhibits more outliers than expected for a normal distribution.

  • The residuals have a non-zero mean, which indicates that the forecasts are biased. In fact, the mean of the residuals is 9.5, so the model is under-predicting actuals on average. A simple improvement to the model would be to add this value to each forecast.

checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 856.11, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
mean(residuals(fc), na.rm = TRUE)
## [1] 9.460661

f. Assess sensitivity of accuracy measures

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

To answer the question, let’s repeat the training/test split multiple times and record the accuracy measures. Since snaive produces 24-month forecasts by default, and we have data from 1982 through 2013, let’s vary the date of the training/test split over a 10-year window from 2002/01 through 2011/12.

# 120-month window starting in 2002/01
start.yr <- 2002
months <- 120

# initialize training/test performance
perf.train <- NULL
perf.test <- NULL

# loop through 120 months of training/test splits
for (j in 1:months) {
  temp.train <- window(myts, end = c(start.yr, j))
  temp.test <- window(myts, start = c(start.yr, j + 1))
  temp.fc <- snaive(temp.train)
  temp.acc <- accuracy(temp.fc, temp.test)
  perf.train <- rbind(perf.train, temp.acc[1, 2:5])
  perf.test <- rbind(perf.test, temp.acc[2, 2:5])
}

# create dataframe for ggplot
cbind(Time = 1:months, Series = "Training", as_tibble(perf.train)) %>%
  pivot_longer(cols = RMSE:MAPE, names_to = "Metric", values_to = "Value") ->
  df.train
cbind(Time = 1:months, Series = "Test", as_tibble(perf.test)) %>%
  pivot_longer(cols = RMSE:MAPE, names_to = "Metric", values_to = "Value") ->
  df.test
df <- rbind(df.train, df.test)

# plot perf metrics
ggplot(df, aes(x = 2002 + (Time-1)/12, y = Value, col = Series)) + 
  geom_line() + 
  geom_vline(xintercept=2011, lty = 2) + 
  facet_wrap(~ Metric, scales = "free") + 
  labs(title = "Performance metrics for training/test splits from 2002/01 to 2011/12", 
       subtitle = "Dashed vertical line = 2011", 
       x = "Month", y = "") + 
  guides(colour=guide_legend(title="Dataset")) + 
  theme(legend.position = "bottom")

The plots demonstrate that the accuracy measures evaluated on the test set are highly sensitive to the choice of training/test split. For instance, for training/test splits ranging from 2002/01 to 2011/12 (while holding the length of the forecast window constant at 24 months), the RMSE varies between 20 and 50+ and the MAPE varies between 4% and 16%+. This sensitivity is explained by the 24-month forecast window on the test data, during which time the actuals may vary substantially from the prior 12-month period used for the seasonal naive model. On the other hand, the accuracy measures evaluated on the training set are not very sensitive to the training/test split, since the fitted values are “forecast” over a 1-month window and also the accuracy measures are averaged over a much longer timeframe.

Finally, the dashed vertical lines in the plots highlight that 2011/01 represents a training/test split that produces stronger performance metrics for the test set than the training set, which confirms our observation in part (d) above.