Rpubs URL: https://rpubs.com/umaisabdullah/DATA624_HW5

Objective

Complete and submit exercises 8.1, 8.5, 8.6, 8.7, 8.8 and 8.9 for Homework 5 from the Hyndman online Forecasting book. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

Homework 3 fit benchmark methods that carry forward a single number or a straight line. This homework moves to exponential smoothing, which is the first family of methods in the book that actually reacts to how the series has been behaving recently rather than just repeating the past.

The core idea across every method in this chapter is a weighted average of past observations, where the weights decay exponentially the further back in time they go. Simple exponential smoothing (SES) does this for a level only, and is appropriate for series with no trend and no season. Holt’s method adds a second smoothing equation for the trend. Holt-Winters adds a third for the season. Each additional component lets the model track one more feature of the data, at the cost of one more smoothing parameter to estimate.

Three ideas run through the exercises below.

The first is that the smoothing parameters, alpha for the level, beta for the trend, and gamma for the season, are not something I pick by eye. ETS() estimates them, together with the initial states, by maximum likelihood. For models with additive errors that is equivalent to minimizing the sum of squared one-step errors, the same way a regression coefficient is estimated; for multiplicative-error models the likelihood weights the errors relative to the level instead. A value of alpha near 1 means the model trusts the most recent observation almost completely and barely remembers anything older. A value near 0 means the opposite, the model barely moves and mostly repeats its previous level.

The second is damping. An un-damped trend, whether additive or multiplicative, gets extended in a straight line or an exponential curve forever, which becomes unrealistic at long horizons. A damped trend flattens out gradually, which is usually a more honest forecast for anything beyond the very near term.

The third is that ETS is not one model but a state space framework covering additive and multiplicative versions of error, trend and season. Letting ETS() search over that space and pick the best model by AICc is convenient, but it is still worth checking that the chosen model matches what the series actually looks like, rather than trusting the automatic selection blindly.

Pre-Requisite

All of the data for this homework comes from packages, so nothing needs to be downloaded into the working directory. Only one package is required:

fpp3 loads the core forecasting functions and all of the datasets used below.

Exercise 8.1

Question. Consider the number of pigs slaughtered in New South Wales (data set aus_livestock).

Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of alpha and l0, and generate forecasts for the next four months.
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.

My thought process. Part (a) is mechanical: filter to the right series, fit ETS(Count ~ error("A") + trend("N") + season("N")), which is the state space equivalent of SES, and read off the parameters. Part (b) asks me to reproduce the prediction interval by hand, which means pulling the residual standard deviation out of the fitted model and comparing my by-hand number against what hilo() returns. The two should be close but not necessarily identical, and the gap is worth explaining rather than ignoring.

(a) Fit SES via ETS and forecast

nsw_pigs <- aus_livestock |>
  filter(State == "New South Wales", Animal == "Pigs")

nsw_pigs |>
  autoplot(Count) +
  labs(title = "Pigs slaughtered in New South Wales",
       x = "Month", y = "Count")

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

report(fit_pigs)

## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.2655419 
## 
##   Initial states:
##      l[0]
##  107861.1
## 
##   sigma^2:  129400339
## 
##      AIC     AICc      BIC 
## 13955.55 13955.59 13968.52

fc_pigs <- fit_pigs |>
  forecast(h = 4)

fc_pigs |>
  autoplot(nsw_pigs |> filter(year(Month) >= 2010)) +
  labs(title = "NSW pig slaughter, 4 month SES forecast",
       x = "Month", y = "Count")

fc_pigs |>
  as_tibble() |>
  select(Month, .mean)

## # A tibble: 4 × 2
##      Month  .mean
##      <mth>  <dbl>
## 1 2019 Jan 73470.
## 2 2019 Feb 73470.
## 3 2019 Mar 73470.
## 4 2019 Apr 73470.

Reading the fit. The optimal alpha is about 0.266 and the estimated initial level l0 is about 107,861. An alpha of 0.27 is moderate-to-low: the model gives real weight to recent months but still keeps a good amount of memory of the earlier level rather than chasing the most recent value alone. Because this is SES with no trend or season term, every one of the four forecasts is the same flat number, about 73,470, which is the final estimated level.

(b) Prediction interval by hand vs. R

resid_pigs <- residuals(fit_pigs)$.resid
s <- sd(resid_pigs, na.rm = TRUE)

yhat1 <- fc_pigs$.mean[1]

by_hand_lower <- yhat1 - 1.96 * s
by_hand_upper <- yhat1 + 1.96 * s

sigma_pigs <- sqrt(glance(fit_pigs)$sigma2)

c(point_forecast = yhat1, s = s, sigma_model = sigma_pigs,
  by_hand_lower = by_hand_lower, by_hand_upper = by_hand_upper)

## point_forecast              s    sigma_model  by_hand_lower  by_hand_upper 
##       73470.44       11362.84       11375.43       51199.28       95741.60

fc_pigs |>
  hilo(level = 95) |>
  as_tibble() |>
  slice(1) |>
  select(Month, .mean, `95%`)

