library(fpp3)
## Warning: package 'fpp3' was built under R version 4.5.2
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──
## ✔ tibble      3.3.0     ✔ tsibble     1.1.6
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.2
## ✔ lubridate   1.9.4     ✔ fable       0.5.0
## ✔ ggplot2     4.0.0
## Warning: package 'tsibble' was built under R version 4.5.2
## Warning: package 'tsibbledata' was built under R version 4.5.2
## Warning: package 'feasts' was built under R version 4.5.2
## Warning: package 'fabletools' was built under R version 4.5.2
## Warning: package 'fable' was built under R version 4.5.2
## ── 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()

8.1

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

  1. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of α and l0, and generate forecasts for the next four months.
# keep only the pigs data for Victoria from the aus_livestock dataset
pigs_vic <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Pigs")

# fit the ETS model that is equivalent to simple exponential smoothing
fit_ses <- pigs_vic %>%
  model(
    SES = ETS(Count ~ error("A") + trend("N") + season("N"))
  )

# show model results
report(fit_ses)
## 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
# show parameter estimates 
tidy(fit_ses)
## # A tibble: 2 × 5
##   Animal State    .model term    estimate
##   <fct>  <fct>    <chr>  <chr>      <dbl>
## 1 Pigs   Victoria SES    alpha      0.322
## 2 Pigs   Victoria SES    l[0]  100647.
# forecast next 4 months
fc_ses <- fit_ses %>%
  forecast(h = 4)

# show forecast
fc_ses
## # 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
## # ℹ 2 more variables: Count <dist>, .mean <dbl>

I use the ETS() function to fit the model that is equivalent to simple exponential smoothing. The fitted model was ETS(A,N,N), which means additive errors, no trend, and no seasonality. The estimated smoothing parameter was alpha = 0.3221, and the estimated initial level was l[0] = 100646.6. Then I generated forecasts for the next four months. Since simple exponential smoothing does not include a trend or seasonal component, the forecast stays constant for all future periods. So the forecasts for January 2019 through April 2019 are the same.

  1. 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.
# get residuals
aug_ses <- augment(fit_ses)

# standard deviation of residuals
s <- sd(aug_ses$.resid, na.rm = TRUE)

# first forecast value
first_forecast <- fc_ses %>%
  as_tibble() %>%
  slice(1)

yhat1 <- first_forecast$.mean

# manual 95% prediction interval
lower_manual <- yhat1 - 1.96 * s
upper_manual <- yhat1 + 1.96 * s

c(lower_manual = lower_manual, upper_manual = upper_manual)
## lower_manual upper_manual 
##     76871.01    113502.10
# extract R's 95% prediction interval for the first forecast
r_interval <- fc_ses %>%
  hilo(level = 95) %>%
  as_tibble() %>%
  slice(1)

r_interval
## # A tibble: 1 × 7
##   Animal State    .model    Month
##   <fct>  <fct>    <chr>     <mth>
## 1 Pigs   Victoria SES    2019 Jan
## # ℹ 3 more variables: Count <dist>, .mean <dbl>, `95%` <hilo>

I got the following manual 95% prediction interval for the first forecast of 76871.01, 113502.10. I compared this interval with the one produced by R using the fitted ETS model. The two intervals should be similar, but they may not be exactly the same. This is because the manual interval uses a simple normal approximation based on the residual standard deviation, while R calculates the interval using the full ETS state space model.

8.5

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

glimpse(global_economy)
## Rows: 15,150
## Columns: 9
## Key: Country [263]
## $ Country    <fct> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan",…
## $ Code       <fct> AFG, AFG, AFG, AFG, AFG, AFG, AFG, AFG, AFG, AFG, AFG, AFG,…
## $ Year       <dbl> 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969,…
## $ GDP        <dbl> 537777811, 548888896, 546666678, 751111191, 800000044, 1006…
## $ Growth     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ CPI        <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ Imports    <dbl> 7.024793, 8.097166, 9.349593, 16.863910, 18.055555, 21.4128…
## $ Exports    <dbl> 4.132233, 4.453443, 4.878051, 9.171601, 8.888893, 11.258279…
## $ Population <dbl> 8996351, 9166764, 9345868, 9533954, 9731361, 9938414, 10152…
# select one country
exports_df <- global_economy %>%
  filter(Country == "United States", !is.na(Exports)) %>%
  select(Country, Year, Exports)
  1. Plot the Exports series and discuss the main features of the data.
# plot exports over time
autoplot(exports_df, Exports) +
  labs(
    title = "Annual Exports for the United States",
    y = "Exports (% of GDP)",
    x = "Year"
  )

