Exercise Week 5: Solutions

library(fpp3)

fpp3 5.10, Ex 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)
  • Bricks (aus_production)
  • NSW Lambs (aus_livestock)
  • Household wealth (hh_budget)
  • Australian takeaway food turnover (aus_retail)

Australian population

global_economy |>
  filter(Country == "Australia") |>
  autoplot(Population)

Data has trend and no seasonality. Random walk with drift model is appropriate.

global_economy |>
  filter(Country == "Australia") |>
  model(RW(Population ~ drift())) |>
  forecast(h = "10 years") |>
  autoplot(global_economy)

Australian clay brick production

aus_production |>
  filter(!is.na(Bricks)) |>
  autoplot(Bricks) +
  labs(title = "Clay brick production")

This data appears to have more seasonality than trend, so of the models available, seasonal naive is most appropriate.

aus_production |>
  filter(!is.na(Bricks)) |>
  model(SNAIVE(Bricks)) |>
  forecast(h = "5 years") |>
  autoplot(aus_production)

NSW Lambs

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

This data appears to have more seasonality than trend, so of the models available, seasonal naive is most appropriate.

nsw_lambs |>
  model(SNAIVE(Count)) |>
  forecast(h = "5 years") |>
  autoplot(nsw_lambs)

Household wealth

hh_budget |>
  autoplot(Wealth)

Annual data with trend upwards, so we can use a random walk with drift.

hh_budget |>
  model(RW(Wealth ~ drift())) |>
  forecast(h = "5 years") |>
  autoplot(hh_budget)

Australian takeaway food turnover

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

This data has strong seasonality and strong trend, so we will use a seasonal naive model with drift.

takeaway |>
  model(SNAIVE(Turnover ~ drift())) |>
  forecast(h = "5 years") |>
  autoplot(takeaway)

This is actually not one of the four benchmark methods discussed in the book, but is sometimes a useful benchmark when there is strong seasonality and strong trend.

The corresponding equation is \hat{y}_{T+h|T} = y_{T+h-m(k+1)} + \frac{h}{T-m}\sum_{t=m+1}^T(y_t - y_{t-m}), where m=12 and k is the integer part of (h-1)/m (i.e., the number of complete years in the forecast period prior to time T+h).

fpp3 5.10, Ex 2

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

  1. Produce a time plot of the series.
fb_stock <- gafa_stock |>
  filter(Symbol == "FB")
fb_stock |>
  autoplot(Close)

An upward trend is evident until mid-2018, after which the closing stock price drops.

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

The data must be made regular before it can be modelled. We will use trading days as our regular index.

fb_stock <- fb_stock |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE)

Time to model a random walk with drift.

fb_stock |>
  model(RW(Close ~ drift())) |>
  forecast(h = 100) |>
  autoplot(fb_stock)

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

Prove drift methods are extrapolations from the first and last observation. First, we will demonstrate it graphically.

fb_stock |>
  model(RW(Close ~ drift())) |>
  forecast(h = 100) |>
  autoplot(fb_stock) +
  geom_line(
    aes(y = Close),
    linetype = "dashed", colour = "blue",
    data = fb_stock |> filter(trading_day %in% range(trading_day))
  )

To prove it algebraically, note that \begin{align*} \hat{y}_{T+h|T} = y_T + h\left(\frac{y_T-y_1}{T-1}\right) \end{align*} which is a straight line with slope (y_T-y_1)/(T-1) that goes through the point (T,y_T).

Therefore, it must also go through the point (1,c) where (y_T-c)/(T-1) = (y_T - y_1) / (T-1), so c=y_1.

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

Use other appropriate benchmark methods. The most appropriate benchmark method is the naive model. The mean forecast is terrible for this type of data, and the data is non-seasonal.

fb_stock |>
  model(NAIVE(Close)) |>
  forecast(h = 100) |>
  autoplot(fb_stock)

The naive method is most appropriate, and will also be best if the efficient market hypothesis holds true.

fpp3 5.10, Ex 3

Apply a seasonal naïve method to the quarterly Australian beer production data from 1992. Check if the residuals look like white noise, and plot the forecasts. The following code will help.

