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

library(fpp3)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──
## ✔ tibble      3.3.1     ✔ tsibble     1.1.6
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.2     ✔ feasts      0.4.2
## ✔ lubridate   1.9.4     ✔ fable       0.5.0
## ✔ ggplot2     4.0.1
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
pigs <- aus_livestock |>
  filter(Animal == "Pigs", State == "Victoria")

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

report(ses_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
ses_fc <- ses_fit |>
  forecast(h = "4 months")

ses_fc
## # A fable: 4 x 6 [1M]
## # Key:     Animal, State, .model [1]
##   Animal State    .model                                                   Month
##   <fct>  <fct>    <chr>                                                    <mth>
## 1 Pigs   Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Jan
## 2 Pigs   Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Feb
## 3 Pigs   Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Mar
## 4 Pigs   Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"… 2019 Apr
## # ℹ 2 more variables: Count <dist>, .mean <dbl>
autoplot(ses_fc, pigs)

I fit a simple exponential smoothing model (ETS(A,N,N)) to the Victorian pigs series, with estimated \(\alpha = 0.3221\) and initial level \(\ell_0 = 100{,}646.6\), and then produced 4-month forecasts.

1b) Compute a 95% prediction interval for the first forecast using \(\hat{y} \pm 1.96s\) where \(s\) is the standard deviation of the residuals. Compare your interval with the interval produced by R.

# residual sd (single number)
s <- ses_fit |>
  augment() |>
  pull(.resid) |>
  sd(na.rm = TRUE)

# manual 95% PI for FIRST forecast
first_fc <- ses_fc |>
  as_tibble() |>
  slice(1)

manual_interval <- first_fc |>
  transmute(
    Month,
    forecast = .mean,
    lower95 = .mean - 1.96 * s,
    upper95 = .mean + 1.96 * s
  )

manual_interval
## # A tibble: 1 × 4
##      Month forecast lower95 upper95
##      <mth>    <dbl>   <dbl>   <dbl>
## 1 2019 Jan   95187.  76871. 113502.
# 95% PI for FIRST forecast (for comparison)
r_interval <- ses_fc |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  slice(1) |>
  transmute(
    Month,
    forecast = .mean,
    lower95 = `95%_lower`,
    upper95 = `95%_upper`
  )

r_interval
## # A tibble: 1 × 4
##      Month forecast lower95 upper95
##      <mth>    <dbl>   <dbl>   <dbl>
## 1 2019 Jan   95187.  76855. 113518.
# compare widths
manual_width <- manual_interval$upper95 - manual_interval$lower95
r_width <- r_interval$upper95 - r_interval$lower95

manual_width
## [1] 36631.09
r_width
## [1] 36663.54

Manual 95% PI for the first forecast in January 2019 was \([76871.01,\ 113502.10]\) (width \(36631.09\)), which is essentially the same as R’s 95% PI \([76854.79,\ 113518.30]\) (width \(36663.54\)).

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

exports_data <- global_economy |>
  filter(Country == "United States")

autoplot(exports_data, Exports)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

The Exports series shows a clear long-run upward trend from the 1960s to the 2010s, with several noticeable swings/dips (especially around the early 1980s and early 2000s), and no repeating seasonal pattern since the data are annual.

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

exports_fit <- exports_data |>
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

report(exports_fit)
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.9998995 
## 
##   Initial states:
##     l[0]
##  4.76618
## 
##   sigma^2:  0.4002
## 
##      AIC     AICc      BIC 
## 186.3596 186.8040 192.5409
exports_fc <- exports_fit |>
  forecast(h = 10)

autoplot(exports_fc, exports_data)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

Using an ETS(A,N,N) model for Exports, the estimated smoothing parameter was \(\alpha = 0.9999\) (with initial level \(\ell_0 = 4.7662\)), and the forecasts are essentially flat around the most recent level with widening 80%/95% prediction intervals over the forecast horizon.

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

exports_fit |>
  accuracy() |>
  select(RMSE)
## # A tibble: 1 × 1
##    RMSE
##   <dbl>
## 1 0.627

The training RMSE for the ETS(A,N,N) model is \(0.6271\).

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

exports_compare <- exports_data |>
  model(
    ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
    AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
  )

report(exports_compare)
## Warning in report.mdl_df(exports_compare): Model reporting is only supported
## for individual models, so a glance will be shown. To see the report for a
## specific model, use `select()` and `filter()` to identify a single model.
## # A tibble: 2 × 10
##   Country       .model sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE   MAE
##   <fct>         <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United States ANN     0.400   -90.2  186.  187.  193. 0.386 0.972    NA
## 2 United States AAN     0.399   -89.0  188.  189.  198. 0.372 0.887    NA
exports_compare |>
  accuracy() |>
  select(.model, RMSE)