## # A tibble: 1 × 3
##      Month  .mean                  `95%`
##      <mth>  <dbl>                 <hilo>
## 1 2019 Jan 73470. [51175.01, 95765.87]95

Comparison. The residual standard deviation is about 11,363, so the by-hand 95% interval for the first forecast comes out to roughly 51,199 to 95,742. R’s own interval from hilo() is roughly 51,175 to 95,766, essentially the same bounds, off by only a few tens of units on each side.

The small gap has a precise source. For SES the one-step forecast variance is exactly the innovation variance \(\sigma^2\), so R’s interval is also \(\hat{y} \pm 1.96\hat{\sigma}\). The only difference is how \(\hat{\sigma}\) is estimated: sd() divides the sum of squared residuals by \(T - 1\), while the model’s sigma^2 divides by \(T\) minus the number of estimated parameters (alpha and l0). That gives \(\hat{\sigma} \approx\) 11,375 instead of 11,363, and plugging that into \(\hat{y} \pm 1.96\hat{\sigma}\) reproduces R’s interval. So the 1.96s shortcut is not an approximation R is quietly correcting; it is how the interval is built, with a slightly different degrees-of-freedom adjustment.

Exercise 8.5

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

Plot the Exports series and discuss the main features of the data.
Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
Compute the RMSE values for the training data.
Compare the results to a model with additive trend, ETS(A,A,N). Was this a good idea? Compare with RMSE.
Compare the forecasts from both models. Which do you think is best?
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.

My thought process. I continue with Australia, the same country I have used in earlier homeworks, so the results stay comparable to the plots I have already discussed. The exercise is really asking me to check whether adding a trend term is worth the extra parameter, and RMSE on the training data plus a look at the forecasts is exactly how to check that.

(a) Plot and discuss

aus_exports <- global_economy |>
  filter(Country == "Australia") |>
  select(Year, Exports)

aus_exports |>
  autoplot(Exports) +
  labs(title = "Australian exports as a percentage of GDP",
       x = "Year", y = "Exports (% of GDP)")

Main features. The series is annual, so there is no seasonality to discuss. There is a clear but noisy upward trend across the full 1960 to 2017 span, from roughly 13 percent of GDP to roughly 21 percent, and it is not a clean straight line. Exports fall through the late 1960s, climb through the 1970s, 80s and 90s, then drop sharply in the early 2000s (from about 22 percent around 2001 to about 17 percent a few years later). They recover to a spike of about 23 percent around 2008, pull back during the financial crisis, drift down through the mid-2010s, and tick up again in the final year. There is no obvious structural break or single outlier, just a rising level with multi-year swings around it.

(b) ETS(A,N,N)

fit_ann <- aus_exports |>
  model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))

report(fit_ann)

## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.5659948 
## 
##   Initial states:
##      l[0]
##  12.98943
## 
##   sigma^2:  1.3621
## 
##      AIC     AICc      BIC 
## 257.3943 257.8387 263.5756

fit_ann |>
  forecast(h = 10) |>
  autoplot(aus_exports) +
  labs(title = "Australian exports, ETS(A,N,N) 10 year forecast",
       x = "Year", y = "Exports (% of GDP)")

The fitted alpha is about 0.566, a fairly high value that leans heavily on recent observations. Since this model has no trend term, the forecast is a flat line at the final estimated level, about 20.6 percent of GDP.

(c) RMSE for ETS(A,N,N)

accuracy(fit_ann) |>
  select(.model, RMSE)

## # A tibble: 1 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ANN     1.15

(d) ETS(A,A,N) and RMSE comparison

fit_aan <- aus_exports |>
  model(AAN = ETS(Exports ~ error("A") + trend("A") + season("N")))

report(fit_aan)

## Series: Exports 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.4407571 
##     beta  = 1e-04 
## 
##   Initial states:
##      l[0]      b[0]
##  12.74993 0.1371607
## 
##   sigma^2:  1.3395
## 
##      AIC     AICc      BIC 
## 258.3124 259.4662 268.6146

accuracy(fit_aan) |>
  select(.model, RMSE)

## # A tibble: 1 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 AAN     1.12

bind_rows(
  accuracy(fit_ann) |> select(.model, RMSE),
  accuracy(fit_aan) |> select(.model, RMSE)
)

## # A tibble: 2 × 2
##   .model  RMSE
##   <chr>  <dbl>
## 1 ANN     1.15
## 2 AAN     1.12

Was adding a trend a good idea? Only marginally, and arguably not once the extra parameters are accounted for. The RMSE drops slightly, from about 1.147 for ANN to about 1.117 for AAN. Because ANN is nested inside AAN (set the slope and beta to zero), adding the trend term will normally lower the training RMSE whether or not the trend is actually earning its place, so a small drop is expected. The AIC, which penalizes the extra parameters, says this is not a clear win: ANN’s AIC (about 257.4) is lower, meaning better, than AAN’s (about 258.3).