# 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()

  • The residuals are not centred around 0 (typically being slightly below it), this is due to the model failing to capture the negative trend in the data.
  • Peaks and troughs in residuals spaced roughly 4 observations apart are apparent leading to a negative spike at lag 4 in the ACF. So they do not resemble white noise. Lags 1 and 3 are also significant, however they are very close to the threshold and are of little concern.
  • The distribution of the residuals does not appear very normal, however it is probably close enough for the accuracy of our intervals (it being not centred on 0 is more concerning).
# Look at some forecasts
fit |>
  forecast() |>
  autoplot(recent_production)

The forecasts look reasonable, although the intervals may be a bit wide. This is likely due to the slight trend not captured by the model (which subsequently violates the assumptions imposed on the residuals).

fpp3 5.10, Ex 4

Repeat the exercise for 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

The data does not contain seasonality, so the naive model is more appropriate.

# Extract data of interest
aus_exports <- filter(global_economy, Country == "Australia")
# Define and estimate a model
fit <- aus_exports |> model(NAIVE(Exports))
# Check residuals
fit |> gg_tsresiduals()

  • The ACF plot reveals that the first lag exceeds the significance threshold.
  • This data may still be white noise, as it is the only lag that exceeds the blue dashed lines (5% of the lines are expected to exceed it). However as it is the first lag, it is probable that there exists some real auto-correlation in the residuals that can be modelled.
  • The distribution appears normal.
  • The residual plot appears mostly random, however more observations appear to be above zero. This again, is due to the model not capturing the trend.
# Look at some forecasts
fit |>
  forecast() |>
  autoplot(aus_exports)

  • The forecasts appear reasonable as the series appears to have flattened in recent years.
  • The intervals are also reasonable — despite the assumptions behind them having been violated.

Australian brick production

The data is seasonal, so the seasonal naive model is more appropriate.

# Remove the missing values at the end of the series
tidy_bricks <- aus_production |>
  filter(!is.na(Bricks))
# Define and estimate a model
fit <- tidy_bricks |>
  model(SNAIVE(Bricks))
# Look at the residuals
fit |> gg_tsresiduals()

  • The residual plot does not appear random. Periods of low production and high production are evident, leading to autocorrelation in the residuals.

  • The residuals from this model are not white noise. The ACF plot shows a strong sinusoidal pattern of decay, indicating that the residuals are auto-correlated.

  • The histogram is also not normally distributed, as it has a long left tail.

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

  • The point forecasts appear reasonable as the series appears to have flattened in recent years.
  • The intervals appear much larger than necessary.

fpp3 5.10, Ex 5

Produce forecasts for the 7 Victorian series in aus_livestock using SNAIVE(). Plot the resulting forecasts including the historical data. Is this a reasonable benchmark for these series?

aus_livestock |>
  filter(State == "Victoria") |>
  model(SNAIVE(Count)) |>
  forecast(h = "5 years") |>
  autoplot(aus_livestock)

  • Most point forecasts look reasonable from the seasonal naive method.
  • Some series are more seasonal than others, and for the series with very weak seasonality it may be better to consider using a naive or drift method.
  • The prediction intervals in some cases go below zero, so perhaps a log transformation would have been better for these series.

fpp3 5.10, Ex 6

Are the following statements true or false? Explain your answer.

  1. Good forecast methods should have normally distributed residuals.
  • False.
  • Although many good forecasting methods produce normally distributed residuals this is not required to produce good forecasts.
  • Other forecasting methods may use other distributions, it is just less common as they can be more difficult to work with.
  1. A model with small residuals will give good forecasts.
  • False.
  • It is possible to produce a model with small residuals by making a highly complicated (overfitted) model that fits the data extremely well.
  • This highly complicated model will often perform very poorly when forecasting new data.
  1. The best measure of forecast accuracy is MAPE.
  • False.
  • There is no single best measure of accuracy - often you would want to see a collection of accuracy measures as they can reveal different things about your residuals. MAPE in particular has some substantial disadvantages - extreme values can result when y_t is close to zero, and it assumes that the unit being measured has a meaningful zero.
  1. If your model doesn’t forecast well, you should make it more complicated.
  • False.
  • There are many reasons why a model may not forecast well, and making the model more complicated can make the forecasts worse.
  • The model specified should capture structures that are evident in the data. Although adding terms that are unrelated to the structures found in the data will improve the model’s residuals, the forecasting performance of the model will not improve.
  • Adding missing features relevant to the data (such as including a seasonal pattern that exists in the data) should improve forecast performance.
  1. Always choose the model with the best forecast accuracy as measured on the test set.
  • False.
  • There are many measures of forecast accuracy, and the appropriate model is the one which is best suited to the forecasting task. For instance, you may be interested in choosing a model which forecasts well for predictions exactly one year ahead. In this case, using cross-validated accuracy could be a more useful approach to evaluating accuracy.