## # A tibble: 2 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ANN    0.627
## 2 AAN    0.615

ETS(A,A,N) fits slightly better than ETS(A,N,N) (training RMSE \(0.6150\) vs \(0.6271\)), but the improvement is small, so the trended model may be preferred only if you believe exports will keep trending rather than staying around the most recent level

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

exports_compare_fc <- exports_compare |>
  forecast(h = 10)

autoplot(exports_compare_fc, exports_data)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

The ETS(A,N,N) forecasts stay roughly flat around the most recent level, while ETS(A,A,N) projects continued growth, so given the long run upward trend in the data (and its slightly lower RMSE), I would choose ETS(A,A,N).

5f. 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_table <- exports_compare |>
  accuracy() |>
  select(.model, RMSE)

first_forecasts <- exports_compare_fc |>
  as_tibble() |>
  group_by(.model) |>
  slice(1) |>
  ungroup()

first_forecasts |>
  left_join(rmse_table, by = ".model") |>
  mutate(
    lower95 = .mean - 1.96 * RMSE,
    upper95 = .mean + 1.96 * RMSE
  ) |>
  select(.model, Year, .mean, lower95, upper95)
## # A tibble: 2 × 5
##   .model  Year .mean lower95 upper95
##   <chr>  <dbl> <dbl>   <dbl>   <dbl>
## 1 AAN     2018  12.1    10.9    13.3
## 2 ANN     2018  11.9    10.7    13.1
exports_compare_fc |>
  hilo(level = 95)
## # A tsibble: 20 x 6 [1Y]
## # Key:       Country, .model [2]
##    Country       .model  Year
##    <fct>         <chr>  <dbl>
##  1 United States ANN     2018
##  2 United States ANN     2019
##  3 United States ANN     2020
##  4 United States ANN     2021
##  5 United States ANN     2022
##  6 United States ANN     2023
##  7 United States ANN     2024
##  8 United States ANN     2025
##  9 United States ANN     2026
## 10 United States ANN     2027
## 11 United States AAN     2018
## 12 United States AAN     2019
## 13 United States AAN     2020
## 14 United States AAN     2021
## 15 United States AAN     2022
## 16 United States AAN     2023
## 17 United States AAN     2024
## 18 United States AAN     2025
## 19 United States AAN     2026
## 20 United States AAN     2027
## # ℹ 3 more variables: Exports <dist>, .mean <dbl>, `95%` <hilo>

For 5f (R Markdown friendly):

Using \(\hat{y}\pm 1.96,\text{RMSE}\), the first-forecast 95% PI for ETS(A,A,N) was \([10.9172,\ 13.3279]\) and for ETS(A,N,N) was \([10.6616,\ 13.1197]\), which are very close to the intervals produced by R.

  1. 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")

lambda_china <- china |>
  features(GDP, guerrero) |>
  pull(lambda_guerrero)

china_fit <- china |>
  model(
    ETS(GDP),
    ETS(GDP ~ trend("Ad")),
    ETS(box_cox(GDP, lambda_china)),
    ETS(box_cox(GDP, lambda_china) ~ trend("Ad"))
  )

report(china_fit)
## Warning in report.mdl_df(china_fit): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 4 × 10
##   Country .model       sigma2 log_lik   AIC  AICc   BIC      MSE     AMSE    MAE
##   <fct>   <chr>         <dbl>   <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>  <dbl>
## 1 China   "ETS(GDP)"  0.0108  -1546.  3102. 3103. 3112. 4.00e+22 1.61e+23 0.0793
## 2 China   "ETS(GDP ~… 0.0113  -1547.  3105. 3107. 3117. 3.98e+22 1.59e+23 0.0803
## 3 China   "ETS(box_c… 0.00143    74.3 -139. -137. -128. 1.33e- 3 3.52e- 3 0.0286
## 4 China   "ETS(box_c… 0.00150    73.4 -135. -133. -123. 1.37e- 3 4.56e- 3 0.0295
china_fc <- china_fit |>
  forecast(h = 20)

autoplot(china_fc, china)

The different ETS options give noticeably different GDP forecasts. The Box–Cox transformed models produce more stable/realistic growth with much smaller errors (lowest MAE), while the untransformed models explode upward and have very large uncertainty, and adding a damped trend slightly tempers the growth compared with the non damped versions.

  1. 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(Quarter, Gas)

