library(fpp3)        # loads tsibble, tsibbledata, feasts, fable, dplyr, ggplot2, ...

# Same helper as the Chapter 5 homework which takes a fitted model (mable) and runs
# the Ljung-Box test on the innovation residuals 
# lag is 10 for non-seasonal data and lag is 2m for seasonal basically 24 months.
lb_test <- function(fit, lag) {
  augment(fit) |>                    # adds .fitted, .resid and .innov columns
    features(.innov, ljung_box, lag = lag)
}

Exercise 8.1 — Pigs slaughtered in Victoria

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 \(\alpha\) and \(\ell_0\), and generate forecasts for the next four months.

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

(a) Simple exponential smoothing with ETS()

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

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

Simple exponential smoothing is ETS(A,N,N) additive errors, no trend, no seasonality so I can spell them out in the formula just like the book

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

report(fit_pigs)      # shows alpha, l[0], sigma^2 and the AIC values
## 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
# tidy() pulls the estimation onto the table so I can grab them directly
tidy(fit_pigs)
## # 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.
fc_pigs <- fit_pigs |> forecast(h = 4)    # next 4 months
fc_pigs
## # 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>
fc_pigs |>
  autoplot(pigs |> filter(year(Month) >= 2014)) +   # zoom in on recent years
  labs(title = "SES forecasts of Victorian pig slaughter, Jan-Apr 2019",
       y = "Number of pigs", x = "Month")

Answer. The optimal values are \(\alpha \approx 0.322\) and \(\ell_0 \approx 100{,}647\) pigs. .

(b) 95% prediction interval for the first forecast

# s = standard deviation of the residuals as one step training
s <- augment(fit_pigs) |> pull(.resid) |> sd()

y_hat <- fc_pigs$.mean[1]       # first forecast which is in Jan 2019

manual_pi <- c(lower = y_hat - 1.96 * s, upper = y_hat + 1.96 * s)

# R's interval hilo() turns the forecast distribution into an interval
r_pi <- fc_pigs |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>          # splits the interval into lower/upper columns
  slice(1)

tibble(
  method = c("Manual (y_hat +/- 1.96 s)", "R (hilo)"),
  sd_used = c(s, sqrt(glance(fit_pigs)$sigma2)),
  lower = c(manual_pi["lower"], r_pi$`95%_lower`),
  upper = c(manual_pi["upper"], r_pi$`95%_upper`)
)
## # A tibble: 2 × 4
##   method                    sd_used  lower   upper
##   <chr>                       <dbl>  <dbl>   <dbl>
## 1 Manual (y_hat +/- 1.96 s)   9345. 76871. 113502.
## 2 R (hilo)                    9353. 76855. 113518.

Answer. My interval is about [76,871, 113,502] and R’s is about [76,855, 113,518], so they’re about the same, R’s is just nearly 16 pigs wider on each side bu this is because of how the standard deviation is calculated, sd() uses \(T-1\) while ETS() uses \(T-k\)

Exercise 8.5 — Exports for one country (global_economy)

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.
  2. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
  3. Compute the RMSE values for the training data.
  4. 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.
  5. Compare the forecasts from both methods. Which do you think is best?
  6. 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.

I picked Germany.

(a) Plot and main features

germany <- global_economy |>
  filter(Country == "Germany", !is.na(Exports))   # Germany starts in 1970

germany |>
  autoplot(Exports) +
  labs(title = "Germany: exports of goods and services",
       y = "% of GDP", x = "Year")

Answer. Germany’s exports went from nearly 15% of GDP in 1970 to 47% in 2017, so overall there’s a strong upward trend. The data is anual so no seasonality. Growth isn’t smooth as theres dips and then a hard run-up from 1995 to 2008. There’s a drop off in 2009 and then it recovers but stays flat after.

(b) ETS(A,N,N) forecasts

fit_ger <- germany |>
  model(
    ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),   # SES
    AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))    # Holt
  )

fit_ger |>
  select(Country, ANN) |>         # keep the key column, just the SES model
  forecast(h = 10) |>
  autoplot(germany) +
  labs(title = "Germany exports: ETS(A,N,N) forecasts",
       y = "% of GDP", x = "Year")

