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

Exercise 8.1

Consider the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.

pigs <- aus_livestock |>
  filter(State == "Victoria", Animal == "Pigs")

pigs |> autoplot(Count) +
  labs(title = "Pigs slaughtered in Victoria", y = "Count")

a. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of alpha and l0, and generate forecasts for the next four months.

fit <- pigs |>
  model(ses = ETS(Count ~ error("A") + trend("N") + season("N")))

report(fit)
## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.3221247 
## 
##   Initial states:
##      l[0]
##  100646.6
## 
##   sigma^2:  87480760
## 
##      AIC     AICc      BIC 
## 13737.10 13737.14 13750.07

The optimal values are alpha = 0.3221 and l0 = 100,646.6. The small alpha means the level is updated slowly, so the forecasts rely on a long stretch of history rather than just the most recent months.

fc <- fit |> forecast(h = 4)
fc
## # A fable: 4 x 6 [1M]
## # Key:     Animal, State, .model [1]
##   Animal State    .model    Month
##   <fct>  <fct>    <chr>     <mth>
## 1 Pigs   Victoria ses    2019 Jan
## 2 Pigs   Victoria ses    2019 Feb
## 3 Pigs   Victoria ses    2019 Mar
## 4 Pigs   Victoria ses    2019 Apr
## # i 2 more variables: Count <dist>, .mean <dbl>
fc |> autoplot(pigs |> filter(year(Month) >= 2010)) +
  labs(title = "Four-month forecasts from simple exponential smoothing", y = "Count")

All four forecasts equal 95,187, which is what simple exponential smoothing always does: the forecast is flat for every horizon.

b. Compute a 95% prediction interval for the first forecast using y +/- 1.96s, where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

s <- sd(augment(fit)$.resid)
s
## [1] 9344.666
y_hat <- fc$.mean[1]
c(lower = y_hat - 1.96 * s, upper = y_hat + 1.96 * s)
##     lower     upper 
##  76871.01 113502.10
fc |> hilo(95) |> pull(`95%`) |> head(1)
## <hilo[1]>
## [1] [76854.79, 113518.3]95

The hand calculation gives [76,871, 113,502] and R gives [76,855, 113,518]. The two are almost identical, differing by only about 16 units at each end, which is less than 0.1% of the forecast. The small gap comes from the fact that R uses the model’s estimated variance, which divides by the number of observations, while sd() divides by n - 1 and also subtracts the residual mean. Either way the interval is very wide, roughly plus or minus 19% of the forecast, which reflects how noisy this series is.

Exercise 8.5

Data set global_economy contains the annual Exports from many countries. Select one country to analyse. I use Australia.

a. Plot the Exports series and discuss the main features of the data.

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

aus_exports |> autoplot(Exports) +
  labs(title = "Australian exports", y = "% of GDP")

Exports rose from about 13% of GDP in the 1960s to about 21% in 2017, so there is a clear long-run upward trend, but it is not steady. The series was flat through the 1960s, jumped in the 1970s, fell back in the late 1980s, and rose again from the mid-1990s. It peaked at 23% in 2009 and has drifted sideways since. The data is annual, so there is no seasonality, and the year-to-year movements look like irregular cycles tied to commodity prices and the exchange rate rather than a fixed pattern.

b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.

fit_ann <- aus_exports |>
  model(ann = ETS(Exports ~ error("A") + trend("N") + season("N")))

report(fit_ann)
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.5659948 
## 
##   Initial states:
##      l[0]
##  12.98943
## 
##   sigma^2:  1.3621
## 
##      AIC     AICc      BIC 
## 257.3943 257.8387 263.5756
fit_ann |> forecast(h = 10) |> autoplot(aus_exports) +
  labs(title = "ETS(A,N,N) forecasts of Australian exports", y = "% of GDP")

c. Compute the RMSE values for the training data.

accuracy(fit_ann) |> select(.model, RMSE, MAE, MASE)
## # A tibble: 1 x 4
##   .model  RMSE   MAE  MASE
##   <chr>  <dbl> <dbl> <dbl>
## 1 ann     1.15 0.914 0.928

The training RMSE is 1.147.

d. Compare the results to those from an ETS(A,A,N) model.