gas_fit <- gas |>
  model(
    ETS(Gas),
    ETS(Gas ~ trend("Ad"))
  )

report(gas_fit)
## Warning in report.mdl_df(gas_fit): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 2 × 9
##   .model                     sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>                       <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 "ETS(Gas)"                0.00324   -831. 1681. 1682. 1711.  21.1  32.2 0.0413
## 2 "ETS(Gas ~ trend(\"Ad\")… 0.00329   -832. 1684. 1685. 1718.  21.1  32.0 0.0417
gas_fc <- gas_fit |>
  forecast(h = "5 years")

autoplot(gas_fc, gas)

Multiplicative seasonality is needed because the seasonal swings get larger as the Gas level increases (the variation scales with the series), and while adding a damped trend gives a very similar fit (MAE \(0.0417\) vs \(0.0413\)), it produces more conservative long-run forecasts with wider uncertainty.

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

set.seed(123456)

my_id <- sample(unique(aus_retail$`Series ID`), 1)

myseries <- aus_retail |>
  filter(`Series ID` == my_id)

myseries |>
  distinct(`Series ID`, State, Industry)
## # A tibble: 1 × 3
##   `Series ID` State      Industry         
##   <chr>       <chr>      <chr>            
## 1 A3349562T   Queensland Department stores
autoplot(myseries, Turnover)

Multiplicative seasonality is necessary because the size of the seasonal spikes increases as turnover increases. The seasonal variation grows proportionally with the level of the series.

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

retail_fit <- myseries |>
  model(
    HW = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

report(retail_fit)
## Warning in report.mdl_df(retail_fit): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 2 × 11
##   State     Industry .model  sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE    MAE
##   <chr>     <chr>    <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Queensla… Departm… HW     0.00242  -2295. 4623. 4625. 4693.  90.2  90.7 0.0378
## 2 Queensla… Departm… HW_da… 0.00246  -2296. 4628. 4630. 4702.  88.9  89.7 0.0380
retail_fc <- retail_fit |>
  forecast(h = "3 years")

autoplot(retail_fc, myseries)

Holt–Winters’ multiplicative model and its damped-trend version give similar forecasts, but the damped version is more conservative, slightly reduces in sample error (MSE \(88.86\) vs \(90.22\)), and avoids pushing the trend upward too aggressively.

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

retail_fit |>
  accuracy() |>
  select(.model, RMSE)
## # A tibble: 2 × 2
##   .model     RMSE
##   <chr>     <dbl>
## 1 HW         9.50
## 2 HW_damped  9.43

The damped Holt–Winters model performs slightly better, with a lower one-step RMSE (\(9.4267\)) than the non-damped version (\(9.4985\)), so I prefer HW_damped.

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

retail_fit |>
  select(HW_damped) |>
  gg_tsresiduals()
## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The residuals from the best model look close to white noise: they are centered around zero with roughly constant variance, the histogram is approximately normal, and most ACF spikes stay within the significance bounds. Only a few small spikes.

8e) Now find the test 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)

test <- myseries |>
  filter(year(Month) >= 2011)

retail_test_fit <- train |>
  model(
    HW = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    SNAIVE = SNAIVE(Turnover)
  )

retail_test_fc <- retail_test_fit |>
  forecast(new_data = test)

retail_test_fc |>
  accuracy(myseries) |>
  select(.model, RMSE)
## # A tibble: 3 × 2
##   .model     RMSE
##   <chr>     <dbl>
## 1 HW         21.5
## 2 HW_damped  18.3
## 3 SNAIVE     12.6

Training through the end of 2010, the seasonal naïve method performs best (test RMSE \(12.6407\)), so neither Holt–Winters model beats SNAIVE (HW: \(21.4596\), HW_damped: \(18.2519\)).

  1. 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 from the test set?
lambda_retail <- train |>
  features(Turnover, guerrero) |>
  pull(lambda_guerrero)

retail_stl_fit <- train |>
  model(
    STL_ETS = decomposition_model(
      STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic")),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))
    )
  )

retail_stl_fc <- retail_stl_fit |>
  forecast(new_data = test)

retail_stl_fc |>
  accuracy(myseries) |>
  select(.model, RMSE)
## # A tibble: 1 × 2
##   .model   RMSE
##   <chr>   <dbl>
## 1 STL_ETS  13.4
autoplot(retail_stl_fc, myseries)

The STL + Box–Cox + ETS approach gives a test RMSE of \(13.3660\), which is worse than the seasonal naïve benchmark from 8e (RMSE \(12.6407\)), so it does not improve on your best test-set forecast.