tidy(fit_ger) |> filter(.model == "ANN")
## # A tibble: 2 × 4
##   Country .model term  estimate
##   <fct>   <chr>  <chr>    <dbl>
## 1 Germany ANN    alpha    1.000
## 2 Germany ANN    l[0]    15.0

(c) Training RMSE

accuracy(fit_ger) |>
  select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 × 5
##   .model  RMSE   MAE  MAPE  MASE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 ANN     1.76  1.31  4.73 0.982
## 2 AAN     1.62  1.11  4.14 0.832

Answer. The training RMSE is about 1.76 percentage points.

(d) Compare with ETS(A,A,N)

tidy(fit_ger) |> filter(.model == "AAN")
## # A tibble: 4 × 4
##   Country .model term   estimate
##   <fct>   <chr>  <chr>     <dbl>
## 1 Germany AAN    alpha  0.965   
## 2 Germany AAN    beta   0.000100
## 3 Germany AAN    l[0]  14.3     
## 4 Germany AAN    b[0]   0.697
# AICc already penalizes the trended model for its extra parameter
glance(fit_ger) |> select(.model, sigma2, AIC, AICc, BIC)
## # A tibble: 2 × 5
##   .model sigma2   AIC  AICc   BIC
##   <chr>   <dbl> <dbl> <dbl> <dbl>
## 1 ANN      3.23  246.  247.  252.
## 2 AAN      2.88  242.  244.  252.
# cross-validation on one-step forecasts, starting with 20 years of data
germany |>
  stretch_tsibble(.init = 20, .step = 1) |>
  model(
    ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
    AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
  ) |>
  forecast(h = 1) |>
  accuracy(germany) |>
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 AAN     2.13  1.47  4.50
## 2 ANN     2.11  1.63  4.90

Answer. The trended model fits a bit better in-sample, RMSE dropped from 1.76 to 1.62, and the ME goes from 0.67 to about 0, so ANN was running low on average because it ignores the trend. But even after being penalized for the extra paramater AIC still chose AAN.

(e) Compare the forecasts

fc_ger <- fit_ger |> forecast(h = 10)

fc_ger |>
  autoplot(germany, level = NULL) +      # level, keeping in mind NULL hides the intervals
  labs(title = "Germany exports: ETS(A,N,N) vs ETS(A,A,N) forecasts",
       y = "% of GDP", x = "Year")

fc_ger |>
  as_tibble() |>
  filter(Year %in% c(2018, 2022, 2027)) |>
  select(.model, Year, .mean) |>
  pivot_wider(names_from = Year, values_from = .mean)
## # A tibble: 2 × 4
##   .model `2018` `2022` `2027`
##   <chr>   <dbl>  <dbl>  <dbl>
## 1 ANN      47.2   47.2   47.2
## 2 AAN      47.9   50.7   54.2

Answer. ANN stays flat at for almost all 10 years while AAN climbs about 0.7 points a year . ANN has better number but for forecasting I would actually go with ANN, the last six years were flat, ANN’s slope is just old growth rate that never updated and exports just can’t keep growing as a share of GDP forever.

(f) 95% prediction intervals for the first forecast

rmse_ger <- accuracy(fit_ger) |> select(.model, RMSE)

#first forecast +- 1.96 * RMSE which are the normal errors
manual_ger <- fc_ger |>
  as_tibble() |>
  filter(Year == 2018) |>
  select(.model, .mean) |>
  left_join(rmse_ger, by = ".model") |>
  mutate(lower = .mean - 1.96 * RMSE,
         upper = .mean + 1.96 * RMSE,
         source = "Manual (RMSE)")

# R's own interval from the forecast distribution
r_ger <- fc_ger |>
  filter(Year == 2018) |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  transmute(.model, .mean, lower = `95%_lower`, upper = `95%_upper`,
            source = "R (hilo)")

bind_rows(manual_ger |> select(-RMSE), r_ger) |>
  arrange(.model, source)
## # A tibble: 4 × 5
##   .model .mean lower upper source       
##   <chr>  <dbl> <dbl> <dbl> <chr>        
## 1 AAN     47.9  44.7  51.1 Manual (RMSE)
## 2 AAN     47.9  44.6  51.2 R (hilo)     
## 3 ANN     47.2  43.8  50.7 Manual (RMSE)
## 4 ANN     47.2  43.7  50.8 R (hilo)

