knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(fpp3)

Exercise 5.1

Produce forecasts for the following series using whichever of NAIVE(y), SNAIVE(y) or RW(y ~ drift()) is more appropriate in each case.

Australian Population (global_economy)

Annual data with a steady upward trend and no seasonality, so the drift method is the right choice.

aus_pop <- global_economy |> filter(Country == "Australia")

aus_pop |>
  model(Drift = RW(Population ~ drift())) |>
  forecast(h = 10) |>
  autoplot(aus_pop) +
  labs(title = "Australian population, drift forecasts", y = "People")

Bricks (aus_production)

Quarterly data with strong seasonality, so the seasonal naive method is appropriate. The series has missing values after 2005 Q3, so I drop them first.

bricks <- aus_production |> filter(!is.na(Bricks))

bricks |>
  model(SNaive = SNAIVE(Bricks)) |>
  forecast(h = "5 years") |>
  autoplot(bricks) +
  labs(title = "Australian clay brick production, seasonal naive forecasts",
       y = "Millions of bricks")

NSW Lambs (aus_livestock)

Monthly data with seasonality, so the seasonal naive method is appropriate.

nsw_lambs <- aus_livestock |>
  filter(State == "New South Wales", Animal == "Lambs")

nsw_lambs |>
  model(SNaive = SNAIVE(Count)) |>
  forecast(h = "3 years") |>
  autoplot(nsw_lambs) +
  labs(title = "NSW lambs slaughtered, seasonal naive forecasts", y = "Count")

Household wealth (hh_budget)

Annual data with an upward trend and no seasonality, so the drift method is appropriate.

hh_budget |>
  model(Drift = RW(Wealth ~ drift())) |>
  forecast(h = 5) |>
  autoplot(hh_budget) +
  labs(title = "Household wealth, drift forecasts", y = "% of net disposable income")

Australian takeaway food turnover (aus_retail)

Monthly data with strong seasonality, so the seasonal naive method is appropriate.

takeaway <- aus_retail |>
  filter(Industry == "Takeaway food services") |>
  summarise(Turnover = sum(Turnover))

takeaway |>
  model(SNaive = SNAIVE(Turnover)) |>
  forecast(h = "3 years") |>
  autoplot(takeaway) +
  labs(title = "Australian takeaway food turnover, seasonal naive forecasts",
       y = "$ million")

The seasonal naive forecasts repeat the last year of data. Since this series is still trending upward, the forecasts are probably a little low. Adding a drift term would fix that, but on its own SNAIVE() is still the standard seasonal benchmark.

Exercise 5.2

Use the Facebook stock price (data set gafa_stock) to do the following.

a. Produce a time plot of the series.

Stock prices are only recorded on trading days, so the dates are irregular. I re-index the data by trading day so that the forecasting functions work.

fb_stock <- gafa_stock |>
  filter(Symbol == "FB") |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)

fb_stock |>
  autoplot(Close) +
  labs(title = "Facebook closing stock price", x = "Trading day", y = "US$")

b. Produce forecasts using the drift method and plot them.

fb_fit <- fb_stock |> model(Drift = RW(Close ~ drift()))

fb_fc <- fb_fit |> forecast(h = 63)

fb_fc |>
  autoplot(fb_stock) +
  labs(title = "Facebook closing price, drift forecasts (about 3 months ahead)",
       x = "Trading day", y = "US$")

c. Show that the forecasts are identical to extending the line drawn between the first and last observations.

first_obs <- fb_stock |> filter(day == min(day))
last_obs <- fb_stock |> filter(day == max(day))

fb_fc |>
  autoplot(fb_stock, level = NULL) +
  geom_segment(aes(x = first_obs$day, y = first_obs$Close,
                   xend = last_obs$day, yend = last_obs$Close),
               colour = "red", linetype = "dashed") +
  labs(title = "Drift forecasts and the line through the first and last observations",
       x = "Trading day", y = "US$")

The drift forecast is a straight line whose slope is (last value - first value) / (number of periods - 1):

n <- nrow(fb_stock)
slope <- (last_obs$Close - first_obs$Close) / (n - 1)
slope
## [1] 0.06076372
# the drift forecast and the extended line give the same values
tibble(h = c(1, 63),
       drift_forecast = fb_fc$.mean[c(1, 63)],
       extended_line = last_obs$Close + slope * c(1, 63))
