All exercises are from Chapter 8 of Hyndman and Athanasopoulos, Forecasting: Principles and Practice (3rd ed.), on exponential smoothing.

Exercise 8.1

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

pigs <- aus_livestock |>
  filter(Animal == "Pigs", State == "Victoria")
autoplot(pigs, 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.

# ETS(A,N,N) is simple exponential smoothing: additive error, no trend, no season
fit_pigs <- pigs |>
  model(ses = ETS(Count ~ error("A") + trend("N") + season("N")))
report(fit_pigs)
## 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
fc_pigs <- fit_pigs |> forecast(h = 4)
fc_pigs |> as_tibble() |> select(Month, .mean)
## # A tibble: 4 × 2
##      Month  .mean
##      <mth>  <dbl>
## 1 2019 Jan 95187.
## 2 2019 Feb 95187.
## 3 2019 Mar 95187.
## 4 2019 Apr 95187.
fc_pigs |>
  autoplot(filter(pigs, year(Month) >= 2010)) +
  geom_line(aes(y = .fitted), 
            data = filter(augment(fit_pigs), 
            year(Month) >= 2010), 
            colour = "darkorange")

The estimated smoothing parameter is alpha = 0.32, with an initial level of about 100,647. Each level update gives about 32% weight to the current observation and 68% to the previous level, which smooths the month-to-month changes. Simple exponential smoothing has no trend or seasonal component, so its four forecasts all equal the final estimated level, about 95,187. This is the model the exercise asks for; the fitted alpha does not show that the series has no trend.

(b) Compute a 95% prediction interval for the first forecast using y-hat plus or minus 1.96s, where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

# s = standard deviation of the innovation residuals
s <- augment(fit_pigs) |> pull(.innov) |> sd(na.rm = TRUE)
yhat <- fc_pigs |> pull(.mean) |> first()
c(lower = yhat - 1.96 * s, upper = yhat + 1.96 * s)
##     lower     upper 
##  76871.01 113502.10
# interval produced by fable for the same first forecast
fc_pigs |> hilo(level = 95) |> slice(1) |> pull(`95%`)
## <hilo[1]>
## [1] [76854.79, 113518.3]95

The two intervals agree to within about 20 on each end, roughly 76,871 to 113,502 by hand against 76,855 to 113,518 from R. The small gap comes from how the residual standard deviation is estimated: sd() and the model’s own sigma^2 divide by slightly different counts. The interval is wide relative to the forecast of about 95,000 because the monthly counts are noisy and the model only has a level to explain them with.

Exercise 8.5

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

# Australia is the country chosen for this exercise
country <- global_economy |> filter(Country == "Australia")

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

autoplot(country, Exports)

Australian exports rise from about 13% of GDP in 1960 to above 20% by the 2000s, so the long-run direction is up. The path is not smooth: the series peaks around 1980, 2001 and 2009, and each peak is followed by a drop. Since the early 2000s it has bounced between roughly 17 and 23 without a clear continued climb. Annual observations do not show within-year seasonality.

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

fit_ann <- country |>
  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(country)

(c) Compute the RMSE values for the training data.

accuracy(fit_ann) |> select(.model, RMSE)
## # A tibble: 1 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ann     1.15

The training RMSE for ETS(A,N,N) is 1.15 percentage points of GDP. The forecast in (b) is a flat line at the last level with intervals that widen with the horizon.

(d) Compare the results to those from an ETS(A,A,N) model. (Remember that the trended model is using one more parameter than the simple model.) Discuss the merits of the two forecasting methods for this data set.

fit_both <- country |>
  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)
## # A tibble: 2 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ann     1.15
## 2 aan     1.12
glance(fit_both) |> select(.model, AICc)
## # A tibble: 2 × 2
##   .model  AICc
##   <chr>  <dbl>
## 1 ann     258.
## 2 aan     259.
fit_both |> forecast(h = 10) |> autoplot(country, level = NULL)