Answer. For ANN my interval is about [43.8, 50.7] and R’s is [43.7, 50.8]. For AAN mine is [44.7, 51.1] and R’s is [44.6, 51.2]. R’s intervals are a little wider in both cases. That’s because RMSE divides the sum of squared errors by \(T\) but ETS() divides by \(T - k\)

Exercise 8.6 — Chinese GDP

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 <- global_economy |>
  filter(Country == "China") |>
  mutate(GDP_b = GDP / 1e9)          # billions of US$ so the axis is readable

china |>
  autoplot(GDP_b) +
  labs(title = "China GDP", y = "Billions of US$", x = "Year")

# Guerrero picks the Box-Cox lambda back from chapter 3 
lambda_china <- china |>
  features(GDP_b, features = guerrero) |>
  pull(lambda_guerrero)
lambda_china
## [1] -0.03446284

GDP grows exponentially, and Guerrero picks \(\lambda \approx -0.03\), which is basically a log.

fit_china <- china |>
  model(
    `ETS auto` = ETS(GDP_b),                     # let AICc choose
    Holt = ETS(GDP_b ~ error("A") + trend("A") + season("N")),
    `Damped Holt` = ETS(GDP_b ~ error("A") + trend("Ad") + season("N")),
    `Box-Cox Holt` = ETS(box_cox(GDP_b, lambda_china) ~
                           error("A") + trend("A") + season("N")),
    `Box-Cox damped Holt` = ETS(box_cox(GDP_b, lambda_china) ~
                                  error("A") + trend("Ad") + season("N"))
  )

fit_china        # the auto model turns out to be ETS(M,A,N)
## # A mable: 1 x 6
## # Key:     Country [1]
##   Country   `ETS auto`         Holt `Damped Holt` `Box-Cox Holt`
##   <fct>        <model>      <model>       <model>        <model>
## 1 China   <ETS(M,A,N)> <ETS(A,A,N)> <ETS(A,Ad,N)>   <ETS(A,A,N)>
## # ℹ 1 more variable: `Box-Cox damped Holt` <model>
# smoothing parameters for each version
tidy(fit_china) |>
  filter(term %in% c("alpha", "beta", "phi")) |>
  select(.model, term, estimate) |>
  pivot_wider(names_from = term, values_from = estimate)
## # A tibble: 5 × 4
##   .model              alpha  beta    phi
##   <chr>               <dbl> <dbl>  <dbl>
## 1 ETS auto            1.000 0.372 NA    
## 2 Holt                1.000 0.552 NA    
## 3 Damped Holt         1.000 0.569  0.980
## 4 Box-Cox Holt        1.000 0.106 NA    
## 5 Box-Cox damped Holt 1.000 0.296  0.959
fc_china <- fit_china |> forecast(h = 20)   # 20 years out, to 2037
# No transformation here just auto, Holt and damped Holt
fc_china |>
  filter(.model %in% c("ETS auto", "Holt", "Damped Holt")) |>
  autoplot(china, level = NULL) +
  labs(title = "China GDP: ETS without a transformation (h = 20)",
       y = "Billions of US$", x = "Year")

# With the Box-Cox transformation, compared with plain Holt
fc_china |>
  filter(.model %in% c("Holt", "Box-Cox Holt", "Box-Cox damped Holt")) |>
  autoplot(china, level = NULL) +
  labs(title = "China GDP: effect of the Box-Cox transformation (h = 20)",
       y = "Billions of US$", x = "Year")

fc_china |>
  as_tibble() |>
  filter(Year %in% c(2018, 2027, 2037)) |>
  select(.model, Year, .mean) |>
  pivot_wider(names_from = Year, values_from = .mean)
## # A tibble: 5 × 4
##   .model              `2018` `2027`  `2037`
##   <chr>                <dbl>  <dbl>   <dbl>
## 1 ETS auto            12982. 19683.  27128.
## 2 Holt                12997. 19835.  27432.
## 3 Damped Holt         12978. 19013.  24549.
## 4 Box-Cox Holt        13891. 48050. 245488.
## 5 Box-Cox damped Holt 13213. 29450.  82048.