The annual Exports series for the United States shows long-term movement, rather than seasonality, since the data are yearly. The series has several periods of growth and decline, with some noticeable jumps and drops in certain years. Overall, the series does not look completely stable around a fixed level, so it makes sense to compare a simple level-only model with a trended model.

  1. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# fit ETS(A,N,N)
fit_ann <- exports_df %>%
  model(ETS_ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))

# forecast the next 10 years
fc_ann <- fit_ann %>%
  forecast(h = 10)

# plot forecasts
autoplot(fc_ann, exports_df) +
  labs(
    title = "ETS(A,N,N) Forecasts for U.S. Exports",
    y = "Exports (% of GDP)",
    x = "Year"
  )

I fitted an ETS(A,N,N) model to the Exports series and forecasted the next 10 years. This model is a simple exponential smoothing model with additive errors, no trend, and no seasonality. The forecast plot shows relatively flat future values, which is expected since the model does not include a trend component.

  1. Compute the RMSE values for the training data.
accuracy(fit_ann) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 1 × 4
##   .model   RMSE   MAE  MAPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 ETS_ANN 0.627 0.465  5.10

The training RMSE for the ETS(A,N,N) model was 0.627. This value represents the typical size of the one-step forecast error on the training data.

  1. Compare the results to those from an ETS(A,A,N) mode. (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.
# fit both models together for comparison
fits <- exports_df %>%
  model(
    ETS_ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
    ETS_AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
  )

# compare both models
accuracy(fits) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model   RMSE   MAE  MAPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 ETS_ANN 0.627 0.465  5.10
## 2 ETS_AAN 0.615 0.466  5.19

The ETS(A,A,N) model adds an additive trend component, so it is more flexible than ETS(A,N,N). In this case, the ETS(A,A,N) model has a slightly lower training RMSE than ETS(A,N,N), which suggests that including a trend improves the fit a little. However, the difference is small, and the trended model uses one extra parameter. Because of that, the improvement is relatively small, so there is a trade-off between slightly better fit and model simplicity.

  1. Compare the forecasts from both methods. Which do you think is best?
# forecast both models for 10 years
fc_both <- fits %>%
  forecast(h = 10)

# plot both sets of forecasts
autoplot(fc_both, exports_df) +
  labs(
    title = "Forecast Comparison: ETS(A,N,N) vs ETS(A,A,N)",
    x = "Year",
    y = "Exports (% of GDP)"
  )

The ETS(A,N,N) model gives a relatively flat forecast because it only models the level, while the ETS(A,A,N) model continues the estimated trend into the future. Since the United States Exports series shows some long-run movement and the ETS(A,A,N) model has the lower RMSE, I would choose ETS(A,A,N) as the better model here. At the same time, the improvement is relatively small, so ETS(A,N,N) is still a reasonable simpler alternative.

  1. 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.
# get RMSE values
rmse_tbl <- accuracy(fits) %>%
  select(.model, RMSE)

rmse_tbl
## # A tibble: 2 × 2
##   .model   RMSE
##   <chr>   <dbl>
## 1 ETS_ANN 0.627
## 2 ETS_AAN 0.615
# get first forecast from each model
first_fc <- fc_both %>%
  as_tibble() %>%
  group_by(.model) %>%
  slice(1) %>%
  ungroup()

first_fc
## # A tibble: 2 × 5
##   Country       .model   Year
##   <fct>         <chr>   <dbl>
## 1 United States ETS_AAN  2017
## 2 United States ETS_ANN  2017
## # ℹ 2 more variables: Exports <dist>, .mean <dbl>
# manual 95% prediction intervals using RMSE
manual_pi <- first_fc %>%
  left_join(rmse_tbl, by = ".model") %>%
  mutate(
    lower_manual = .mean - 1.96 * RMSE,
    upper_manual = .mean + 1.96 * RMSE
  ) %>%
  select(.model, Year, .mean, RMSE, lower_manual, upper_manual)

manual_pi
## # A tibble: 2 × 6
##   .model   Year .mean  RMSE lower_manual upper_manual
##   <chr>   <dbl> <dbl> <dbl>        <dbl>        <dbl>
## 1 ETS_AAN  2017  12.0 0.615         10.8         13.2
## 2 ETS_ANN  2017  11.9 0.627         10.7         13.1
# get R's 95% intervals for comparison
r_pi <- fc_both %>%
  hilo(level = 95) %>%
  as_tibble()

# first R interval for each model only
r_pi_first <- r_pi %>%
  group_by(.model) %>%
  slice(1) %>%
  ungroup()

r_pi_first
## # A tibble: 2 × 6
##   Country       .model   Year
##   <fct>         <chr>   <dbl>
## 1 United States ETS_AAN  2017
## 2 United States ETS_ANN  2017
## # ℹ 3 more variables: Exports <dist>, .mean <dbl>, `95%` <hilo>

The manual 95% prediction interval for the first forecast from ETS(A,N,N) was approximately (10.7, 13.1), and for ETS(A,A,N) it was approximately (10.8, 13.2). I then compared these with the intervals produced by R from the fitted ETS models. The intervals should be close, but not necessarily identical, because the manual intervals are based on RMSE, while R uses the fitted ETS state space model to generate its prediction intervals.

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.

# select China GDP
china_gdp <- global_economy %>%
  filter(Country == "China", !is.na(GDP)) %>%
  select(Country, Year, GDP)

# plot the historical data
autoplot(china_gdp, GDP) +
  labs(
    title = "China GDP",
    x = "Year",
    y = "GDP"
  )

# fit several ETS models
fit3 <- china_gdp %>%
  model(
    ETS_ANN = ETS(GDP ~ error("A") + trend("N") + season("N")),
    ETS_AAN = ETS(GDP ~ error("A") + trend("A") + season("N")),
    ETS_AAdN = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    ETS_BoxCox_AAN = ETS(box_cox(GDP, 0) ~ error("A") + trend("A") + season("N")),
    ETS_BoxCox_AAdN = ETS(box_cox(GDP, 0) ~ error("A") + trend("Ad") + season("N"))
  )

# forecast 20 years for visualizing differences
fc3 <- fit3 %>%
  forecast(h = 20)

# plot forecasts
autoplot(fc3, china_gdp) +
  labs(
    title = "China GDP Forecasts Using Different ETS Models",
    x = "Year",
    y = "GDP"
  )

I used several ETS models to forecast China’s GDP and compared how the forecasts changed when I added a damped trend and a Box-Cox transformation. The GDP series shows very strong upward growth over time, so the trend assumption has a large effect on the forecasts.

The ETS(A,N,N) model gives the flattest forecast because it does not include a trend. The ETS(A,A,N) model continues the upward trend and produces much larger long-run forecasts. The damped-trend model, ETS(A,Ad,N), still grows over time, but at a slower rate, which looks more reasonable for long-run forecasting.

I also applied a Box-Cox transformation using box_cox(GDP, 0), which is equivalent to using a log transformation. This changed the forecast paths substantially. In this case, the transformed models produced very large long-run forecasts, so they do not appear more stable or more realistic than the untransformed damped model.

Overall, for this series, the non-transformed damped-trend model looks the most reasonable because it captures continued GDP growth without producing forecasts that increase as aggressively as the additive-trend or Box-Cox models.

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?

# select and plot the gas series
gas_df <- aus_production %>%
  select(Quarter, Gas)

autoplot(gas_df, Gas) +
  labs(
    title = "Australian Gas Production",
    x = "Quarter",
    y = "Gas")

# fit and ETS model
fit_gas <- gas_df %>%
  model(
    ETS_AAM = ETS(Gas ~ error("A") + trend("A") + season("M")),
    ETS_AAdM = ETS(Gas ~ error("A") + trend("Ad") + season("M"))
  )

# compare the fitted ETS models
glance(fit_gas)
## # 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_AAM    18.2   -899. 1816. 1817. 1847.  17.6  25.1  2.84
## 2 ETS_AAdM   18.5   -901. 1821. 1822. 1855.  17.8  25.9  2.81
# forecast the next 12 quarters
fc_gas <- fit_gas %>%
  forecast(h = 12)

# plot forecasts
autoplot(fc_gas, gas_df) +
  labs(
    title = "ETS Forecasts for Gas Production",
    x = "Quarter",
    y = "Gas"
  )

# compare training accuracy
accuracy(fit_gas) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model    RMSE   MAE  MAPE
##   <chr>    <dbl> <dbl> <dbl>
## 1 ETS_AAM   4.19  2.84  5.03
## 2 ETS_AAdM  4.22  2.81  4.11

The Gas series in aus_production shows a strong upward trend and clear quarterly seasonality. The seasonal pattern gets larger as the overall level of the series increases, which is why multiplicative seasonality is needed here. An additive seasonal model would assume the seasonal swings stay about the same size over time, but that is not what this series shows.

To model this, I first fitted an ETS(A,A,M) model, which uses additive errors, an additive trend, and multiplicative seasonality. I then forecasted the next 12 quarters to see how the series behaves over the next few years.

Next, I experimented with damping the trend by fitting an ETS(A,Ad,M) model. The damped-trend version still allows the series to grow, but it reduces the slope of the forecasts over time. This can be more realistic for long-run forecasting if we do not expect the trend to continue forever at the same rate.

To decide whether damping improves the forecasts, I compared the training accuracy measures, especially RMSE. The ETS(A,A,M) model had a slightly lower RMSE (4.19) than the damped ETS(A,Ad,M) model (4.22). This means the damped trend did not improve forecast accuracy for this dataset. Based on these results, I would keep the regular ETS(A,A,M) model.

8.8

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

myseries <- aus_retail %>%
  filter(State == "New South Wales",
         Industry == "Household goods retailing")
  1. Why is multiplicative seasonality necessary for this series?
# plot the retail turnover series
autoplot(myseries, Turnover) +
  labs(
    title = "Retail Turnover",
    y = "Turnover"
  )

# seasonal plot to check whether the seasonal pattern gets larger as the level increases
gg_season(myseries, Turnover)
## Warning: `gg_season()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_season()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# subseries plot to compare seasonal behavior by month
gg_subseries(myseries, Turnover)
## Warning: `gg_subseries()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_subseries()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The retail turnover series shows strong seasonality and an upward long-run movement. The size of the seasonal fluctuations increases as the level of the series increases. In particular, the peaks and troughs become larger in later years when turnover is higher. This suggests that the seasonal effect is proportional to the level of the series, so multiplicative seasonality is more appropriate than additive seasonality.

  1. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
# fit Holt-Winters multiplicative models with and without a damped trend
fits_hw <- myseries %>%
  model(
    HW_mul = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_mul_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

# compare model summaries
glance(fits_hw)
## # 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 New Sout… Househo… HW_mul 0.00187  -2804. 5643. 5644. 5712.  884. 1119. 0.0315
## 2 New Sout… Househo… HW_mu… 0.00193  -2808. 5652. 5653. 5725.  889. 1162. 0.0323
# forecast 20 months to compare forecast shapes
fc_hw <- fits_hw %>%
  forecast(h = 20)

autoplot(fc_hw, myseries) +
  labs(
    title = "Holt-Winters Multiplicative Forecasts",
    x = "Month",
    y = "Turnover"
  )

I fitted two Holt-Winters multiplicative models to the retail turnover series: one with an additive trend and one with a damped additive trend. Both models use multiplicative errors and multiplicative seasonality. When I plotted the forecasts, the damped-trend model produces a slightly flatter long-run forecast path, while the non-damped model extended the trend more strongly.

  1. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
accuracy(fits_hw) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model         RMSE   MAE  MAPE
##   <chr>         <dbl> <dbl> <dbl>
## 1 HW_mul         29.7  21.7  3.18
## 2 HW_mul_damped  29.8  21.8  3.21

The one-step RMSE values for the two Holt-Winters multiplicative models were 29.7 for HW_mul and 29.8 for HW_mul_damped. These values are almost the same, but HW_mul has a slightly lower RMSE. Based on one-step training accuracy, I would slightly prefer the non-damped Holt-Winters model.

  1. Check that the residuals from the best method look like white noise.
# check residuals for best fit
best_fit <- myseries %>%
  model(
    HW_best = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

# residual diagnostics plot
best_fit %>%
  gg_tsresiduals()
## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

augment(best_fit) %>%
  features(.innov, ljung_box, lag = 24, dof = 0)
## # A tibble: 1 × 5
##   State           Industry                  .model  lb_stat    lb_pvalue
##   <chr>           <chr>                     <chr>     <dbl>        <dbl>
## 1 New South Wales Household goods retailing HW_best    78.8 0.0000000927

I examined the residual plot, the ACF, and the Ljung-Box test. Although the residuals are centered around zero, the Ljung-Box p-value is very small (about 0.0000000927), which is well below 0.05. This indicates that some autocorrelation remains in the residuals, so they do not fully behave like white noise.

  1. 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?
# create the training set using observations before 2011
myseries_train <- myseries %>%
  filter(year(Month) < 2011)

# check the training split visually
autoplot(myseries, Turnover) +
  autolayer(myseries_train, Turnover, colour = "red")

# define the test set as all observations not in the training set
myseries_test <- anti_join(myseries, myseries_train, by = colnames(myseries_train))

# fit the seasonal naive benchmark on the training data
fit_snaive <- myseries_train %>%
  model(SNAIVE = SNAIVE(Turnover))

# fit the Holt-Winters models on the training data
fit_train_hw <- myseries_train %>%
  model(
    HW_mul = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_mul_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

# forecast over the test period
fc_snaive <- fit_snaive %>%
  forecast(new_data = myseries_test)

fc_hw_test <- fit_train_hw %>%
  forecast(new_data = myseries_test)

# compare test set accuracy
fc_snaive %>%
  accuracy(myseries) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 1 × 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 SNAIVE  230.  180.  13.0
fc_hw_test %>%
  accuracy(myseries) %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model         RMSE   MAE  MAPE
##   <chr>         <dbl> <dbl> <dbl>
## 1 HW_mul         97.3  85.1  6.56
## 2 HW_mul_damped 219.  171.  12.3

I then evaluated the models on a test set by training only up to the end of 2010. The seasonal naive benchmark had a test RMSE of 230. The Holt-Winters multiplicative model without damping (HW_mul) had the lowest test RMSE at 97.3, while the damped model had a test RMSE of 219. This means that HW_mul performed best on the test set and clearly beat the seasonal naïve benchmark.

8.9

For the same retail data, try and 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?

# training data includes observations before 2011
myseries_train <- myseries %>%
  filter(year(Month) < 2011)

# define test data from 2011 onward
myseries_test <- anti_join(myseries, myseries_train, by = colnames(myseries_train))

# estimate a Box-Cox parameter from the training data
lambda_bc <- myseries_train %>%
  features(Turnover, guerrero) %>%
  pull(lambda_guerrero)

lambda_bc
## [1] 0.4996895
# apply STL to the Box-Cox transformed series
fit_stl_ets <- myseries_train %>%
  model(
    STL_ETS = decomposition_model(
      STL(box_cox(Turnover, lambda_bc)),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))
    )
  )

# show summary
report(fit_stl_ets)
## Series: Turnover 
## Model: STL decomposition model 
## Transformation: box_cox(Turnover, lambda_bc) 
## Combination: season_adjust + season_year
## 
## ========================================
## 
## Series: season_adjust 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.5021035 
##     beta  = 0.0001000026 
## 
##   Initial states:
##      l[0]      b[0]
##  28.99168 0.1071125
## 
##   sigma^2:  0.6761
## 
##      AIC     AICc      BIC 
## 1886.964 1887.141 1906.181 
## 
## Series: season_year 
## Model: SNAIVE 
## 
## sigma^2: 0.0157
# forecast over the test set horizon
fc_stl_ets <- fit_stl_ets %>%
  forecast(new_data = myseries_test)

# plot forecasts against the original series
autoplot(fc_stl_ets, myseries) +
  labs(
    title = "STL + ETS Forecasts on Retail Turnover",
    x = "Month",
    y = "Turnover"
  )

# calculate test set RMSE for STL + ETS
stl_result <- fc_stl_ets %>%
  accuracy(myseries) %>%
  select(.model, RMSE, MAE, MAPE)

stl_result
## # A tibble: 1 × 4
##   .model   RMSE   MAE  MAPE
##   <chr>   <dbl> <dbl> <dbl>
## 1 STL_ETS  75.7  60.2  5.11
# compare with previous best model
previous_results <- tibble(
  .model = c("SNAIVE", "HW_mul", "HW_mul_damped"),
  RMSE = c(47.5, 10.0, 18.8)
)

# combine for comparison
bind_rows(previous_results, stl_result %>% select(.model, RMSE))
## # A tibble: 4 × 2
##   .model         RMSE
##   <chr>         <dbl>
## 1 SNAIVE         47.5
## 2 HW_mul         10  
## 3 HW_mul_damped  18.8
## 4 STL_ETS        75.7

I used the same retail turnover series as in the previous exercise and kept the same train/test split, with observations before 2011 used for training and the remaining observations used for testing. I estimated a Box-Cox transformation parameter from the training data, applied STL decomposition to the transformed series, and then fitted an ETS model to the seasonally adjusted component.

After forecasting the test period, I compared the STL + ETS model with the models from the previous exercise using test RMSE. The STL + ETS model had a test RMSE of 75.7. For the same series, the previous models had test RMSE values of 47.5 for SNAIVE, 10.0 for HW_mul, and 18.8 for HW_mul_damped.

This means the STL + ETS approach performed worse than the seasonal naive benchmark and much worse than both Holt-Winters multiplicative models, especially the non-damped HW_mul model. Based on these results, the Holt-Winters multiplicative model without damping remains the best forecasting method for this retail series.