The fitted beta for AAN is essentially zero (0.0001). That means the slope is never updated; AAN is effectively ANN plus a fixed drift of about 0.137 percentage points per year estimated from the start of the sample. So the small RMSE improvement on the training set is not enough to justify the added complexity.

(e) Comparing the forecasts

aus_exports |>
  model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
        AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))) |>
  forecast(h = 15) |>
  autoplot(aus_exports, level = NULL) +
  labs(title = "Australian exports, ANN vs AAN 15 year forecasts",
       x = "Year", y = "Exports (% of GDP)")

Which is best? ANN forecasts a flat line at roughly 20.6 percent forever. AAN forecasts a gentle continued rise, starting around 20.8 percent and climbing steadily, because it carries forward the fixed drift it estimated from the long-run history and extends it in a straight line indefinitely.

I prefer ANN here. The training RMSE gap is small and the AIC favors ANN once AAN’s extra parameters are penalized, so AAN is not clearly earning its extra complexity. An un-damped linear trend extended 15 years into the future is also exactly the kind of forecast that looks overconfident once you get past the first few years. Exports as a share of GDP cannot rise forever, and the series itself has already shown it can reverse direction for several years at a time, as it did in the late 1960s and the early 2000s. A flat forecast is the more conservative and, I think, the more defensible choice at this horizon. If I did want to keep a trend, a damped version, ETS(A,Ad,N), would be the better compromise, since it would capture the upward drift without committing to it running in a straight line for 15 years.

(f) Prediction intervals by hand vs. R

rmse_ann <- accuracy(fit_ann)$RMSE
rmse_aan <- accuracy(fit_aan)$RMSE

fc_ann1 <- fit_ann |> forecast(h = 1)
fc_aan1 <- fit_aan |> forecast(h = 1)

by_hand <- tibble(
  model = c("ANN", "AAN"),
  point_forecast = c(fc_ann1$.mean, fc_aan1$.mean),
  rmse = c(rmse_ann, rmse_aan)
) |>
  mutate(lower = point_forecast - 1.96 * rmse,
         upper = point_forecast + 1.96 * rmse)

by_hand

## # A tibble: 2 × 5
##   model point_forecast  rmse lower upper
##   <chr>          <dbl> <dbl> <dbl> <dbl>
## 1 ANN             20.6  1.15  18.4  22.9
## 2 AAN             20.8  1.12  18.6  23.0

bind_rows(
  fc_ann1 |> hilo(level = 95) |> as_tibble() |> mutate(model = "ANN"),
  fc_aan1 |> hilo(level = 95) |> as_tibble() |> mutate(model = "AAN")
) |>
  select(model, Year, .mean, `95%`)

## # A tibble: 2 × 4
##   model  Year .mean                  `95%`
##   <chr> <dbl> <dbl>                 <hilo>
## 1 ANN    2018  20.6 [18.31970, 22.89462]95
## 2 AAN    2018  20.8 [18.57028, 23.10700]95

Comparison. For ANN, the by-hand interval is roughly 18.36 to 22.85, and R’s interval is roughly 18.32 to 22.89. For AAN the by-hand interval is roughly 18.65 to 23.03, against R’s roughly 18.57 to 23.11. Both pairs are close but not identical, and R’s are slightly wider in both cases. The gap comes from the same source as in Exercise 8.1. RMSE divides the sum of squared errors by the number of observations, while the model’s sigma^2 divides by the number of observations minus the number of estimated parameters, so it is a little larger. The gap is bigger for AAN because it estimates more parameters. The agreement is close enough to confirm the 1.96 times RMSE shortcut is a reasonable approximation, with R applying a degrees-of-freedom correction on top.

Exercise 8.6

Question. Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the trend and damped parameters, and take a Box-Cox transformation of the data first. Model the logarithmic growth rate. Compare your forecasts with the actual growth rate.

My thought process. China’s GDP is the textbook case for why a Box-Cox transformation matters: it grows so fast, and the growth compounds so strongly, that additive-trend methods on the raw series cannot follow its shape. The exercise wants me to try several ETS specifications side by side and see which handles that acceleration most sensibly. I fit an un-damped trend, a damped trend, and versions on both the raw series and a Box-Cox transformed series so the effect of the transformation is visible on its own.

china_gdp <- global_economy |>
  filter(Country == "China") |>
  select(Year, GDP)

china_gdp |>
  autoplot(GDP) +
  labs(title = "Chinese GDP, current US dollars",
       x = "Year", y = "GDP")

The series is not just trending up, it is accelerating: growth from 2000 onward dwarfs anything before it. That curvature is exactly what a plain additive trend struggles with.

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

lambda_china

## [1] -0.03446284

fit_china <- china_gdp |>
  model(
    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, lambda_china) ~
                           error("A") + trend("A") + season("N")),
    ETS_BoxCox_AAdN = ETS(box_cox(GDP, lambda_china) ~
                            error("A") + trend("Ad") + season("N"))
  )

glance(fit_china) |>
  select(.model, AIC, AICc, BIC)