## # A tibble: 2 x 3
##       h drift_forecast extended_line
##   <dbl>          <dbl>         <dbl>
## 1     1           131.          131.
## 2    63           135.          135.

The dashed red line passes exactly through the forecasts, and the two columns above match, so the drift forecasts are the extension of the line joining the first and last observations.

d. Try using some of the other benchmark functions to forecast the same data set. Which do you think is best? Why?

fb_stock |>
  model(Mean = MEAN(Close), Naive = NAIVE(Close), Drift = RW(Close ~ drift())) |>
  forecast(h = 63) |>
  autoplot(fb_stock, level = NULL) +
  labs(title = "Benchmark forecasts for Facebook closing price",
       x = "Trading day", y = "US$")

The mean method is clearly the worst. It forecasts about $120, which ignores everything the price has done recently.

The naive method is the best choice here. Stock prices behave like a random walk, so the most recent price is the best guess for tomorrow’s price. The drift method extends the overall upward slope of the whole series, but the price had been falling for the last few months of the data, so continuing to forecast growth is not justified. Drift also gives wider forecast intervals.

Exercise 5.3

Apply a seasonal naive method to the quarterly Australian beer production data from 1992. Check if the residuals look like white noise, and plot the forecasts.

# Extract data of interest
recent_production <- aus_production |>
  filter(year(Quarter) >= 1992)

# Define and estimate a model
fit <- recent_production |> model(SNAIVE(Beer))

# Look at the residuals
fit |> gg_tsresiduals()

# Look at some forecasts
fit |> forecast() |> autoplot(recent_production)

fit |> augment() |> features(.innov, ljung_box, lag = 8)
## # A tibble: 1 x 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 SNAIVE(Beer)    32.3 0.0000834

What do you conclude?

The residuals are close to white noise, but not quite. The ACF plot has one significant spike at lag 4, which says that some of the year-to-year pattern is still not captured. The Ljung-Box test gives a p-value of about 0.0001, so we reject the idea that the residuals are completely uncorrelated. On the other hand, the residuals have a mean near zero, roughly constant variance and a histogram that is close to normal, so the prediction intervals should be reasonable. Overall the seasonal naive method is a decent benchmark for this series, but a better model could still extract the remaining information at lag 4.

Exercise 5.4

Repeat the previous exercise using the Australian Exports series from global_economy and the Bricks series from aus_production. Use whichever of NAIVE() or SNAIVE() is more appropriate in each case.

Australian Exports

This is annual data, so there is no seasonality and NAIVE() is the appropriate method.

aus_exports <- global_economy |> filter(Country == "Australia")

exports_fit <- aus_exports |> model(NAIVE(Exports))

exports_fit |> gg_tsresiduals()

exports_fit |> forecast(h = 10) |> autoplot(aus_exports) +
  labs(title = "Australian exports, naive forecasts", y = "% of GDP")

exports_fit |> augment() |> features(.innov, ljung_box, lag = 10)
## # A tibble: 1 x 4
##   Country   .model         lb_stat lb_pvalue
##   <fct>     <chr>            <dbl>     <dbl>
## 1 Australia NAIVE(Exports)    16.4    0.0896

The residuals look close to white noise. Almost all ACF spikes are inside the bounds, and the Ljung-Box test gives a p-value of about 0.09, so there is not enough evidence of autocorrelation. The residual mean is near zero and the histogram is roughly normal, though slightly skewed. The naive method is a reasonable benchmark here.

Bricks

This is quarterly data with clear seasonality, so SNAIVE() is appropriate.

bricks_fit <- bricks |> model(SNAIVE(Bricks))

bricks_fit |> gg_tsresiduals()

bricks_fit |> forecast(h = 12) |> autoplot(bricks) +
  labs(title = "Brick production, seasonal naive forecasts", y = "Millions of bricks")

bricks_fit |> augment() |> features(.innov, ljung_box, lag = 8)
## # A tibble: 1 x 3
##   .model         lb_stat lb_pvalue
##   <chr>            <dbl>     <dbl>
## 1 SNAIVE(Bricks)    274.         0