Answer. No transformation which is auto extend the trend in a straight line. The auto model pick ETS(M,A,N) and it lands almost on top of plain Holt. A straight line can’t keep up with something that grows exponentially. Damped trend bend the straight line flat with time predictin 24.5 trillion by 2037 instead of 27.4 is a big gap and gap keeps getting bigger as time does.Box-Cox transformation takes on a log-ish scale and has a constant slope meaning constant growth and when it gets transformed back it curves exponetially which would mean you would have to keep up the historical growth.Mixed Box-Cox and damped, the growth rate fades instead of the dollar increasing. Basically its a middle ground. So essentially transformation decides the kind of growth while damping decides how long it lasts.

Exercise 8.7 — Gas production (aus_production)

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 |>
  autoplot(Gas) +
  labs(title = "Australian gas production", y = "Petajoules", x = "Quarter")

# size of the seasonal swing vs the level
gas |>
  as_tibble() |>
  mutate(Year = year(Quarter)) |>
  group_by(Year) |>
  filter(n() == 4) |>                     # full years only
  summarise(seasonal_range = max(Gas) - min(Gas), mean_level = mean(Gas)) |>
  filter(Year %in% c(1960, 1970, 1980, 1990, 2000, 2009))
## # A tibble: 6 × 3
##    Year seasonal_range mean_level
##   <dbl>          <dbl>      <dbl>
## 1  1960              2       7.25
## 2  1970             11      18.2 
## 3  1980             35      86.2 
## 4  1990             38     140.  
## 5  2000             49     187.  
## 6  2009             56     224

multiplicative seasonnality because swings grow with the series, additive assumes that swings are the same size the whole time which is clearly on the case here. Multiplicative sees it as each quarter is a percentage above or below so the swings scale up as production grows.

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