## # A tibble: 4 × 4
##   .model            AIC  AICc   BIC
##   <chr>           <dbl> <dbl> <dbl>
## 1 ETS_AAN         3258. 3259. 3268.
## 2 ETS_AAdN        3260. 3262. 3273.
## 3 ETS_BoxCox_AAN  -139. -137. -128.
## 4 ETS_BoxCox_AAdN -135. -133. -123.

The information criteria can only be compared within a scale. The two raw-GDP models are fitted to dollars and the two Box-Cox models to transformed values, so the gap between roughly 3,260 and roughly -137 says nothing about which pair is better. Within each pair, the un-damped model has the slightly lower AICc.

fc_china <- fit_china |>
  forecast(h = 20)

fc_china |>
  autoplot(china_gdp, level = NULL) +
  labs(title = "Chinese GDP, 20 year ETS forecasts",
       x = "Year", y = "GDP")

last_gdp <- last(china_gdp$GDP)

implied_growth <- fc_china |>
  as_tibble() |>
  group_by(.model) |>
  summarise(gdp_2037 = last(.mean)) |>
  mutate(multiple_of_2017 = gdp_2037 / last_gdp,
         avg_log_growth = (log(gdp_2037) - log(last_gdp)) / 20)

implied_growth

## # A tibble: 4 × 4
##   .model          gdp_2037 multiple_of_2017 avg_log_growth
##   <chr>              <dbl>            <dbl>          <dbl>
## 1 ETS_AAN          2.74e13             2.24         0.0404
## 2 ETS_AAdN         2.45e13             2.00         0.0347
## 3 ETS_BoxCox_AAN   2.24e14            18.3          0.145 
## 4 ETS_BoxCox_AAdN  5.11e13             4.18         0.0715

g_boxcox_aan  <- implied_growth |> filter(.model == "ETS_BoxCox_AAN")  |> pull(avg_log_growth)
g_boxcox_aadn <- implied_growth |> filter(.model == "ETS_BoxCox_AAdN") |> pull(avg_log_growth)
g_raw_aan     <- implied_growth |> filter(.model == "ETS_AAN")         |> pull(avg_log_growth)
m_boxcox_aan  <- implied_growth |> filter(.model == "ETS_BoxCox_AAN")  |> pull(multiple_of_2017)

Effect of damping and the Box-Cox transformation. The four forecasts fan out in a way that is the opposite of what I first expected.

The two raw-scale models are the most conservative. A linear trend on the dollar scale adds a fixed dollar amount every year, so from 2018 onward both lines grow in a straight line and sit almost on top of each other. The damping parameter barely changes anything over 20 years. Because GDP compounds, a straight line in dollars means the implied growth rate keeps shrinking. The raw AAN path averages only about 4 percent log growth a year, well below both the full-sample average and the average since 2008 (shown below). So these models under-forecast rather than over-forecast.

The Box-Cox transformation changes the picture completely. With a lambda of -0.034, the transformation is very close to a log, and an additive trend on the log scale is a constant percentage growth rate on the original scale. That matches how GDP behaves, because economic growth compounds. But the un-damped Box-Cox AAN model takes the rapid growth of the last few decades and compounds it forward for 20 years without limit. It averages about 14.5 percent log growth a year and ends at roughly 18 times the 2017 level. That is not a believable forecast for an economy that is already one of the largest in the world.

The damped Box-Cox model sits in between. It keeps the compounding shape but lets the growth rate fade, averaging about 7.1 percent a year. Of the four, this is the forecast I find most defensible: the transformation gets the shape right and damping stops the growth rate from being projected forward unchanged forever.

china_growth <- china_gdp |>
  mutate(log_gdp = log(GDP),
         growth_rate = difference(log_gdp)) |>
  filter(!is.na(growth_rate))

china_growth |>
  autoplot(growth_rate) +
  labs(title = "Chinese GDP, log growth rate (first difference of log GDP)",
       x = "Year", y = "Log growth rate")

fit_growth <- china_growth |>
  model(ETS_growth = ETS(growth_rate ~ error("A") + trend("N") + season("N")))

report(fit_growth)

## Series: growth_rate 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.1088357 
## 
##   Initial states:
##        l[0]
##  0.04310116
## 
##   sigma^2:  0.0087
## 
##       AIC      AICc       BIC 
## -36.33972 -35.88689 -30.21057

fc_growth <- fit_growth |>
  forecast(h = 10)

fc_growth |>
  autoplot(china_growth) +
  labs(title = "Forecast of Chinese log GDP growth rate",
       x = "Year", y = "Log growth rate")

g_forecast <- fc_growth$.mean[1]
g_hist     <- mean(china_growth$growth_rate)
g_recent   <- china_growth |> filter(Year >= 2008) |> pull(growth_rate) |> mean()

c(growth_forecast = g_forecast, historical_mean = g_hist, mean_since_2008 = g_recent)

## growth_forecast historical_mean mean_since_2008 
##      0.11552989      0.09338016      0.12369592