fit_both <- aus_exports |>
  model(ann = ETS(Exports ~ error("A") + trend("N") + season("N")),
        aan = ETS(Exports ~ error("A") + trend("A") + season("N")))

accuracy(fit_both) |> select(.model, RMSE, MAE, MASE)
## # A tibble: 2 x 4
##   .model  RMSE   MAE  MASE
##   <chr>  <dbl> <dbl> <dbl>
## 1 ann     1.15 0.914 0.928
## 2 aan     1.12 0.893 0.907
glance(fit_both) |> select(.model, sigma2, AIC, AICc, BIC)
## # A tibble: 2 x 5
##   .model sigma2   AIC  AICc   BIC
##   <chr>   <dbl> <dbl> <dbl> <dbl>
## 1 ann      1.36  257.  258.  264.
## 2 aan      1.34  258.  259.  269.
tidy(fit_both)
## # A tibble: 6 x 4
##   Country   .model term   estimate
##   <fct>     <chr>  <chr>     <dbl>
## 1 Australia ann    alpha  0.566   
## 2 Australia ann    l[0]  13.0     
## 3 Australia aan    alpha  0.441   
## 4 Australia aan    beta   0.000100
## 5 Australia aan    l[0]  12.7     
## 6 Australia aan    b[0]   0.137

The trended model fits very slightly better on the training data (RMSE 1.117 against 1.147), but it uses two extra parameters (beta and b0) to get there, and the estimated beta is tiny (0.0001), so the slope barely adapts and is essentially fixed at its initial value of 0.14 percentage points per year. The information criteria, which penalise the extra parameters, prefer the simpler model: AICc is 257.8 for ETS(A,N,N) against 259.5 for ETS(A,A,N). A 3% improvement in RMSE is not worth two extra parameters here, especially since the upward movement in this series comes in irregular bursts rather than as a constant trend.

e. Compare the forecasts from both methods. Which do you think is best?

fit_both |> forecast(h = 10) |> autoplot(aus_exports, level = NULL) +
  labs(title = "ETS(A,N,N) and ETS(A,A,N) forecasts", y = "% of GDP")

ETS(A,N,N) forecasts a flat 20.6%, while ETS(A,A,N) adds a steady upward slope of about 0.14 percentage points per year, reaching 22.1% after ten years. The gap widens the further ahead you look, so the choice does matter at longer horizons. I prefer ETS(A,N,N). It has the better AICc, it is simpler, and extrapolating a trend ten years into the future is risky for a series driven by commodity cycles, which have gone up and down repeatedly over this period.

f. Calculate a 95% prediction interval for the first forecast for each model, using the RMSE values and assuming normal errors. Compare your intervals with those produced using R.

rmse <- accuracy(fit_both) |> select(.model, RMSE)
first_fc <- fit_both |> forecast(h = 1) |> as_tibble() |> select(.model, .mean)

left_join(first_fc, rmse, by = ".model") |>
  mutate(lower = .mean - 1.96 * RMSE,
         upper = .mean + 1.96 * RMSE)
## # A tibble: 2 x 5
##   .model .mean  RMSE lower upper
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 ann     20.6  1.15  18.4  22.9
## 2 aan     20.8  1.12  18.6  23.0
fit_both |> forecast(h = 1) |> hilo(95) |> select(.model, Year, .mean, `95%`)
## # A tsibble: 2 x 4 [1Y]
## # Key:       .model [2]
##   .model  Year .mean                  `95%`
##   <chr>  <dbl> <dbl>                 <hilo>
## 1 ann     2018  20.6 [18.31970, 22.89462]95
## 2 aan     2018  20.8 [18.57028, 23.10700]95

The hand-calculated intervals are [18.36, 22.86] for ETS(A,N,N) and [18.65, 23.03] for ETS(A,A,N). R gives [18.32, 22.89] and [18.57, 23.11]. The intervals match closely, with R’s being very slightly wider in both cases, because the RMSE divides the sum of squared errors by the number of observations while R’s variance estimate accounts for the parameters that were estimated.

Exercise 8.6

Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation.

china <- global_economy |> filter(Country == "China")

china |> autoplot(GDP) +
  labs(title = "Chinese GDP", y = "US$")