fpp3 5.10, Ex 8

Consider the number of pigs slaughtered in New South Wales (data set aus_livestock).

  1. Produce some plots of the data in order to become familiar with it.
nsw_pigs <- aus_livestock |>
  filter(State == "New South Wales", Animal == "Pigs")
nsw_pigs |>
  autoplot(Count)

Data generally follows a downward trend, however there are some periods where the amount of pigs slaughtered changes rapidly.

nsw_pigs |> gg_season(Count, labels = "right")

nsw_pigs |> gg_subseries(Count)

Some seasonality is apparent, with notable increases in December and decreases during January, February and April.

  1. Create a training set of 486 observations, withholding a test set of 72 observations (6 years).
nsw_pigs_train <- nsw_pigs |> slice(1:486)
  1. Try using various benchmark methods to forecast the training set and compare the results on the test set. Which method did best?
fit <- nsw_pigs_train |>
  model(
    mean = MEAN(Count),
    naive = NAIVE(Count),
    snaive = SNAIVE(Count),
    drift = RW(Count ~ drift())
  )
fit |>
  forecast(h = "6 years") |>
  accuracy(nsw_pigs)
# A tibble: 4 × 12
  .model Animal State       .type      ME   RMSE    MAE    MPE  MAPE  MASE RMSSE
  <chr>  <fct>  <fct>       <chr>   <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
1 drift  Pigs   New South … Test   -4685.  8091.  6967.  -7.36  10.1 0.657 0.557
2 mean   Pigs   New South … Test  -39360. 39894. 39360. -55.9   55.9 3.71  2.75 
3 naive  Pigs   New South … Test   -6138.  8941.  7840.  -9.39  11.4 0.740 0.615
4 snaive Pigs   New South … Test   -5838. 10111.  8174.  -8.81  11.9 0.771 0.696
# ℹ 1 more variable: ACF1 <dbl>

The drift method performed best for all measures of accuracy (although it had a larger first order auto-correlation)

  1. Check the residuals of your preferred method. Do they resemble white noise?
fit |>
  select(drift) |>
  gg_tsresiduals()

The residuals do not appear to be white noise as the ACF plot contains many significant lags. It is also clear that the seasonal component is not captured by the drift method, as there exists a strong positive auto-correlation at lag 12 (1 year). The histogram appears to have a slightly long left tail.

fpp3 5.10, Ex 11

We will use the bricks data from aus_production (Australian quarterly clay brick production 1956–2005) for this exercise.

  1. Use an STL decomposition to calculate the trend-cycle and seasonal indices. (Experiment with having fixed or changing seasonality.)
tidy_bricks <- aus_production |>
  filter(!is.na(Bricks))
tidy_bricks |>
  model(STL(Bricks)) |>
  components() |>
  autoplot()

Data is multiplicative, and so a transformation should be used.

dcmp <- tidy_bricks |>
  model(STL(log(Bricks))) |>
  components()
dcmp |>
  autoplot()

Seasonality varies slightly.

dcmp <- tidy_bricks |>
  model(stl = STL(log(Bricks) ~ season(window = "periodic"))) |>
  components()
dcmp |> autoplot()

The seasonality looks fairly stable, so I’ve used a periodic season (window). The decomposition still performs well when the seasonal component is fixed. The remainder term does not appear to contain a substantial amount of seasonality.

  1. Compute and plot the seasonally adjusted data.
dcmp |>
  as_tsibble() |>
  autoplot(season_adjust)

  1. Use a naïve method to produce forecasts of the seasonally adjusted data.
fit <- dcmp |>
  select(-.model) |>
  model(naive = NAIVE(season_adjust)) |>
  forecast(h = "5 years")
dcmp |>
  as_tsibble() |>
  autoplot(season_adjust) + autolayer(fit)

  1. Use decomposition_model() to reseasonalise the results, giving forecasts for the original data.