fit_gas |> select(`ETS 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 ETS auto 1682.
## 2 M,A,M    1682.
## 3 M,Ad,M   1685.
## 4 A,A,A    1873.
accuracy(fit_gas) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 4 × 5
##   .model    RMSE   MAE  MAPE  MASE
##   <chr>    <dbl> <dbl> <dbl> <dbl>
## 1 ETS auto  4.60  3.02  4.08 0.542
## 2 M,A,M     4.60  3.02  4.08 0.542
## 3 M,Ad,M    4.59  3.03  4.10 0.544
## 4 A,A,A     4.76  3.35 10.9  0.600

ETS() picked ETS(M,A,M) on its own, which agrees with the plot.

fc_gas <- fit_gas |>
  select(`M,A,M`, `M,Ad,M`) |>
  forecast(h = "3 years")

fc_gas |>
  autoplot(gas |> filter(year(Quarter) >= 2000), level = 80, alpha = 0.5) +
  labs(title = "Gas production: ETS(M,A,M) vs damped ETS(M,Ad,M)",
       y = "Petajoules", x = "Quarter")

Three checks to see if damping helps .

# 1. Hold  the last 3 years and a half  as a test set
gas_train <- gas |> filter(year(Quarter) <= 2006)

gas_train |>
  model(
    `M,A,M` = ETS(Gas ~ error("M") + trend("A") + season("M")),
    `M,Ad,M` = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
  ) |>
  forecast(h = 14) |>
  accuracy(gas) |>
  select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 × 5
##   .model  RMSE   MAE  MAPE  MASE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 M,A,M   8.92  7.86  3.67  1.44
## 2 M,Ad,M  7.27  5.79  2.74  1.06
# 2. cross validation on 2 years ahead
gas |>
  filter(year(Quarter) >= 1990) |>
  stretch_tsibble(.init = 40, .step = 4) |>
  model(
    `M,A,M` = ETS(Gas ~ error("M") + trend("A") + season("M")),
    `M,Ad,M` = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
  ) |>
  forecast(h = 8) |>
  accuracy(gas) |>
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 2 × 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 M,A,M   9.69  7.76  3.86
## 2 M,Ad,M 10.5   8.35  4.14

Answer. Not really. Results go both way AICc prefers the undamped model and the training RMSE stays the same. On the last 3 years and a half the test set damped wins (RMSE 7.27 vs 8.92). In cross-validation over many origins since 1990, undamped wins (RMSE 9.69 vs 10.5). Damping is barely doing anything, over a 3 year range, both forecasts are within about 2 petajoules of each other. I’d just with stick with ETS(M,A,M) which is what ETS() ended up picking . Damping didn’t give a consistent improvement, it just helped on one particular stretch of test data.

Exercise 8.8 — Retail series: Holt-Winters

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

  1. Why is multiplicative seasonality necessary for this series?
  2. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
  3. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
  4. Check that the residuals from the best method look like white noise.
  5. 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(12345678)                          # same seed as my earlier homeworks
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

# Same series as before
myseries |> as_tibble() |> distinct(State, Industry, `Series ID`)
## # A tibble: 1 × 3
##   State              Industry                                        `Series ID`
##   <chr>              <chr>                                           <chr>      
## 1 Northern Territory Clothing, footwear and personal accessory reta… A3349767W

(a) Why multiplicative seasonality?

myseries |>
  autoplot(Turnover) +
  labs(title = "Northern Territory clothing, footwear and accessory turnover",
       y = "$ million AUD", x = "Month")

# seasonal swing vs level, a few years apart
myseries |>
  as_tibble() |>
  mutate(Year = year(Month)) |>
  group_by(Year) |>
  filter(n() == 12) |>
  summarise(seasonal_range = max(Turnover) - min(Turnover),
            mean_level = mean(Turnover)) |>
  filter(Year %in% c(1989, 1995, 2005, 2010, 2015, 2018))
## # A tibble: 6 × 3
##    Year seasonal_range mean_level
##   <dbl>          <dbl>      <dbl>
## 1  1989            2.9       3.26
## 2  1995            7.6       8.37
## 3  2005            5.2       9.12
## 4  2010            8.7      12.4 
## 5  2015           10.5      12.9 
## 6  2018            9.6      14.4

Answer. The seasonal swings get bigger as the turnover gets bigger

(b) Holt-Winters’ multiplicative, with and without damping

fit_retail <- myseries |>
  model(
    `HW multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    `HW damped` = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

tidy(fit_retail) |>
  filter(term %in% c("alpha", "beta", "gamma", "phi")) |>
  select(.model, term, estimate) |>
  pivot_wider(names_from = term, values_from = estimate)
## # A tibble: 2 × 5
##   .model            alpha     beta    gamma    phi
##   <chr>             <dbl>    <dbl>    <dbl>  <dbl>
## 1 HW multiplicative 0.526 0.000103 0.000104 NA    
## 2 HW damped         0.557 0.000238 0.111     0.979
fit_retail |>
  forecast(h = "3 years") |>
  autoplot(myseries |> filter(year(Month) >= 2012), level = NULL) +
  labs(title = "Holt-Winters multiplicative vs damped: 3 year forecasts",
       y = "$ million AUD", x = "Month")

(c) Compare the one-step RMSE

accuracy(fit_retail) |>
  select(.model, ME, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 × 6
##   .model                 ME  RMSE   MAE  MAPE  MASE
##   <chr>               <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HW multiplicative -0.0128 0.613 0.450  5.15 0.513
## 2 HW damped          0.0398 0.616 0.444  5.06 0.507
glance(fit_retail) |> select(.model, AICc)
## # A tibble: 2 × 2
##   .model             AICc
##   <chr>             <dbl>
## 1 HW multiplicative 1776.
## 2 HW damped         1783.

Answer. The one-step RMSEs are basically the same,0.613 for regular Holt-Winters and 0.616 for damped. The damped version is a tiny bit better on MAE and MAPE, and AICc leans toward the regular one.

(d) Residuals from the damped method

fit_retail |>
  select(`HW damped`) |>
  gg_tsresiduals()

lb_test(fit_retail, lag = 24)               # monthly -> lag 2m = 24
## # A tibble: 2 × 5
##   State              Industry                           .model lb_stat lb_pvalue
##   <chr>              <chr>                              <chr>    <dbl>     <dbl>
## 1 Northern Territory Clothing, footwear and personal a… HW da…    23.8    0.472 
## 2 Northern Territory Clothing, footwear and personal a… HW mu…    34.5    0.0765

Answer. Yeah, the damped model’s residuals look like white noise for the most part. The Ljung-Box p-value is about 0.47 leaving no evidence of leftover autocorrelation, and almost all the ACF spikes are inside the blue bounds.

(e) Test set - training through 2010

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

fit_retail_train <- retail_train |>
  model(
    `HW multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    `HW damped` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
    SNAIVE = SNAIVE(Turnover)               # the Chapter 5 benchmark
  )

fc_retail <- fit_retail_train |> forecast(new_data = retail_test)

fc_retail |>
  autoplot(myseries |> filter(year(Month) >= 2007), level = NULL) +
  labs(title = "Forecasts for 2011-2018 ",
       y = "$ million AUD", x = "Month")

accuracy(fc_retail, myseries) |>
  select(.model, ME, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)
## # A tibble: 3 × 6
##   .model                ME  RMSE   MAE  MAPE  MASE
##   <chr>              <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HW damped          0.163  1.15 0.878  6.45 0.960
## 2 SNAIVE             0.836  1.55 1.24   9.06 1.36 
## 3 HW multiplicative -1.54   1.78 1.60  12.6  1.75

Answer. Yes, the damped Holt-Winters model beats seasonal naive. Test RMSE is 1.15 vs 1.55 for SNAIVE while its MAPE is 6.5% vs 9.1%. The regular Holt-Winters model actually did worse than seasonal naive with if RMSE 1.78. Its ME of -1.54 means it over forecasted because it kept projecting the 2000s growth into the later years when turnover actually leveled off for a couple of years.

Exercise 8.9 — STL + ETS on the Box-Cox transformed retail series

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?

decomposition_model() from chapter 5 does the STL step and fits ETS to the seasonally adjusted part then forecasts the seasonal part with seasonal naive and then adds them back together which undoes the Box-Cox.

# lambda from the training data only so the test set isn't looked at beforehand
lambda_retail <- retail_train |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)
lambda_retail
## [1] 0.1971431
fit_stl <- retail_train |>
  model(
    # ETS chooses the model for the seasonally adjusted series by AICc
    `STL + ETS` = decomposition_model(
      STL(box_cox(Turnover, lambda_retail)),
      ETS(season_adjust ~ season("N"))
    ),
    # same thing but forcing a damped trend on the seasonally adjusted series
    `STL + damped ETS` = decomposition_model(
      STL(box_cox(Turnover, lambda_retail)),
      ETS(season_adjust ~ trend("Ad") + season("N"))
    ),
    `HW damped` = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

fit_stl |> select(`STL + ETS`) |> report()    # shows the ETS model it picked
## 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.5934852 
##     beta  = 0.0001000023 
## 
##   Initial states:
##     l[0]        b[0]
##  1.07674 0.008107742
## 
##   sigma^2:  0.008
## 
##      AIC     AICc      BIC 
## 220.6875 220.9123 238.7349 
## 
## Series: season_year 
## Model: SNAIVE 
## 
## sigma^2: 1e-04
fc_stl <- fit_stl |> forecast(new_data = retail_test)

fc_stl |>
  autoplot(myseries |> filter(year(Month) >= 2007), level = NULL) +
  labs(title = "STL + ETS vs damped Holt-Winters, 2011-2018",
       y = "$ million AUD", x = "Month")

accuracy(fc_stl, myseries) |>
  select(.model, ME, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)
## # A tibble: 3 × 6
##   .model               ME  RMSE   MAE  MAPE  MASE
##   <chr>             <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HW damped         0.163  1.15 0.878  6.45 0.960
## 2 STL + damped ETS -0.483  1.16 0.979  7.68 1.07 
## 3 STL + ETS        -3.59   4.01 3.59  27.5  3.93

Answer. The STL and ETS approach doesn’t beat my best model from 8.8. When ETS() picks the model on its own, the test RMSE is 4.01 in comparison to 1.15 for my damped Holt-Winters, but that’s because its trend keeps growing and the forecasts end up way too high. When I add a damped trend, the RMSE drops to 1.16, which basically ties damped Holt-Winters. Since the damped Holt-Winters model is just as accurate and simpler, I’d just stick with it.