lambda <- china |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.03446284
fit_china <- china |>
  model(ets = ETS(GDP),
        damped = ETS(GDP ~ trend("Ad")),
        log = ETS(log(GDP)),
        box_cox = ETS(box_cox(GDP, lambda)))

fit_china
## # A mable: 1 x 5
## # Key:     Country [1]
##   Country          ets        damped          log      box_cox
##   <fct>        <model>       <model>      <model>      <model>
## 1 China   <ETS(M,A,N)> <ETS(M,Ad,N)> <ETS(A,A,N)> <ETS(A,A,N)>
glance(fit_china) |> select(.model, AICc)
## # A tibble: 4 x 2
##   .model    AICc
##   <chr>    <dbl>
## 1 ets     3103. 
## 2 damped  3107. 
## 3 log      -31.9
## 4 box_cox -137.
fit_china |>
  forecast(h = 20) |>
  autoplot(china, level = NULL) +
  labs(title = "Twenty-year forecasts of Chinese GDP under four ETS options", y = "US$")

fit_china |>
  forecast(h = 20) |>
  as_tibble() |>
  filter(Year %in% c(2027, 2037)) |>
  mutate(forecast_trillions = round(.mean / 1e12, 1)) |>
  select(.model, Year, forecast_trillions) |>
  arrange(.model, Year)
## # A tibble: 8 x 3
##   .model   Year forecast_trillions
##   <chr>   <dbl>              <dbl>
## 1 box_cox  2027               47.6
## 2 box_cox  2037              224. 
## 3 damped   2027               18.7
## 4 damped   2037               23.9
## 5 ets      2027               19.7
## 6 ets      2037               27.3
## 7 log      2027               42.8
## 8 log      2037              172.

GDP was about $12.2 trillion in 2017, the last observation. The four options give wildly different answers for 2037:

The intuition: a transformation decides whether growth is treated as constant in dollars (no transformation) or constant in percent (log or Box-Cox), and damping decides whether growth is allowed to continue forever or is gradually choked off. The AICc values cannot be compared across transformations, since the models are fitted to different data. Judging by plausibility, the damped model is the most believable of the four. China’s growth rate has already slowed, and the untransformed exponential forecasts imply an economy many times the size of the world economy today.

Exercise 8.7

Find an ETS model for the Gas data from aus_production and forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?

gas <- aus_production |> select(Gas)

gas |> autoplot(Gas) +
  labs(title = "Australian gas production", y = "Petajoules")

fit_gas <- gas |>
  model(auto = ETS(Gas),
        damped = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
        additive = ETS(Gas ~ error("A") + trend("A") + season("A")))

fit_gas
## # A mable: 1 x 3
##           auto        damped     additive
##        <model>       <model>      <model>
## 1 <ETS(M,A,M)> <ETS(M,Ad,M)> <ETS(A,A,A)>
report(fit_gas |> select(auto))
## Series: Gas 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.6528545 
##     beta  = 0.1441675 
##     gamma = 0.09784922 
## 
##   Initial states:
##      l[0]       b[0]      s[0]    s[-1]    s[-2]     s[-3]
##  5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
## 
##   sigma^2:  0.0032
## 
##      AIC     AICc      BIC 
## 1680.929 1681.794 1711.389

Why multiplicative seasonality is necessary: the seasonal swings grow as the level of the series grows. In the 1950s production was around 5 petajoules and the swings were barely visible, while by the 2000s production was around 220 and the peaks and troughs are about 30 petajoules apart. Additive seasonality assumes the seasonal effect is the same size in every year, which would make the early years far too variable and the later years far too smooth. The automatic selection confirms this by choosing ETS(M,A,M).

accuracy(fit_gas) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 3 x 5
##   .model    RMSE   MAE  MAPE  MASE
##   <chr>    <dbl> <dbl> <dbl> <dbl>
## 1 auto      4.60  3.02  4.08 0.542
## 2 damped    4.59  3.03  4.10 0.544
## 3 additive  4.76  3.35 10.9  0.600
glance(fit_gas) |> select(.model, sigma2, AICc)
## # A tibble: 3 x 3
##   .model     sigma2  AICc
##   <chr>       <dbl> <dbl>
## 1 auto      0.00324 1682.
## 2 damped    0.00329 1685.
## 3 additive 23.6     1873.