Comparing forecasts with the actual growth rate. The historical log growth rate is noisy and trendless, swinging between roughly -0.18 and 0.26 with no clear direction across the sample. It includes sharp drops around 1961 and 1978 and a softening in the mid-2010s. An SES model on the growth rate forecasts a flat line at about 0.116, against a full-sample average of about 0.093, rather than extrapolating a trend that is not really there in the growth rate itself. The fitted alpha is small, about 0.109, so the model updates its estimate of the average growth rate only slowly in response to new data. That is appropriate for a quantity that fluctuates around a roughly stable mean rather than trending.

Lining this up against the level models makes the comparison concrete:

The growth-rate model expects about 11.6 percent a year.
The damped Box-Cox model implies about 7.1 percent.
The un-damped Box-Cox model implies about 14.5 percent.
The raw-scale linear model implies only about 4 percent.

The actual growth rate since 2008 has averaged about 12.4 percent and has been slowing. That is another reason to distrust the un-damped Box-Cox forecast: it projects forward a growth rate the recent data no longer supports.

That is the real value of modelling the growth rate instead of the level. On the level scale, GDP looks like it has an ever-steepening trend, which tempts a model into extrapolating something explosive or, on the raw scale, into flattening it out. On the log-difference scale, the trend disappears and what is left looks close to a stable process fluctuating around an average growth rate, which is a far more defensible basis for a long-run forecast. One caveat: this series is nominal GDP in current US dollars, so its growth rate mixes real growth with inflation and exchange-rate movements.

Exercise 8.7

Question. Find an ETS model for the Gas data from aus_production. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?

My thought process. Gas production has a seasonal swing that grows as the overall level of the series grows, which is exactly the signature that calls for multiplicative rather than additive seasonality. I check that visually first, then let ETS() confirm it, then compare damped and un-damped trend versions on accuracy and on how sensible the long-run forecast looks.

gas <- aus_production |>
  select(Quarter, Gas)

gas |>
  autoplot(Gas) +
  labs(title = "Australian quarterly gas production",
       x = "Quarter", y = "Gas production (petajoules)")

Why multiplicative seasonality. The seasonal swings are small and roughly constant in the early part of the series, then grow visibly larger as production rises through the 1970s onward. When the size of the seasonal fluctuation scales with the level of the series rather than staying a fixed absolute amount, additive seasonality cannot represent it, because additive seasonal indices are constant amounts added on top of the trend regardless of how high the series has climbed. Multiplicative seasonality instead expresses the season as a percentage of the current level, so the swing automatically grows as the series grows. That is the pattern visible in this plot.

fit_gas <- gas |>
  model(
    ETS_auto  = ETS(Gas),
    ETS_MAM   = ETS(Gas ~ error("M") + trend("A") + season("M")),
    ETS_MAdM  = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
  )

glance(fit_gas) |>
  select(.model, AIC, AICc, BIC)

## # A tibble: 3 × 4
##   .model     AIC  AICc   BIC
##   <chr>    <dbl> <dbl> <dbl>
## 1 ETS_auto 1681. 1682. 1711.
## 2 ETS_MAM  1681. 1682. 1711.
## 3 ETS_MAdM 1684. 1685. 1718.

accuracy(fit_gas) |>
  select(.model, RMSE, MAPE, MASE)

## # A tibble: 3 × 4
##   .model    RMSE  MAPE  MASE
##   <chr>    <dbl> <dbl> <dbl>
## 1 ETS_auto  4.60  4.08 0.542
## 2 ETS_MAM   4.60  4.08 0.542
## 3 ETS_MAdM  4.59  4.10 0.544

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

ETS() run with no restrictions selects ETS(M,A,M) on its own, confirming that multiplicative error and multiplicative season are what the automatic search prefers, matching the visual read above.

fit_gas |>
  select(ETS_MAdM) |>
  report()

## 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

phi_gas <- tidy(fit_gas) |>
  filter(.model == "ETS_MAdM", term == "phi") |>
  pull(estimate)

fit_gas |>
  forecast(h = 12) |>
  autoplot(gas |> filter(year(Quarter) >= 2000), level = NULL) +
  labs(title = "Gas production, damped vs un-damped ETS forecasts",
       x = "Quarter", y = "Gas production (petajoules)")

Does damping improve the forecasts? Not on this series. In-sample accuracy is essentially tied: RMSE, MAPE and MASE differ only in the second or third decimal place. The information criteria actually lean the other way. MAM has the lower AIC, AICc and BIC by a few points, so once the extra damping parameter is penalized, the un-damped model is preferred.

The estimated damping parameter explains why the two forecasts are so similar: phi is about 0.98. That is consistent with the plot of the series, which does not show the trend slowing down through the sample, so there is little in the data pushing phi away from 1. When phi is that close to 1, the damped model behaves almost exactly like its un-damped counterpart over a short horizon, and the two forecast lines in the plot nearly overlap.

Damping would matter more at a longer horizon, where an un-damped trend keeps extending while a damped one gradually flattens. For this data and a 12-quarter horizon, though, damping does not improve the forecasts, and the AICc favors the un-damped MAM model that ETS() chose automatically. I would only reach for the damped version as a deliberate hedge if I needed forecasts many years out.

