Forecasting: Principles & Practice - Chapter 8 Exercises: 8.1, 8.5, 8.6, 8.7, 8.8, 8.9

  1. Consider the 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 ℓ0, and generate forecasts for the next four months.
# Pigs slaughtered in Victoria
pigs_vic <- aus_livestock %>%
  filter(State == "Victoria", Animal == "Pigs")

autoplot(pigs_vic, Count)

# Simple exponential smoothing equivalent: ETS(A,N,N)
fit_pigs <- pigs_vic %>%
  model(ets = 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
# Forecast next 4 months
fc_pigs <- fit_pigs %>%
  forecast(h = "4 months")

fc_pigs
autoplot(fc_pigs, pigs_vic)

  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.
#standard deviation of the residuals
aug_pigs <- augment(fit_pigs)

s <- sd(aug_pigs$.resid, na.rm = TRUE)
s
## [1] 9344.666
#first forecast value
first_forecast <- fc_pigs %>%
  as_tibble() %>%
  slice(1)

first_forecast
#manual 95% prediction interval
lower <- first_forecast$.mean - 1.96 * s
upper <- first_forecast$.mean + 1.96 * s

lower
## [1] 76871.01
upper
## [1] 113502.1
#Compare with R's interval
hilo(fc_pigs, 95)
  1. Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
  1. Plot the Exports series and discuss the main features of the data.
#Country: United States
us_exports <- global_economy %>%
  filter(Country == "United States")

autoplot(us_exports, Exports) +
  labs(
    title = "Exports for the United States",
    y = "Exports (% of GDP)",
    x = "Year")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

The exports series for the United States shows fluctuations over time with noticeable increases and decreases across different years. Overall, the series shows variability with possible trend changes over time.

  1. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# Country: United States
us_exports <- global_economy %>%
  filter(Country == "United States")

# ETS(A,N,N) model
fit_exports <- us_exports %>%
  model(ETS(Exports ~ error("A") + trend("N") + season("N")))

report(fit_exports)
## 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
# Generate forecasts
fc_exports <- fit_exports %>%
  forecast(h = 10)

# Plot forecasts
autoplot(us_exports, Exports) +
  autolayer(fc_exports) +
  labs(
    title = "ETS(A,N,N) Forecast of U.S. Exports",
    y = "Exports (% of GDP)",
    x = "Year")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

  1. Compute the RMSE values for the training data.
accuracy(fit_exports)
accuracy(fit_exports) %>%
  select(.model, RMSE)
  1. 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.
# Fit both models
fit_compare <- us_exports %>%
  model(
    ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
    AAN = ETS(Exports ~ error("A") + trend("A") + season("N")))

# Compare model accuracy
accuracy(fit_compare) %>%
  select(.model, RMSE)
report(fit_compare)
## Warning in report.mdl_df(fit_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.
  1. Compare the forecasts from both methods. Which do you think is best?
fc_compare <- fit_compare %>%
  forecast(h = 10)

autoplot(us_exports, Exports) +
  autolayer(fc_compare) +
  labs(
    title = "Forecast Comparison: ETS(A,N,N) vs ETS(A,A,N)",
    y = "Exports (% of GDP)",
    x = "Year"
  )
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

The forecasts from the ETS(A,N,N) model remain relatively flat because the model assumes no trend in the data. In contrast, the ETS(A,A,N) model includes a trend component, so its forecasts continue to increase over time.The ETS(A,A,N) model provides the more realistic forecast.

  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
acc <- accuracy(fit_compare)

rmse_ann <- acc %>%
  filter(.model == "ANN") %>%
  pull(RMSE)

rmse_aan <- acc %>%
  filter(.model == "AAN") %>%
  pull(RMSE)

rmse_ann
## [1] 0.6270672
rmse_aan
## [1] 0.6149566
#First forecast from each model
fc_tbl <- fc_compare %>%
  as_tibble()

first_ann <- fc_tbl %>%
  filter(.model == "ANN") %>%
  slice(1)

first_aan <- fc_tbl %>%
  filter(.model == "AAN") %>%
  slice(1)

first_ann
first_aan
#Manual 95% prediction intervals
ann_lower <- first_ann$.mean - 1.96 * rmse_ann
ann_upper <- first_ann$.mean + 1.96 * rmse_ann

aan_lower <- first_aan$.mean - 1.96 * rmse_aan
aan_upper <- first_aan$.mean + 1.96 * rmse_aan

ann_lower
## [1] 10.66163
ann_upper
## [1] 13.11974
aan_lower
## [1] 10.91722
aan_upper
## [1] 13.32785
#Compare with R’s intervals
hilo(fc_compare, 95)
  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.

[Hint: use a relatively large value of h when forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]

#China GDP
china_gdp <- global_economy %>%
  filter(Country == "China")

autoplot(china_gdp, GDP) +
  labs(title = "China GDP", y = "GDP", x = "Year")

#Various options in the ETS() function
fit_china <- china_gdp %>%
  model(
    default = ETS(GDP),
    trend = ETS(GDP ~ error("A") + trend("A") + season("N")),
    damped = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    boxcox = ETS(box_cox(GDP, 0))
  )

report(fit_china)
## Warning in report.mdl_df(fit_china): 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.
#Forecast
fc_china <- fit_china %>%
  forecast(h = 30)

autoplot(fc_china, china_gdp) +
  labs(title = "China GDP Forecast Comparison")

China’s GDP shows a strong upward trend over time. The standard ETS model and the additive trend model both produce forecasts that continue increasing. When a damped trend is used, the forecasts still increase but the growth rate gradually slows over time.

  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_data <- aus_production %>%
  select(Quarter, Gas)

autoplot(gas_data, Gas) +
  labs(title = "Australian Gas Production", y = "Gas Production")

#ETS Model
fit_gas <- gas_data %>%
  model(ETS(Gas))

report(fit_gas)
## 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
#Forecast next few years
fc_gas <- fit_gas %>%
  forecast(h = "5 years")

autoplot(fc_gas, gas_data) +
  labs(title = "Gas Production Forecast")

#Damped Trend
fit_gas_damped <- gas_data %>%
  model(ETS(Gas ~ error("M") + trend("Ad") + season("M")))

report(fit_gas_damped)
## Series: Gas 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.6489044 
##     beta  = 0.1551275 
##     gamma = 0.09369372 
##     phi   = 0.98 
## 
##   Initial states:
##      l[0]       b[0]      s[0]    s[-1]   s[-2]     s[-3]
##  5.858941 0.09944006 0.9281912 1.177903 1.07678 0.8171255
## 
##   sigma^2:  0.0033
## 
##      AIC     AICc      BIC 
## 1684.028 1685.091 1717.873
fc_gas_damped <- fit_gas_damped %>%
  forecast(h = "5 years")

autoplot(gas_data, Gas) +
  autolayer(fc_gas, series = "ETS") +
  autolayer(fc_gas_damped, series = "Damped ETS")
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Warning in ggdist::geom_lineribbon(without(intvl_mapping, "colour_ramp"), :
## Ignoring unknown parameters: `series`
## Warning in geom_line(mapping = without(mapping, "shape"), data =
## unpack_data(object[single_row[["FALSE"]], : Ignoring unknown parameters:
## `series`
## Scale for fill_ramp is already present.
## Adding another scale for fill_ramp, which will replace the existing scale.

accuracy(fit_gas, fit_gas_damped)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `fit = map(fit, accuracy, measures = measures, ...)`.
## Caused by warning:
## ! 2875 errors (11 unique) encountered
## [2862] 'what' must be a function or character string
## [1] 5 arguments passed to 'gamma' which requires 1
## [1] could not find function "A"
## [1] could not find function "b"
## [1] could not find function "l"
## [2] could not find function "M"
## [1] could not find function "phi"
## [1] could not find function "s0"
## [1] could not find function "s1"
## [1] could not find function "s2"
## [3] unused arguments (.resid = c(0.0887019577390822, -0.635541451706318, 0.181848944105302, 0.4958404364164, -0.201447712972539, 0.330052059101432, -0.76574381018493, 0.20718614708589, -0.181328812837202, 0.334202457614468, 0.404194679517922, -0.437227100352282, -0.335092537654099, 0.24175567746742, 0.41838316814583, -0.317993705252993, 0.71810434631144, 0.151867354090209, -0.904131062635319, 0.394648349454513, -0.120430796419829, -1.09249717469163, 0.0684944391226265, -0.348482592462237, 0.745664155380026, 0.455300309540561, 
## -0.726713604027282, 0.424544674451222, -0.14628453727688, 0.0994223209241891, 0.30248813612535, -0.27556928420397, -0.293791850289678, 0.0506642640809609, 0.322673217664057, -0.174841903610843, -0.231470531453827, 0.0599954937442195, 0.323013100633004, -0.113474316689153, -0.19645600245198, 0.0456758575526601, 1.29972072378184, -0.70628718534188, -0.385055487817186, 0.942645548194687, 0.286209610044773, -0.786259380470792, -0.40459598340701, 0.829646684736847, 1.38419482850271, -0.235498119033231, 
## -0.277924383497541, 0.184127492557002, 1.86893863560563, 1.37672930278689, 2.24342453393183, 1.51863250511909, 1.99788933377637, 1.95681827544784, 0.39128806747463, -4.52721026667563, -0.985780096448998, 1.69504720661791, 2.23107548265966, 1.22617839113938, -0.40252683270937, 2.78277873814032, 1.99859423628947, -5.87031685332145, -2.83999428453362, 0.015200382963144, 2.31016725090544, -2.31616157573146, -5.41393889966384, 0.551432475100718, 0.0868993649890228, -2.14900060003532, -2.19394425188715, 
## 3.07023776899567, 2.40972782066834, 1.36097020617685, -1.9805163808668, 1.00597966452512, 0.183372845307801, 0.0529409411833868, 0.162822256209125, -1.96958340089971, 3.00427570719002, -0.848212538041068, -4.21728426259169, 1.26890807348039, 0.224595643565216, -1.97513451533737, 2.95129291903953, 0.565293815251266, 3.57164858529192, 4.53549839813883, -1.15956358003709, -1.36494571142734, 11.1873432275973, -0.631015322861458, -8.57662832627588, 0.654389057624144, -1.51473294542249, -3.40110941732118, 
## -6.59823860371469, -6.73939302889977, 0.266162219601014, 3.14206737868689, -7.40947031982429, 2.95807019638208, 2.24082131452113, -4.86963933503505, -0.340106900620697, 6.14170075663216, 2.38092583761541, -8.14545793813869, -4.72598656942191, -0.999450609478657, 0.339184538341414, 7.07495741180307, -2.35503414095885, -0.196176820186494, 3.31142359198716, -0.428091496311538, 13.2915528214586, -3.96663783563802, 0.398139769186827, -7.59593204317244, 0.318993188402089, -1.83461877600533, 0.899773299675758, 
## 5.37328455482111, 8.35237074691531, -1.31699787649626, -5.52550593202756, -5.23289284249267, -12.3734670635969, -7.96729796373239, -2.80493382870131, 6.06192540178506, 6.69078974233304, -0.00790072227259486, 2.47019284121299, 7.81075303238779, 8.53984330779085, -9.1996306956556, -6.34423955035732, 3.20118747414639, -1.77697991535811, 10.0694487842319, 1.73539371817458, -0.28215874554067, 5.40700704665116, 1.01706523927433, -5.53146995883924, 0.00774288554879377, -6.81506524398299, -0.139583123077813, 
## 0.0774541032258469, -4.82717082396425, 2.30020331280573, 0.8455765962307, 1.65254423634909, 2.8055413911427, -7.78529711699835, 4.11342851615498, 3.74162479559638, -1.58073555393344, -3.84923478508307, 3.29953253924, 1.62486543199319, 1.0563203550675, -5.01810920765772, 4.47424448808093, 13.4208321586183, -2.97229514514564, -11.9775652153685, -3.34092743442531, 15.797247143719, -10.1938017010449, -14.6201900135658, 4.31085530514204, -5.81243154299312, -0.280762311602388, 3.11111088464978, 0.604968652209095, 
## -0.506721351387171, -6.35775045057756, 5.0058708727218, 0.306692775167392, -1.09622579524995, 4.14623776020179, -10.5104777486533, 3.18557198989947, -4.51073233622543, -2.13411314815809, -5.39236133625883, -0.447826201334266, 6.78150147615685, 18.9965307816004, -8.24086416463936, -1.91652292510784, 8.45465516631648, 6.26908898753041, -9.9887457066684, -0.360227889620603, 2.08001221692516, -9.42658794950387, -2.25579228078391, -2.99106138163407, 6.31496873890768, 1.13646453521167, -6.28195899468415, 
## 0.560142761192168, 8.61910474510941, -11.0602888683393), .train = c(5, 6, 7, 6, 5, 7, 7, 6, 5, 7, 8, 6, 5, 7, 8, 6, 6, 8, 8, 7, 6, 7, 8, 6, 6, 8, 8, 7, 6, 8, 9, 7, 6, 8, 9, 7, 6, 8, 9, 7, 6, 8, 10, 7, 6, 9, 10, 7, 6, 9, 11, 8, 7, 10, 13, 11, 12, 18, 23, 20, 19, 23, 28, 24, 24, 34, 40, 35, 34, 41, 48, 39, 38, 48, 52, 43, 39, 49, 55, 47, 44, 58, 65, 54, 49, 64, 74, 58, 56, 71, 78, 65, 59, 74, 88, 71, 68, 91, 103, 83, 88, 110, 120, 101, 93, 116, 127, 98, 92, 118, 124, 104, 97, 116, 130, 111, 103, 120, 
## 132, 107, 98, 127, 137, 112, 106, 130, 158, 123, 116, 137, 157, 125, 117, 149, 175, 139, 125, 152, 161, 123, 111, 141, 160, 125, 117, 151, 175, 129, 116, 149, 163, 138, 127, 159, 184, 147, 131, 167, 181, 145, 133, 162, 184, 146, 135, 171, 183, 151, 141, 174, 191, 157, 145, 182, 198, 165, 164, 199, 213, 173, 177, 205, 218, 185, 166, 204, 228, 186, 172, 204, 232, 188, 173, 215, 227, 190, 170, 206, 221, 180, 171, 224, 233, 192, 187, 234, 245, 205, 194, 229, 249, 203, 196, 238, 252, 210, 205, 236), .actual = c(5, 
## 6, 7, 6, 5, 7, 7, 6, 5, 7, 8, 6, 5, 7, 8, 6, 6, 8, 8, 7, 6, 7, 8, 6, 6, 8, 8, 7, 6, 8, 9, 7, 6, 8, 9, 7, 6, 8, 9, 7, 6, 8, 10, 7, 6, 9, 10, 7, 6, 9, 11, 8, 7, 10, 13, 11, 12, 18, 23, 20, 19, 23, 28, 24, 24, 34, 40, 35, 34, 41, 48, 39, 38, 48, 52, 43, 39, 49, 55, 47, 44, 58, 65, 54, 49, 64, 74, 58, 56, 71, 78, 65, 59, 74, 88, 71, 68, 91, 103, 83, 88, 110, 120, 101, 93, 116, 127, 98, 92, 118, 124, 104, 97, 116, 130, 111, 103, 120, 132, 107, 98, 127, 137, 112, 106, 130, 158, 123, 116, 137, 157, 125, 
## 117, 149, 175, 139, 125, 152, 161, 123, 111, 141, 160, 125, 117, 151, 175, 129, 116, 149, 163, 138, 127, 159, 184, 147, 131, 167, 181, 145, 133, 162, 184, 146, 135, 171, 183, 151, 141, 174, 191, 157, 145, 182, 198, 165, 164, 199, 213, 173, 177, 205, 218, 185, 166, 204, 228, 186, 172, 204, 232, 188, 173, 215, 227, 190, 170, 206, 221, 180, 171, 224, 233, 192, 187, 234, 245, 205, 194, 229, 249, 203, 196, 238, 252, 210, 205, 236), .fc = c(4.91129804226092, 6.63554145170632, 6.8181510558947, 5.5041595635836, 
## 5.20144771297254, 6.66994794089857, 7.76574381018493, 5.79281385291411, 5.1813288128372, 6.66579754238553, 7.59580532048208, 6.43722710035228, 5.3350925376541, 6.75824432253258, 7.58161683185417, 6.31799370525299, 5.28189565368856, 7.84813264590979, 8.90413106263532, 6.60535165054549, 6.12043079641983, 8.09249717469163, 7.93150556087737, 6.34848259246224, 5.25433584461997, 7.54469969045944, 8.72671360402728, 6.57545532554878, 6.14628453727688, 7.90057767907581, 8.69751186387465, 7.27556928420397, 
## 6.29379185028968, 7.94933573591904, 8.67732678233594, 7.17484190361084, 6.23147053145383, 7.94000450625578, 8.676986899367, 7.11347431668915, 6.19645600245198, 7.95432414244734, 8.70027927621816, 7.70628718534188, 6.38505548781719, 8.05735445180531, 9.71379038995523, 7.78625938047079, 6.40459598340701, 8.17035331526315, 9.61580517149729, 8.23549811903323, 7.27792438349754, 9.815872507443, 11.1310613643944, 9.62327069721311, 9.75657546606817, 16.4813674948809, 21.0021106662236, 18.0431817245522, 18.6087119325254, 
## 27.5272102666756, 28.985780096449, 22.3049527933821, 21.7689245173403, 32.7738216088606, 40.4025268327094, 32.2172212618597, 32.0014057637105, 46.8703168533214, 50.8399942845336, 38.9847996170369, 35.6898327490946, 50.3161615757315, 57.4139388996638, 42.4485675248993, 38.913100635011, 51.1490006000353, 57.1939442518871, 43.9297622310043, 41.5902721793317, 56.6390297938231, 66.9805163808668, 52.9940203354749, 48.8166271546922, 63.9470590588166, 73.8371777437909, 59.9695834008997, 52.99572429281, 71.8482125380411, 
## 82.2172842625917, 63.7310919265196, 58.7754043564348, 75.9751345153374, 85.0487070809605, 70.4347061847487, 64.4283514147081, 86.4645016018612, 104.159563580037, 84.3649457114273, 76.8126567724027, 110.631015322861, 128.576628326276, 100.345610942376, 94.5147329454225, 119.401109417321, 133.598238603715, 104.7393930289, 91.733837780399, 114.857932621313, 131.409470319824, 101.041929803618, 94.7591786854789, 120.869639335035, 130.340106900621, 104.858299243368, 100.619074162385, 128.145457938139, 
## 136.725986569422, 107.999450609479
  1. Recall your retail time series data (from Exercise 7 in Section 2.10).
  1. Why is multiplicative seasonality necessary for this series?

Multiplicative seasonality is necessary because the seasonal fluctuations increase as the level of the retail series increases. This indicates that the seasonal pattern is proportional to the series level rather than constant, making a multiplicative seasonal model more appropriate.

  1. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
set.seed(1234)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))