The additive model is clearly worse on every measure (MAPE of 10.9% against 4.1%), which confirms the argument above.

fit_gas |>
  select(auto, damped) |>
  forecast(h = "3 years") |>
  autoplot(gas |> filter(year(Quarter) >= 1990), level = NULL) +
  labs(title = "ETS(M,A,M) and damped forecasts of gas production", y = "Petajoules")

Does damping improve the forecasts? Barely, and not enough to prefer it. The damped model has almost the same training RMSE (4.59 against 4.60) but a slightly worse AICc (1685 against 1682). A fairer test is time series cross-validation:

gas_cv <- gas |>
  stretch_tsibble(.init = 60, .step = 1) |>
  model(auto = ETS(Gas ~ error("M") + trend("A") + season("M")),
        damped = ETS(Gas ~ error("M") + trend("Ad") + season("M"))) |>
  forecast(h = 4) |>
  accuracy(gas)

gas_cv |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
##   .model  RMSE   MAE  MAPE  MASE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 auto    7.60  5.69  4.75  1.02
## 2 damped  7.54  5.67  4.76  1.02

Over four-step-ahead forecasts the damped model is marginally better (RMSE 7.54 against 7.60), a difference of less than 1%. Damping is a sensible safeguard for long horizons, since it stops the trend growing without limit, but for this series it makes almost no practical difference.

Exercise 8.8

Recall your retail time series data (from Exercise 7 in Section 2.10).

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 |> autoplot(Turnover) +
  labs(title = "Retail turnover", y = "$ million")

a. Why is multiplicative seasonality necessary for this series?

The size of the seasonal swings grows with the level of the series. In the early 1980s turnover was around $80 million and the December peak stood only a few million above the rest of the year. By 2018 turnover was around $600 million and the December peak is tens of millions above the other months. The seasonal effect is roughly a constant percentage of the level, not a constant number of dollars, which is exactly what multiplicative seasonality models.

b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