The residuals are clearly not white noise. The ACF has large significant spikes at many lags, especially the early ones, which shows that the trend and the cycles in brick production are not captured by the seasonal naive method. The Ljung-Box p-value is essentially 0. The residuals also do not have constant variance: they are much more variable in the middle of the series. There is still a lot of information left in the residuals, so this benchmark could be beaten easily.

Exercise 5.7

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

a. Create a training dataset consisting of observations before 2011.

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
distinct(as_tibble(myseries), State, Industry)
## # A tibble: 1 x 2
##   State           Industry              
##   <chr>           <chr>                 
## 1 New South Wales Takeaway food services
myseries_train <- myseries |>
  filter(year(Month) < 2011)

b. Check that your data have been split appropriately.

autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red") +
  labs(title = "Retail turnover: training data in red, full series in black",
       y = "$ million")

c. Fit a seasonal naive model to the training data.

fit <- myseries_train |>
  model(SNAIVE(Turnover))

d. Check the residuals. Do the residuals appear to be uncorrelated and normally distributed?

fit |> gg_tsresiduals()

fit |> augment() |> features(.innov, ljung_box, lag = 24)
## # A tibble: 1 x 5
##   State           Industry               .model           lb_stat lb_pvalue
##   <chr>           <chr>                  <chr>              <dbl>     <dbl>
## 1 New South Wales Takeaway food services SNAIVE(Turnover)   1101.         0

No. The residuals are strongly correlated. The ACF has large significant spikes at almost every lag, and they decay only slowly, which is what happens when a trend is left in the residuals. The Ljung-Box p-value is essentially 0. The residuals are also not centred on zero: most are positive, because the series keeps growing and the seasonal naive method always forecasts last year’s value. The histogram is right-skewed rather than normal, and the variance grows over time. The prediction intervals from this model will not be reliable.

e. Produce forecasts for the test data.

fc <- fit |>
  forecast(new_data = anti_join(myseries, myseries_train))

fc |> autoplot(myseries) +
  labs(title = "Seasonal naive forecasts against the actual data", y = "$ million")

f. Compare the accuracy of your forecasts against the actual values.

fit |> accuracy()
## # A tibble: 1 x 12
##   State    Industry .model .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>    <chr>    <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 New Sou~ Takeawa~ SNAIV~ Trai~  11.5  26.1  19.2  4.81  9.59     1     1 0.890
fc |> accuracy(myseries)
## # A tibble: 1 x 12
##   .model    State Industry .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>     <chr> <chr>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(T~ New ~ Takeawa~ Test   48.6  96.8  79.5  7.67  16.3  4.14  3.71 0.964

The test set errors are much larger than the training errors. The RMSE rises from about 26 to about 97, and the MAPE from about 10% to about 16%. The MASE of about 4 means the forecasts are roughly four times worse than a one-step seasonal naive forecast on the training data. The forecast plot shows why: turnover kept growing after 2010, but the seasonal naive forecasts just repeat the 2010 values, so they fall further behind as the horizon grows.

g. How sensitive are the accuracy measures to the amount of training data used?

train_test <- function(cutoff) {
  train <- myseries |> filter(year(Month) < cutoff)
  train |>
    model(SNAIVE(Turnover)) |>
    forecast(new_data = anti_join(myseries, train)) |>
    accuracy(myseries) |>
    mutate(training_ends = cutoff - 1) |>
    select(training_ends, RMSE, MAE, MAPE, MASE)
}

bind_rows(train_test(2005), train_test(2011), train_test(2015))
## # A tibble: 3 x 5
##   training_ends  RMSE   MAE  MAPE  MASE
##           <dbl> <dbl> <dbl> <dbl> <dbl>
## 1          2004 166.  133.   29.7  7.93
## 2          2010  96.8  79.5  16.3  4.14
## 3          2014 118.  112.   20.8  5.30

The accuracy measures are quite sensitive, but the length of the training set is not the only thing that matters. Using data up to 2004 gives by far the worst results (MAPE of about 30%), because the test period is 14 years long and the series grows a great deal over it. Training up to 2010 is the best of the three (MAPE of about 16%). Training up to 2014 leaves only a four-year test set, yet it is worse again (MAPE of about 21%), because turnover grew very quickly in 2016 to 2018 and the seasonal naive forecasts cannot follow it.

So the errors depend on both the forecast horizon and on what happens during the test period. Longer horizons generally give bigger errors, and any method that ignores the trend will look worse when the test period contains rapid growth.