#Holt-Winters multiplicative method
fit_retail <- myseries %>%
  model(
    hw_mul = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    hw_mul_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))

#Forecasts
fc_retail <- fit_retail %>%
  forecast(h = "3 years")

autoplot(myseries, Turnover) +
  autolayer(fc_retail) +
  labs(
    title = "Holt-Winters Multiplicative Forecast",
    y = "Retail Turnover",
    x = "Year")

  1. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
accuracy(fit_retail) %>%
  select(.model, RMSE)

I prefer the Holt-Winters multiplicative model (hw_mul) because it indicates better forecasting performance.

  1. Check that the residuals from the best method look like white noise.
gg_tsresiduals(fit_retail %>% select(hw_mul))
## 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.

  1. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?
set.seed(1234)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

#Training data: up to end of 2010
myseries_train <- myseries %>%
  filter(year(Month) <= 2010)

#Test data: after 2010
myseries_test <- myseries %>%
  filter(year(Month) > 2010)

#Fit ETS models on training set
fit_train <- myseries_train %>%
  model(
    hw_mul = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    hw_mul_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    snaive = SNAIVE(Turnover)
  )

# Forecast the test period
fc_test <- fit_train %>%
  forecast(new_data = myseries_test)

# Compare test RMSE
accuracy(fc_test, myseries_test) %>%
  select(.model, RMSE)
  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 on the test set?
set.seed(1234)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

#Training and test split
myseries_train <- myseries %>%
  filter(year(Month) <= 2010)

myseries_test <- myseries %>%
  filter(year(Month) > 2010)

#STL decomposition with Box-Cox transformation + ETS
fit_stl <- myseries_train %>%
  model(
    stl_ets = decomposition_model(
      STL(box_cox(Turnover, 0) ~ season(window = "periodic")),
      ETS(season_adjust ~ error("A") + trend("A") + season("N"))))

#Forecast 
fc_stl <- fit_stl %>%
  forecast(new_data = myseries_test)

#Test RMSE
accuracy(fc_stl, myseries_test) %>%
  select(.model, RMSE)

The STL + Box-Cox + ETS model was applied to the seasonally adjusted retail series and evaluated on the test set. When comparing the RMSE values,the hw_mul model produced better forecasts.