Exercise 8.8

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

Why is multiplicative seasonality necessary for this series?
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
Check that the residuals from the best method look like white noise.
Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naive approach from Exercise 5.7?

My thought process. This is the same NSW takeaway food series from Homeworks 1, 2 and 3, drawn with the same seed so the results stay comparable. The exercise walks through the same damped-vs-not comparison as 8.7 but now on a series I already know has both a strong trend and a growing seasonal swing, and it adds a genuine test set comparison against the SNAIVE benchmark from Homework 3, which is the sharpest check of all: does exponential smoothing actually forecast better, not just fit better.

(a) Why multiplicative seasonality

set.seed(624)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`, 1))

myseries |>
  as_tibble() |>
  distinct(State, Industry, `Series ID`)

## # A tibble: 1 × 3
##   State           Industry               `Series ID`
##   <chr>           <chr>                  <chr>      
## 1 New South Wales Takeaway food services A3349792X

myseries |>
  autoplot(Turnover) +
  labs(title = "NSW takeaway food services turnover",
       x = "Month", y = "Turnover ($ million)")

Same reasoning as the gas series in Exercise 8.7. Turnover has both a strong upward trend and a December peak that grows in absolute size as the overall level of the series rises: the seasonal swing in the 1980s is a small fraction of the swing visible by the 2010s. Since the seasonal amplitude scales with the level rather than staying fixed, only multiplicative seasonality can represent it correctly.

(b) Holt-Winters multiplicative, damped and un-damped

fit_hw <- myseries |>
  model(
    HW_MAM  = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_MAdM = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

fit_hw |>
  select(HW_MAM) |>
  report()

## Series: Turnover 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.7150896 
##     beta  = 0.03088603 
##     gamma = 0.1055336 
## 
##   Initial states:
##      l[0]       b[0]     s[0]     s[-1]    s[-2]    s[-3]     s[-4]     s[-5]
##  82.38423 -0.8137261 1.006568 0.9563809 1.061494 1.080698 0.9912616 0.9966824
##     s[-6]     s[-7]     s[-8]     s[-9]   s[-10]    s[-11]
##  0.961185 0.9961148 0.9945169 0.9640538 1.001246 0.9897988
## 
##   sigma^2:  0.0017
## 
##      AIC     AICc      BIC 
## 4634.521 4635.968 4704.035

fit_hw |>
  select(HW_MAdM) |>
  report()

## Series: Turnover 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.7751114 
##     beta  = 0.007524738 
##     gamma = 0.09682012 
##     phi   = 0.9799997 
## 
##   Initial states:
##      l[0]      b[0]     s[0]     s[-1]    s[-2]    s[-3]     s[-4]    s[-5]
##  83.44206 0.2885412 1.011822 0.9520354 1.055695 1.079874 0.9936296 1.000433
##      s[-6]     s[-7]     s[-8]     s[-9]    s[-10]    s[-11]
##  0.9745973 0.9915392 0.9945614 0.9683761 0.9932686 0.9841682
## 
##   sigma^2:  0.0017
## 
##      AIC     AICc      BIC 
## 4625.238 4626.859 4698.841

fit_hw |>
  forecast(h = "3 years") |>
  autoplot(myseries |> filter(year(Month) >= 2005), level = NULL) +
  labs(title = "NSW takeaway turnover, Holt-Winters damped vs un-damped",
       x = "Month", y = "Turnover ($ million)")

(c) Compare one-step RMSE

accuracy(fit_hw) |>
  select(.model, RMSE, MAPE, MASE)

## # A tibble: 2 × 4
##   .model   RMSE  MAPE  MASE
##   <chr>   <dbl> <dbl> <dbl>
## 1 HW_MAM   9.74  3.08 0.308
## 2 HW_MAdM  9.73  3.03 0.307

Which do I prefer? The two are essentially tied on the training data: RMSE of about 9.74 for HW_MAM against 9.73 for HW_MAdM, with MAPE and MASE equally close. Unlike the gas series, the information criteria here favor the damped model: HW_MAdM’s AIC is about 4625 against about 4634 for HW_MAM, so the damping parameter earns its place even after the penalty. I prefer HW_MAdM. It fits at least as well, the AIC supports it, and it protects against overconfident straight-line extrapolation at longer horizons.

(d) Residual check on the preferred model

fit_hw |>
  select(HW_MAdM) |>
  gg_tsresiduals()

augment(fit_hw) |>
  filter(.model == "HW_MAdM") |>
  features(.innov, ljung_box, lag = 24)

## # A tibble: 1 × 5
##   State           Industry               .model  lb_stat lb_pvalue
##   <chr>           <chr>                  <chr>     <dbl>     <dbl>
## 1 New South Wales Takeaway food services HW_MAdM    25.3     0.392

Do the residuals look like white noise? Yes. The residual time plot shows roughly constant spread with no strong pattern, and the ACF sits almost entirely inside the significance bounds, with only an isolated spike or two of the kind expected by chance across 24 lags. The Ljung-Box test at lag 24 gives a statistic of about 25.3 with a p-value of about 0.39, well above 0.05. That means there is no significant evidence against the residuals being uncorrelated. This is a real contrast with the seasonal naive residuals from Homework 3, which showed strong, undecaying autocorrelation at almost every lag and a Ljung-Box p-value indistinguishable from zero. Holt-Winters, by explicitly modelling level, trend and multiplicative season together, has captured essentially all of the structure that seasonal naive left behind in its residuals.

(e) Test set RMSE vs. seasonal naive

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

fit_compare <- myseries_train |>
  model(
    SNAIVE  = SNAIVE(Turnover),
    HW_MAM  = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    HW_MAdM = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
  )

fc_compare <- fit_compare |>
  forecast(new_data = myseries_test)

fc_compare |>
  autoplot(myseries, level = NULL) +
  labs(title = "NSW takeaway turnover, test period: SNAIVE vs Holt-Winters",
       x = "Month", y = "Turnover ($ million)")

accuracy(fc_compare, myseries) |>
  select(.model, RMSE, MAE, MAPE, MASE)

## # A tibble: 3 × 5
##   .model   RMSE   MAE  MAPE  MASE
##   <chr>   <dbl> <dbl> <dbl> <dbl>
## 1 HW_MAM   70.4  60.9  15.3  3.17
## 2 HW_MAdM  69.0  62.8  14.7  3.27
## 3 SNAIVE   96.8  79.5  16.3  4.14

Can Holt-Winters beat seasonal naive? Yes, though the margin on the test set is more modest than the training fit alone would suggest. In Homework 3, SNAIVE trained through 2010 gave a test RMSE of 96.8 on this exact split. HW_MAM comes in at 70.4 and HW_MAdM at 69.0, roughly 28 to 29 percent lower RMSE. MASE is around 3.2 for both against SNAIVE’s 4.14, and MAE is about 61 to 63 against 79.5. MAPE improves too, but only modestly (14.7 and 15.3 against 16.3). The two Holt-Winters models trade places depending on the measure: the un-damped model has the lower MAE (60.9 vs 62.8), while the damped model has the lower RMSE and MAPE (69.0 and 14.7 vs 70.4 and 15.3). Both beat the benchmark on every measure, clearly on RMSE, MAE and MASE.

This is the payoff for the extra complexity in this chapter, even if it is smaller than the near-perfect training residuals might suggest. SNAIVE only ever repeats last year’s same month, so on a series with a strong ongoing trend it falls further and further behind as the forecast horizon grows, which is exactly the failure documented in Homework 3. Holt-Winters explicitly models the trend, so it keeps climbing along with the actual series instead of forecasting a flat repeat of 2010, and it models the multiplicative season on top of that climbing level, keeping the December peaks in roughly the right proportion. The gap between the excellent training fit (RMSE under 10) and the much larger test RMSE (around 70) is also a useful reminder from Homework 3’s lesson: training accuracy is not test accuracy, and an eight year-ahead forecast is a genuinely harder problem than the one-step-ahead errors the model was fitted on, even for a model this much better specified than a naive benchmark.

Exercise 8.9

Question. 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 this compare with your best previous forecasts on the test set?

My thought process. This combines two ideas from earlier in the book: the Box-Cox transformation from Homework 2 and STL decomposition, together with the ETS approach just used in Exercise 8.8. decomposition_model() handles the mechanics: it STL-decomposes the Box-Cox transformed series, fits an ETS model to only the seasonally adjusted component, and forecasts the seasonal component separately (with a seasonal naive method) before adding it back in. I fit it on the same pre-2011 training window so the test set RMSE lines up directly against the SNAIVE and Holt-Winters numbers from part (e). I first let ETS() choose the error and trend for the seasonally adjusted series, then try a damped trend version for comparison.

lambda_retail <- myseries_train |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

lambda_retail

## [1] -0.1011142

stl_spec <- STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic"),
                robust = TRUE)

fit_stl_ets <- myseries_train |>
  model(
    STL_ETS = decomposition_model(
      stl_spec,
      ETS(season_adjust ~ season("N"))
    ),
    STL_ETS_damped = decomposition_model(
      stl_spec,
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
    )
  )

fit_stl_ets |>
  select(STL_ETS) |>
  report()

## Series: Turnover 
## Model: STL decomposition model 
## Transformation: box_cox(Turnover, lambda_retail) 
## Combination: ETS(A,A,N) + SNAIVE
## 
## ================================
## 
## Series: season_adjust 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.7369909 
##     beta  = 0.0001000001 
## 
##   Initial states:
##      l[0]        b[0]
##  3.587272 0.002469821
## 
##   sigma^2:  6e-04
## 
##       AIC      AICc       BIC 
## -533.0937 -532.9167 -513.8760 
## 
## Series: season_year 
## Model: SNAIVE 
## 
## sigma^2: 0

stl_report <- capture.output(fit_stl_ets |> select(STL_ETS) |> report())
stl_has_trend  <- any(grepl("beta", stl_report))
stl_has_damped <- any(grepl("phi", stl_report))

fc_stl_ets <- fit_stl_ets |>
  forecast(new_data = myseries_test)

fc_stl_ets |>
  autoplot(myseries |> filter(year(Month) >= 2000), level = NULL) +
  labs(title = "NSW takeaway turnover, STL + Box-Cox + ETS on the test period",
       x = "Month", y = "Turnover ($ million)")

acc_all <- bind_rows(
  accuracy(fc_compare, myseries),
  accuracy(fc_stl_ets, myseries)
) |>
  select(.model, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)

acc_all

## # A tibble: 5 × 5
##   .model          RMSE   MAE  MAPE  MASE
##   <chr>          <dbl> <dbl> <dbl> <dbl>
## 1 HW_MAdM         69.0  62.8  14.7  3.27
## 2 HW_MAM          70.4  60.9  15.3  3.17
## 3 STL_ETS_damped  76.5  67.3  14.6  3.50
## 4 STL_ETS         93.9  90.2  20.9  4.69
## 5 SNAIVE          96.8  79.5  16.3  4.14

get_acc <- function(m, col) acc_all |> filter(.model == m) |> pull({{ col }})

stl_rmse <- get_acc("STL_ETS", RMSE)
stl_mae  <- get_acc("STL_ETS", MAE)
stl_mape <- get_acc("STL_ETS", MAPE)
stl_mase <- get_acc("STL_ETS", MASE)
dmp_rmse <- get_acc("STL_ETS_damped", RMSE)
dmp_mase <- get_acc("STL_ETS_damped", MASE)
hw_rmse  <- get_acc("HW_MAdM", RMSE)
sn_rmse  <- get_acc("SNAIVE", RMSE)
sn_mase  <- get_acc("SNAIVE", MASE)

How does it compare? Not well, and this is the more interesting and slightly humbling result of the two exercises. The STL plus Box-Cox plus ETS model gives a test RMSE of about 93.9, an MAE of about 90.2, a MAPE of about 20.9 and a MASE of about 4.69. That is only marginally better than SNAIVE on RMSE (96.8) and clearly worse than either Holt-Winters model from Exercise 8.8. On MAE, MAPE and MASE it is actually worse than SNAIVE. Because MASE scales every model’s errors by the same in-sample seasonal naive MAE, it always ranks models the same way MAE does, so the MASE and MAE results tell the same story: on a typical month this model misses by more than seasonal naive does. It only edges ahead on RMSE because SNAIVE’s errors include a few larger misses that RMSE punishes more heavily.

This is not the outcome I expected going in, and it is worth being honest about why rather than explaining it away. The plot shows the problem directly: the forecast does not lag behind the series, it overshoots it. The STL forecast climbs well above the actual turnover through the whole test period and ends far above where the series actually finished in 2018.

The report() output above shows the ETS model chosen for the seasonally adjusted series, which includes an un-damped trend. The Box-Cox lambda of -0.101 is close to zero, so the transformed series behaves almost like a log. A linear trend on that scale becomes compounding percentage growth back on the original scale. The model takes the fast growth of the 2000s, including the jump in 2010 right at the end of the training window, and compounds it forward for eight years. Meanwhile the actual series dipped from 2011 to about 2014 before recovering, so the forecast starts too high and the gap widens with the horizon. A second, smaller factor is that season(window = "periodic") fixes one seasonal shape for the whole horizon and cannot let it drift.

The damped version tests whether reining in the trend helps. It gives a test RMSE of about 76.5 and a MASE of about 3.5, an improvement over the automatically selected model. It still does not beat the best Holt-Winters model on RMSE. Unlike the original STL model, it beats SNAIVE on every measure, and its MAPE is essentially level with the Holt-Winters models. The best forecast on this test set by RMSE remains HW_MAdM, with an RMSE of 69. The jump from the un-damped to the damped STL model is the clearest evidence in this homework that the problem was the over-extrapolated trend, not the decomposition step itself.

Overall verdict across the whole homework. The two approaches in 8.8 and 8.9 are not interchangeable, and the honest conclusion is that neither exponential smoothing method is automatically better just because it is more sophisticated than a benchmark. Holt-Winters, which estimates level, trend and multiplicative season jointly in one set of smoothing equations, clearly beat the seasonal naive floor from Homework 3, cutting test RMSE by close to 30 percent. The STL-plus-ETS pipeline with an automatically chosen ETS model barely beat the benchmark on RMSE and lost to it on MAE, MAPE and MASE, because the growth it extrapolated on a near-log scale ran well ahead of what actually happened. Damping that trend recovered most of the gap and beat the benchmark on every measure, though it still trailed Holt-Winters on RMSE. The lesson is the same one that ran through Homework 3: a model only outperforms its benchmark to the extent that it captures the real structure in the series without over-reading it, and simply layering more machinery into the forecasting pipeline does not guarantee that happens. Getting the specification right, here especially how the trend is allowed to behave over a long horizon, matters more than which general family of methods is used.

Homework 5: Chapter 8 Exercises From FPP3

Umais Siddiqui

2026-10-10