fit_hw <- myseries |>
  model(hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

tidy(fit_hw) |> filter(term %in% c("alpha", "beta", "gamma", "phi"))
## # A tibble: 7 x 5
##   State           Industry               .model    term  estimate
##   <chr>           <chr>                  <chr>     <chr>    <dbl>
## 1 New South Wales Takeaway food services hw        alpha  0.715  
## 2 New South Wales Takeaway food services hw        beta   0.0309 
## 3 New South Wales Takeaway food services hw        gamma  0.106  
## 4 New South Wales Takeaway food services hw_damped alpha  0.775  
## 5 New South Wales Takeaway food services hw_damped beta   0.00752
## 6 New South Wales Takeaway food services hw_damped gamma  0.0968 
## 7 New South Wales Takeaway food services hw_damped phi    0.980
fit_hw |> forecast(h = "3 years") |> autoplot(myseries, level = NULL) +
  labs(title = "Holt-Winters' multiplicative forecasts, with and without damping",
       y = "$ million")

The damping parameter is estimated at phi = 0.98, which is close to 1, so the damped forecasts bend only slightly below the undamped ones over three years.

c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

accuracy(fit_hw) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
##   .model     RMSE   MAE  MAPE  MASE
##   <chr>     <dbl> <dbl> <dbl> <dbl>
## 1 hw         9.74  7.13  3.08 0.308
## 2 hw_damped  9.73  7.10  3.03 0.307
glance(fit_hw) |> select(.model, sigma2, AICc)
## # A tibble: 2 x 3
##   .model     sigma2  AICc
##   <chr>       <dbl> <dbl>
## 1 hw        0.00174 4636.
## 2 hw_damped 0.00171 4627.

The one-step RMSE values are effectively tied: 9.74 for Holt-Winters and 9.73 for the damped version. The AICc separates them more clearly, favouring the damped model (4627 against 4636), so I prefer the damped version. It fits just as well, it is more cautious at longer horizons, and the information criterion supports it.

d. Check that the residuals from the best method look like white noise.

best <- myseries |>
  model(ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

best |> gg_tsresiduals()

best |> 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 "ETS(Turnover ~ erro~    25.3     0.392

The residuals do look like white noise. Almost every ACF spike sits inside the significance bounds, and the Ljung-Box test gives a p-value of 0.39, so there is no evidence of leftover autocorrelation. The residuals are centred on zero with roughly constant variance, and the histogram is close to normal with slightly heavy tails. This is a big improvement on the seasonal naive model from Homework 3, whose residuals failed the same test badly.

e. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naive approach from Exercise 7 in Section 5.11?

train <- myseries |> filter(year(Month) < 2011)

fits <- train |>
  model(snaive = SNAIVE(Turnover),
        hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
        hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

fc <- fits |> forecast(h = "8 years")

fc |> autoplot(myseries, level = NULL) +
  labs(title = "Forecasts from the 2010 training set against the actual data",
       y = "$ million")

accuracy(fc, myseries) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 3 x 5
##   .model     RMSE   MAE  MAPE  MASE
##   <chr>     <dbl> <dbl> <dbl> <dbl>
## 1 hw         70.4  60.9  15.3  3.17
## 2 hw_damped  69.0  62.8  14.7  3.27
## 3 snaive     96.8  79.5  16.3  4.14

Yes, comfortably. The seasonal naive method has a test RMSE of 96.8, while Holt-Winters gets 70.4 and the damped version 69.0. That is an improvement of about 29%. The reason is visible in the plot: the seasonal naive forecasts simply repeat 2010 and fall further behind every year, while the Holt-Winters models keep following the upward trend.

Exercise 8.9

For the same retail data, try an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data. How does that compare with your best previous forecasts on the test set?

lambda <- train |> features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.1011142
fits_stl <- train |>
  model(
    stl_ets = decomposition_model(
      STL(box_cox(Turnover, lambda)),
      ETS(season_adjust ~ season("N"))),
    stl_ets_damped = decomposition_model(
      STL(box_cox(Turnover, lambda)),
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N")))
  )

fc_stl <- fits_stl |> forecast(h = "8 years")

fc_stl |> autoplot(myseries, level = NULL) +
  labs(title = "STL plus ETS forecasts on the Box-Cox transformed series",
       y = "$ million")

bind_rows(accuracy(fc, myseries), accuracy(fc_stl, myseries)) |>
  select(.model, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)
## # A tibble: 5 x 5
##   .model          RMSE   MAE  MAPE  MASE
##   <chr>          <dbl> <dbl> <dbl> <dbl>
## 1 hw_damped       69.0  62.8  14.7  3.27
## 2 hw              70.4  60.9  15.3  3.17
## 3 stl_ets_damped  74.8  66.6  14.6  3.47
## 4 snaive          96.8  79.5  16.3  4.14
## 5 stl_ets        110.  106.   24.1  5.52
bind_rows(accuracy(fits), accuracy(fits_stl)) |>
  select(.model, RMSE, MASE) |>
  arrange(RMSE)
## # A tibble: 5 x 3
##   .model          RMSE  MASE
##   <chr>          <dbl> <dbl>
## 1 stl_ets         7.23 0.274
## 2 stl_ets_damped  7.30 0.276
## 3 hw              8.33 0.310
## 4 hw_damped       8.34 0.311
## 5 snaive         26.1  1

The Guerrero method gives lambda = -0.10, close to a log transformation, which again reflects the seasonal variation growing with the level.

On the training data the STL approach fits best of all (RMSE 7.23 against 8.33 for Holt-Winters), because the STL decomposition can follow a seasonal pattern that changes shape over time, which the Holt-Winters seasonal component cannot.

On the test set it is a different story. Plain STL plus ETS gives an RMSE of 109.9, which is worse even than the seasonal naive benchmark, while the damped version gives 74.8. Both are beaten by the Holt-Winters models at 70.4 and 69.0.

So the STL approach does not beat my best previous forecasts. The explanation is visible in the forecast plot: ETS applied to the seasonally adjusted series picks up the steep growth of the late 2000s and extends it in a straight line on the transformed scale, which back-transforms into forecasts that climb too quickly. Damping the trend fixes most of that problem but still does not match Holt-Winters. This is also a good reminder that the best in-sample fit is not the best forecaster: the STL model was first on the training data and last on the test data.