The trended model reduces the training RMSE from 1.15 to 1.12, but its AICc is slightly higher (259 against 258). The extra trend terms improve the fit by very little, so I prefer the simpler ETS(A,N,N) model, although an AICc gap this small is only a modest preference. Its flat forecast also matches the fluctuations around a broadly similar level since the early 2000s.

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

The ETS(A,N,N) forecast stays flat at the most recent level, while the ETS(A,A,N) forecast keeps climbing at the slope it estimated, reaching about 22 by the end of the horizon, above anything the series has done in the last decade except the 2009 peak. Given how the last twenty years look, I find the flat forecast more believable, and the AICc points the same way, though only modestly. I would pick ETS(A,N,N).

(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.

fc_both <- fit_both |> forecast(h = 1)
rmse_both <- accuracy(fit_both) |> select(.model, RMSE)

fc_both |>
  as_tibble() |>
  left_join(rmse_both, by = ".model") |>
  transmute(.model, by_hand_lower = .mean - 1.96 * RMSE, by_hand_upper = .mean + 1.96 * RMSE)
## # A tibble: 2 × 3
##   .model by_hand_lower by_hand_upper
##   <chr>          <dbl>         <dbl>
## 1 ann             18.4          22.9
## 2 aan             18.6          23.0
fc_both |> hilo(level = 95) |> select(.model, `95%`)
## # A tsibble: 2 x 3 [1Y]
## # Key:       .model [2]
##   .model                  `95%`  Year
##   <chr>                  <hilo> <dbl>
## 1 ann    [18.31970, 22.89462]95  2018
## 2 aan    [18.57028, 23.10700]95  2018

The by-hand intervals using 1.96 times the RMSE are close to those from R: 18.4 to 22.9 against 18.3 to 22.9 for ETS(A,N,N), and 18.6 to 23.0 against 18.6 to 23.1 for ETS(A,A,N). The multiplier 1.96 is the 97.5th percentile of the standard normal distribution, which gives a two-sided 95% interval, and for a one-step forecast the forecast standard deviation of both models is the residual standard deviation, which the RMSE approximates. The small differences come from RMSE and the model’s sigma estimate being computed with slightly different divisors. This shortcut only works for the first forecast, because the true forecast variance grows with the horizon.

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. Try to develop an intuition of what each is doing to the forecasts.

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

lambda_china <- china |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda_china
## [1] -0.03446284
fit_china <- china |>
  model(
    ets_auto      = ETS(GDP),
    holt          = ETS(GDP ~ error("A") + trend("A") + season("N")),
    damped        = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    boxcox        = ETS(box_cox(GDP, lambda_china)),
    boxcox_damped = ETS(box_cox(GDP, lambda_china) ~
                          error("A") + trend("Ad") + season("N"))
  )
fit_china |> select(ets_auto) |> report()
## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9998998 
##     beta  = 0.3119984 
## 
##   Initial states:
##         l[0]       b[0]
##  45713434615 3288256682
## 
##   sigma^2:  0.0108
## 
##      AIC     AICc      BIC 
## 3102.064 3103.218 3112.366
# model form chosen by ETS() for each variant, and the smoothing parameters
fit_china
## # A mable: 1 x 6
## # Key:     Country [1]
##   Country     ets_auto         holt        damped       boxcox boxcox_damped
##   <fct>        <model>      <model>       <model>      <model>       <model>
## 1 China   <ETS(M,A,N)> <ETS(A,A,N)> <ETS(A,Ad,N)> <ETS(A,A,N)> <ETS(A,Ad,N)>
fit_china |>
  tidy() |>
  filter(term %in% c("alpha", "beta", "phi")) |>
  select(.model, term, estimate)
## # A tibble: 12 × 3
##    .model        term  estimate
##    <chr>         <chr>    <dbl>
##  1 ets_auto      alpha   1.000 
##  2 ets_auto      beta    0.312 
##  3 holt          alpha   1.000 
##  4 holt          beta    0.552 
##  5 damped        alpha   1.000 
##  6 damped        beta    0.563 
##  7 damped        phi     0.980 
##  8 boxcox        alpha   1.000 
##  9 boxcox        beta    0.0726
## 10 boxcox_damped alpha   1.000 
## 11 boxcox_damped beta    0.355 
## 12 boxcox_damped phi     0.920
fit_china |> forecast(h = 15) |> autoplot(china, level = NULL)

# forecast for the last year of the horizon, in trillions of US dollars
fit_china |>
  forecast(h = 15) |>
  as_tibble() |>
  filter(Year == max(Year)) |>
  transmute(.model, mean = .mean / 1e12, median = median(GDP) / 1e12)
## # A tibble: 5 × 3
##   .model         mean median
##   <chr>         <dbl>  <dbl>
## 1 ets_auto       23.5   23.5
## 2 holt           23.6   23.6
## 3 damped         21.9   21.9
## 4 boxcox        101.    81.3
## 5 boxcox_damped  36.3   20.7

The automatic search picked ETS(M,A,N) with an alpha of about 1, so the level follows the data almost exactly and the forecast is a straight line from the last value. Holt’s linear method (additive trend) and the damped version give nearly the same lines, around 2.2 to 2.4e13 by 2032. The damping parameter phi came out at 0.98, which is the upper end of the default estimation range, so damping only bends the line slightly. China’s GDP level is still rising steeply at the end of the series, so the damped and undamped trends end up close together.

The Box-Cox transformation changes the answer a lot. The estimated lambda is -0.034, which is close to a log transformation but is not exactly a log. The Box-Cox model is an ETS(A,A,N) on the transformed scale, and a straight trend there becomes a curved, rapidly increasing forecast after back-transformation. The forecasts reported as means include the bias adjustment fable applies when it back-transforms, so the size of the mean is affected by the forecast uncertainty as well as by the trend. The table shows how much that matters: the 2032 mean is about 101 trillion dollars without damping and 36 trillion with damping, far above the raw-scale Holt forecast of about 24 trillion, and the medians are noticeably lower than the means for both Box-Cox models. Damping reduces the extrapolated trend (its phi is 0.92, below the 0.98 of the damped model on the raw scale), but these numbers alone do not show which model will forecast Chinese GDP best. AICc values cannot be compared between the raw and transformed models because the response variable is on a different scale.

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?

autoplot(aus_production, Gas)

gg_season(aus_production, Gas)

fit_gas <- aus_production |>
  model(
    additive       = ETS(Gas ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
    mult_damped    = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
    auto           = ETS(Gas)
  )
# model form chosen by the automatic search
fit_gas |> select(auto) |> report()
## 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
glance(fit_gas) |> select(.model, AICc)
## # A tibble: 4 × 2
##   .model          AICc
##   <chr>          <dbl>
## 1 additive       1873.
## 2 multiplicative 1682.
## 3 mult_damped    1685.
## 4 auto           1682.
accuracy(fit_gas) |> select(.model, RMSE, MAE)
## # A tibble: 4 × 3
##   .model          RMSE   MAE
##   <chr>          <dbl> <dbl>
## 1 additive        4.76  3.35
## 2 multiplicative  4.60  3.02
## 3 mult_damped     4.59  3.03
## 4 auto            4.60  3.02
# plot only recent history so the differences between the forecasts are visible
fit_gas |>
  forecast(h = "5 years") |>
  autoplot(filter(aus_production, year(Quarter) >= 1995), level = NULL)

Gas production (in petajoules) starts near 5 and ends above 200, and the size of the seasonal swing grows along with it: the early wiggles are barely visible and the recent ones span about 50 units. In gg_season() the year lines fan out as the level rises. When the seasonal swing scales with the level, additive seasonality (a fixed number of units) cannot fit both ends of the series, and multiplicative seasonality (a percentage of the level) can. The AICc agrees: 1873 for additive against 1682 for multiplicative. This comparison changes the error form as well as the seasonal form, so it supports the full ETS(M,A,M) specification rather than isolating seasonality alone. The automatic search also settled on ETS(M,A,M).

The undamped multiplicative model has the lower AICc, about 1682 against 1685 for the damped one. The training RMSE values, 4.60 and 4.59, are almost identical, so that small gain in fit does not outweigh the AICc penalty, and I prefer the undamped model on this evidence. The forecasts have similar seasonal timing and levels, with damping gradually lowering the projected trend.

Exercise 8.8

Recall your retail time series data (from Exercise 7 in Section 2.10). Turnover is measured in millions of Australian dollars.

set.seed(219562)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> distinct(State, Industry)
## # A tibble: 1 × 2
##   State      Industry                                     
##   <chr>      <chr>                                        
## 1 Queensland Cafes, restaurants and takeaway food services
autoplot(myseries, Turnover)

(a) Why is multiplicative seasonality necessary for this series?

gg_season(myseries, Turnover)

In the time plot, the up-and-down swings from month to month widen as turnover grows. The season plot shows it more clearly: the gap between the lowest month and the highest month is small in the early years and gets much larger in the recent ones, while the shape (a dip in February and a high in December) stays about the same. The seasonal effect scales with the level of the series, which is what multiplicative seasonality describes.

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

fit_retail <- myseries |>
  model(
    hw        = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )
report(fit_retail |> select(hw))
## Series: Turnover 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.6579904 
##     beta  = 0.01278694 
##     gamma = 0.1657436 
## 
##   Initial states:
##     l[0]      b[0]     s[0]     s[-1]    s[-2]   s[-3]     s[-4]    s[-5]
##  47.3486 0.3464185 0.993717 0.9543249 1.077358 1.13534 0.9889655 1.000461
##      s[-6]     s[-7]     s[-8]     s[-9]    s[-10]    s[-11]
##  0.9826625 0.9954632 0.9832827 0.9613547 0.9771835 0.9498865
## 
##   sigma^2:  0.0018
## 
##      AIC     AICc      BIC 
## 4824.534 4825.981 4894.048
report(fit_retail |> select(hw_damped))
## Series: Turnover 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.6655221 
##     beta  = 0.02105753 
##     gamma = 0.0001069884 
##     phi   = 0.9799988 
## 
##   Initial states:
##      l[0]      b[0]      s[0]     s[-1]    s[-2]   s[-3]    s[-4]   s[-5]
##  46.54249 0.1436165 0.9788615 0.9039255 1.031767 1.10979 1.004905 1.03803
##     s[-6]   s[-7]    s[-8]     s[-9]    s[-10]    s[-11]
##  1.001666 1.01403 1.015002 0.9586309 0.9804664 0.9629266
## 
##   sigma^2:  0.0017
## 
##      AIC     AICc      BIC 
## 4811.053 4812.674 4884.656
fit_retail |> forecast(h = "3 years") |> autoplot(myseries, level = NULL)

The undamped model has alpha = 0.66, beta = 0.013 and gamma = 0.17. The damped model has a phi near 0.98, the upper end of the default estimation range, so damping is mild each month but builds up over the forecast horizon. Its forecasts generally show less upward growth, although they are not below the undamped forecasts in every month, because the other smoothing parameters and the seasonal states are re-estimated as well.

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

accuracy(fit_retail) |> select(.model, RMSE)
## # A tibble: 2 × 2
##   .model     RMSE
##   <chr>     <dbl>
## 1 hw         13.4
## 2 hw_damped  13.1
glance(fit_retail) |> select(.model, AICc)
## # A tibble: 2 × 2
##   .model     AICc
##   <chr>     <dbl>
## 1 hw        4826.
## 2 hw_damped 4813.

The damped model has the lower one-step RMSE, 13.1 against 13.4, and the lower AICc, 4813 against 4826. The RMSE gap is only about 2%, but both measures point the same way, so on training fit I prefer the damped model. The test set in (e) tells a different story.

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

fit_retail |> select(hw_damped) |> gg_tsresiduals()

augment(fit_retail) |>
  features(.innov, ljung_box, lag = 24) |>
  select(.model, lb_stat, lb_pvalue) |>
  knitr::kable(digits = c(0, 1, 6))
.model lb_stat lb_pvalue
hw 47.0 0.003317
hw_damped 65.7 0.000010

The damped model’s residuals do not look like white noise. The Ljung-Box p-value is far below 0.05, and the ACF has several bars outside the blue limits. The time plot has one large negative residual in late 1988, and the histogram has a long left tail from it. The undamped model fails the same test with a larger p-value, though it still rejects white noise. The remaining autocorrelation means the models have not captured all of the predictable structure, so prediction intervals based on their assumptions may be unreliable.

(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 through 2010, test on everything after, as the exercise specifies
train <- myseries |> filter(year(Month) <= 2010)
test <- anti_join(myseries, train, by = c("State", "Industry", "Month"))
fit_train <- 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_test <- fit_train |> forecast(new_data = test)
accuracy(fc_test, myseries) |> select(.model, ME, RMSE, MAE)
## # A tibble: 3 × 4
##   .model       ME  RMSE   MAE
##   <chr>     <dbl> <dbl> <dbl>
## 1 hw         96.4  111.  97.8
## 2 hw_damped 141.   158. 142. 
## 3 snaive    163.   180. 164.

The snaive row is the seasonal naive model from Exercise 7 of Section 5.11, fitted to the same training data through 2010 and scored on the same test set. Both Holt-Winters models beat seasonal naive on the common 2011 to 2018 test set. Undamped Holt-Winters has the lower test RMSE, about 111 against 158 for damped Holt-Winters and 180 for seasonal naive. The positive ME means all three under-forecast on average. The comparison in (c) used models fitted to the full series and one-step errors, which favoured damping; this evaluation uses models fitted only through 2010 and an eight-year forecast horizon. A better one-step fit on the full series therefore does not guarantee better forecasts here.

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_retail <- train |> features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda_retail
## [1] 0.08563282
# the same STL step with two choices for the seasonally adjusted part: a damped
# trend, and whatever ETS() picks on its own. season(window = "periodic") fixes the
# seasonal pattern on the transformed scale, and robust = TRUE limits the influence of
# unusual observations on the STL fit
fit_stl <- train |>
  model(
    stl_ets_damped = decomposition_model(
      STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic"), robust = TRUE),
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
    ),
    stl_ets_auto = decomposition_model(
      STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic"), robust = TRUE),
      ETS(season_adjust)
    )
  )

fit_stl |> select(stl_ets_auto) |> report()
## Series: Turnover 
## Model: STL decomposition model 
## Transformation: box_cox(Turnover, lambda_retail) 
## Combination: ETS(A,A,N) + SNAIVE
## 
## ================================
## 
## Series: season_adjust 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.6926966 
##     beta  = 0.0001000008 
## 
##   Initial states:
##      l[0]       b[0]
##  4.555833 0.01061611
## 
##   sigma^2:  0.0047
## 
##      AIC     AICc      BIC 
## 173.3643 173.5413 192.5820 
## 
## Series: season_year 
## Model: SNAIVE 
## 
## sigma^2: 0
fc_stl <- fit_stl |> forecast(new_data = test)
bind_rows(
  accuracy(fc_stl, myseries),
  accuracy(fc_test, myseries)
) |>
  select(.model, RMSE, MAE) |>
  knitr::kable(digits = 1)
.model RMSE MAE
stl_ets_auto 102.6 75.0
stl_ets_damped 143.9 129.5
hw 111.5 97.8
hw_damped 158.0 142.4
snaive 180.0 163.8

The Box-Cox lambda is about 0.09. STL on the transformed series separates the seasonal pattern, and ETS then models the seasonally adjusted part. The seasonal component is forecast by seasonal naive, as the decomposition report above shows. I fitted two versions of the model for the seasonally adjusted part: a damped additive trend that I specified myself, and the model ETS() selects on its own, which is ETS(A,A,N), an undamped additive trend.

The damped version has a test RMSE of 144, better than damped Holt-Winters (158) and seasonal naive (180) but worse than undamped Holt-Winters (111). That repeats the pattern from 8.8: the damped version has the higher RMSE over this test period. The automatically selected version has a test RMSE of 103, the lowest of the five models, about 8% below undamped Holt-Winters. So decomposing first and letting ETS keep the full trend does give the best test-set forecast here, though the margin over plain Holt-Winters is small compared with the gap to seasonal naive.