fit <- tidy_bricks |>
  model(stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)))
fit |>
  forecast(h = "5 years") |>
  autoplot(tidy_bricks)

  1. Do the residuals look uncorrelated?
fit |> gg_tsresiduals()

The residuals do not appear uncorrelated as there are several lags of the ACF which exceed the significance threshold.

  1. Repeat with a robust STL decomposition. Does it make much difference?
fit_robust <- tidy_bricks |>
  model(stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)))

fit_robust |> gg_tsresiduals()

The residuals appear slightly less auto-correlated, however there is still significant auto-correlation at lag 8.

  1. Compare forecasts from decomposition_model() with those from SNAIVE(), using a test set comprising the last 2 years of data. Which is better?
tidy_bricks_train <- tidy_bricks |>
  slice(1:(n() - 8))
fit <- tidy_bricks_train |>
  model(
    stl_mdl = decomposition_model(STL(log(Bricks)), NAIVE(season_adjust)),
    snaive = SNAIVE(Bricks)
  )

fc <- fit |>
  forecast(h = "2 years")
fc |>
  autoplot(tidy_bricks, level = NULL)

The decomposition forecasts appear to more closely follow the actual future data.

fc |>
  accuracy(tidy_bricks)
# A tibble: 2 × 10
  .model  .type    ME  RMSE   MAE     MPE  MAPE  MASE RMSSE    ACF1
  <chr>   <chr> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>   <dbl>
1 snaive  Test  2.75   20    18.2  0.395   4.52 0.504 0.407 -0.0503
2 stl_mdl Test  0.368  18.1  15.1 -0.0679  3.76 0.418 0.368  0.115

The STL decomposition forecasts are more accurate than the seasonal naive forecasts across all accuracy measures.

fpp3 5.10, Ex 12

tourism contains quarterly visitor nights (in thousands) from 1998 to 2017 for 76 regions of Australia.

  1. Extract data from the Gold Coast region using filter() and aggregate total overnight trips (sum over Purpose) using summarise(). Call this new dataset gc_tourism.
gc_tourism <- tourism |>
  filter(Region == "Gold Coast") |>
  summarise(Trips = sum(Trips))
gc_tourism
# A tsibble: 80 x 2 [1Q]
   Quarter Trips
     <qtr> <dbl>
 1 1998 Q1  827.
 2 1998 Q2  681.
 3 1998 Q3  839.
 4 1998 Q4  820.
 5 1999 Q1  987.
 6 1999 Q2  751.
 7 1999 Q3  822.
 8 1999 Q4  914.
 9 2000 Q1  871.
10 2000 Q2  780.
# ℹ 70 more rows
  1. Using slice() or filter(), create three training sets for this data excluding the last 1, 2 and 3 years. For example, gc_train_1 <- gc_tourism |> slice(1:(n()-4)).
gc_train_1 <- gc_tourism |> slice(1:(n() - 4))
gc_train_2 <- gc_tourism |> slice(1:(n() - 8))
gc_train_3 <- gc_tourism |> slice(1:(n() - 12))
  1. Compute one year of forecasts for each training set using the seasonal naïve (SNAIVE()) method. Call these gc_fc_1, gc_fc_2 and gc_fc_3, respectively.
gc_fc <- bind_cols(
  gc_train_1 |> model(gc_fc_1 = SNAIVE(Trips)),
  gc_train_2 |> model(gc_fc_2 = SNAIVE(Trips)),
  gc_train_3 |> model(gc_fc_3 = SNAIVE(Trips))
) |> forecast(h = "1 year")
gc_fc |> autoplot(gc_tourism)

  1. Use accuracy() to compare the test set forecast accuracy using MAPE. Comment on these.
gc_fc |> accuracy(gc_tourism)
# A tibble: 3 × 10
  .model  .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
  <chr>   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
1 gc_fc_1 Test   75.1 167.  154.   6.36 15.1  2.66  2.36  -0.410
2 gc_fc_2 Test   12.0  43.1  39.5  1.14  4.32 0.670 0.599 -0.792
3 gc_fc_3 Test   35.8  91.4  83.9  3.56  9.07 1.46  1.30   0.239

The second set of forecasts are most accurate (as can be seen in the previous plot), however this is likely due to chance.