This assignment covers Exercises 8.1, 8.5, 8.6, 8.7, 8.8, and 8.9 from Forecasting: Principles and Practice (3rd edition).

Exercise 8.1: Pigs Slaughtered in Victoria

Exercise 8.1(a): Fit Simple Exponential Smoothing

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

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

report(pigs_fit)
## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.3221247 
## 
##   Initial states:
##      l[0]
##  100646.6
## 
##   sigma^2:  87480760
## 
##      AIC     AICc      BIC 
## 13737.10 13737.14 13750.07

Simple exponential smoothing is represented by an ETS(A,N,N) model: additive errors, no trend, and no seasonality.

Exercise 8.1(a): Four-Month Forecasts

The estimated smoothing parameter is \(\alpha = 0.3221247\), and the estimated initial level is \(\ell_0 = 100646.6\) pigs.

pigs_fc <- pigs_fit |>
  forecast(h = 4)

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

Exercise 8.1(a): Forecast Plot

pigs_fc |>
  autoplot(victoria_pigs) +
  geom_line(
    aes(y = .fitted),
    colour = "#D55E00",
    data = augment(pigs_fit)
  ) +
  labs(
    title = "Victoria Pig Slaughter: Simple Exponential Smoothing",
    x = "Month",
    y = "Number of pigs slaughtered"
  ) +
  guides(colour = "none")

The point forecast is approximately 95,187 pigs for each month from January through April 2019. Simple exponential smoothing produces flat point forecasts equal to the final estimated level. The prediction intervals widen with the forecast horizon.

Exercise 8.1(b): First Forecast’s 95% Prediction Interval

The manual interval is calculated as \(\hat{y} \pm 1.96s\), where \(s\) is the sample standard deviation of the residuals.

pigs_residuals <- augment(pigs_fit)

pigs_s <- sd(pigs_residuals$.resid, na.rm = TRUE)
pigs_first_mean <- pigs_fc$.mean[1]

pigs_manual_interval <- data.frame(
  Forecast = pigs_first_mean,
  ResidualSD = pigs_s,
  Lower95 = pigs_first_mean - 1.96 * pigs_s,
  Upper95 = pigs_first_mean + 1.96 * pigs_s
)

knitr::kable(pigs_manual_interval, digits = 2)
Forecast ResidualSD Lower95 Upper95
95186.56 9344.67 76871.01 113502.1 
pigs_fc |>
  slice_head(n = 1) |>
  hilo(level = 95)
## # A tsibble: 1 x 7 [1M]
## # Key:       Animal, State, .model [1]
##   Animal State    .model    Month
##   <fct>  <fct>    <chr>     <mth>
## 1 Pigs   Victoria SES    2019 Jan
## # ℹ 3 more variables: Count <dist>, .mean <dbl>, `95%` <hilo>

The manual 95% prediction interval for January 2019 is [76,871.01, 113,502.10] pigs, calculated using a forecast of 95,186.56 and a residual standard deviation of 9,344.67.

R produces [76,854.79, 113,518.30], which is slightly wider. The small difference reflects the model’s variance estimate versus the sample standard deviation used in the manual calculation, as well as using 1.96 instead of the exact normal quantile. Both calculations give essentially the same interval.

Exercise 8.5: Exports from Algeria

Exercise 8.5(a): Plot the Exports Series

algeria_economy <- global_economy |>
  filter(Country == "Algeria")

algeria_economy |>
  autoplot(Exports) +
  labs(
    y = "% of GDP",
    title = "Exports: Algeria"
  )

Algeria’s exports as a percentage of GDP fluctuate substantially from 1960 to 2017, with no consistent upward or downward trend across the entire period. The series shows sharp changes, including a spike in the early 1960s, a low in the late 1980s, and relatively high values in the mid-2000s.

Exports decline markedly toward the end of the series, followed by a small recovery in 2017. Because the observations are annual, within-year seasonality cannot be assessed. The changing level suggests simple exponential smoothing is a useful baseline, although it will not extrapolate the recent decline.

Exercise 8.5(b): Simple Exponential Smoothing Forecasts

algeria_ses_fit <- algeria_economy |>
  model(
    SES = ETS(
      Exports ~ error("A") + trend("N") + season("N")
    )
  )

algeria_ses_fc <- algeria_ses_fit |>
  forecast(h = 5)

algeria_ses_fc |>
  autoplot(algeria_economy) +
  geom_line(
    aes(y = .fitted),
    colour = "#D55E00",
    data = augment(algeria_ses_fit)
  ) +
  labs(
    title = "Algerian Exports: ETS(A,N,N) Forecasts",
    x = "Year",
    y = "Exports (% of GDP)"
  ) +
  guides(colour = "none")

Exercise 8.5(c): Training RMSE

algeria_ses_accuracy <- algeria_ses_fit |>
  accuracy()

algeria_ses_accuracy |>
  select(.model, .type, RMSE) |>
  knitr::kable(digits = 4)
.model .type RMSE
SES Training 5.8653

The ETS(A,N,N) model has a training RMSE of 5.8653 percentage points of GDP. This measures the size of the model’s one-step-ahead errors on the training data. It does not measure accuracy on future, unseen observations.

Exercise 8.5(d): Compare ETS(A,N,N) and ETS(A,A,N)

algeria_models <- algeria_economy |>
  model(
    SES = ETS(
      Exports ~ error("A") + trend("N") + season("N")
    ),
    Holt = ETS(
      Exports ~ error("A") + trend("A") + season("N")
    )
  )

algeria_models |>
  accuracy() |>
  select(.model, .type, RMSE) |>
  knitr::kable(digits = 4)
.model .type RMSE
SES Training 5.8653
Holt Training 5.8886

The ETS(A,N,N) model has a training RMSE of 5.8653, compared with 5.8886 for ETS(A,A,N). The difference is small, but adding a trend does not improve the training RMSE here.

SES adapts to changes in level and produces flat forecasts. Holt’s method estimates both a level and a trend, allowing forecasts to rise or fall. This flexibility introduces an additional trend smoothing parameter and an initial trend state.

Because Algeria’s exports do not show a consistent long-term trend, the simpler SES model is a reasonable choice. However, training RMSE alone does not establish which model will forecast unseen data more accurately.

Exercise 8.5(e): Compare Forecasts

algeria_fc <- algeria_models |>
  forecast(h = 5)

algeria_fc |>
  autoplot(algeria_economy, level = NULL) +
  labs(
    title = "Algerian Exports: SES and Holt Forecasts",
    x = "Year",
    y = "Exports (% of GDP)"
  ) +
  guides(colour = guide_legend(title = "Model"))

SES forecasts exports at a constant level of approximately 22.4% of GDP. Holt’s method produces gradually declining forecasts by extending its estimated downward trend.

I prefer SES for this series because exports have fluctuated considerably without a persistent long-term trend. The small recovery in the final observation also makes continued decline uncertain. SES is simpler and has a slightly lower training RMSE (5.8653 versus 5.8886). However, this preference should be checked using out-of-sample forecast accuracy.

Exercise 8.5(f): First Forecast’s 95% Prediction Intervals

algeria_training_rmse <- algeria_models |>
  accuracy() |>
  select(.model, RMSE)

algeria_manual_intervals <- algeria_fc |>
  as_tibble() |>
  filter(Year == min(Year)) |>
  select(.model, Year, .mean) |>
  left_join(algeria_training_rmse, by = ".model") |>
  mutate(
    Lower95 = .mean - 1.96 * RMSE,
    Upper95 = .mean + 1.96 * RMSE
  )

knitr::kable(algeria_manual_intervals, digits = 4)
.model Year .mean RMSE Lower95 Upper95
SES 2018 22.4447 5.8653 10.9487 33.9406
Holt 2018 21.9781 5.8886 10.4365 33.5198 
algeria_fc |>
  filter(Year == min(Year)) |>
  hilo(level = 95)
## # A tsibble: 2 x 6 [1Y]
## # Key:       Country, .model [2]
##   Country .model  Year
##   <fct>   <chr>  <dbl>
## 1 Algeria SES     2018
## 2 Algeria Holt    2018
## # ℹ 3 more variables: Exports <dist>, .mean <dbl>, `95%` <hilo>

For the first forecast year, 2018, the manual 95% prediction intervals are [10.9487, 33.9406] for SES and [10.4365, 33.5198] for Holt, measured as exports as a percentage of GDP.

R produces [10.7455, 34.1439] for SES and [10.0169, 33.9394] for Holt. Both intervals are slightly wider than the manual intervals because the fitted models use parameter-adjusted estimates of error variance, while the manual calculation uses training RMSE. R also uses the exact normal quantile rather than the rounded multiplier 1.96.

The intervals overlap substantially, indicating considerable uncertainty relative to the small difference between the two point forecasts.

Exercise 8.6: Forecasting Chinese GDP

Examine the GDP Series

china_economy <- global_economy |>
  filter(Country == "China")

china_economy |>
  autoplot(GDP) +
  labs(
    title = "China: Annual GDP",
    x = "Year",
    y = "GDP (current US dollars)"
  )

China’s GDP shows a strong, nonlinear upward trend from 1960 to 2017. Annual increases are relatively small early in the series and become much larger after 2000. Growth slows near the end before another increase in the final year.

The data are annual, so within-year seasonality cannot be assessed. A Box–Cox transformation may make the growth pattern easier to model, while a damped trend can reduce the extent to which the estimated trend is extrapolated into the distant future. Both options will be demonstrated in the following forecasts.

Exercise 8.6: Fit Models With and Without Damping and Transformation

china_models <- china_economy |>
  model(
    Holt = ETS(
      GDP ~ error("A") + trend("A") + season("N")
    ),
    Damped = ETS(
      GDP ~ error("A") + trend("Ad") + season("N")
    ),
    BoxCox = ETS(
      box_cox(GDP, 0) ~ error("A") + trend("A") + season("N")
    ),
    BoxCox_Damped = ETS(
      box_cox(GDP, 0) ~ error("A") + trend("Ad") + season("N")
    )
  )

china_models
## # A mable: 1 x 5
## # Key:     Country [1]
##   Country         Holt        Damped       BoxCox BoxCox_Damped
##   <fct>        <model>       <model>      <model>       <model>
## 1 China   <ETS(A,A,N)> <ETS(A,Ad,N)> <ETS(A,A,N)> <ETS(A,Ad,N)>

These models compare an additive linear trend with a damped trend, both on the original GDP scale and after a Box–Cox transformation with \(\lambda = 0\).

Exercise 8.6: Compare Twenty-Year Forecasts

china_fc <- china_models |>
  forecast(h = 20)

china_fc |>
  autoplot(china_economy, level = NULL) +
  labs(
    title = "China GDP: Damping and Box–Cox Forecast Comparison",
    x = "Year",
    y = "GDP (current US dollars)"
  ) +
  guides(colour = guide_legend(title = "Model"))

The forecasts differ substantially over the twenty-year horizon. Holt’s method on the original GDP scale projects approximately constant annual dollar increases. Its damped version progressively reduces those increases and produces a lower forecast.

With Box–Cox λ = 0, the models fit log GDP. A linear trend on the log scale corresponds to approximately constant percentage growth, producing an upward-curving forecast after transformation back to dollars. The undamped Box–Cox model produces the highest forecasts, while damping reduces its projected growth.

Damping limits trend extrapolation, whereas transformation changes how growth is represented. The large differences show how strongly long-term forecasts depend on these assumptions. This plot alone cannot establish which model is most accurate.

Exercise 8.7: Australian Gas Production

Examine Trend and Seasonality

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

gas_data |>
  autoplot(Gas) +
  labs(
    title = "Australian Quarterly Gas Production",
    x = "Year",
    y = "Gas production (petajoules)"
  )

Australian gas production shows a strong upward trend and a repeating annual seasonal pattern in the quarterly observations. Seasonal fluctuations are small when production is low and become much larger as production increases.

Multiplicative seasonality is appropriate because the seasonal effect changes with the level of the series. It represents seasonal variation proportionally, whereas additive seasonality assumes approximately constant-sized seasonal fluctuations.

Exercise 8.7: Fit Multiplicative Seasonal Models

gas_models <- gas_data |>
  model(
    HoltWinters = ETS(
      Gas ~ error("M") + trend("A") + season("M")
    ),
    Damped = ETS(
      Gas ~ error("M") + trend("Ad") + season("M")
    )
  )

gas_models
## # A mable: 1 x 2
##    HoltWinters        Damped
##        <model>       <model>
## 1 <ETS(M,A,M)> <ETS(M,Ad,M)>

Exercise 8.7: Compare Three-Year Forecasts

gas_fc <- gas_models |>
  forecast(h = "3 years")

gas_fc |>
  autoplot(gas_data, level = NULL) +
  labs(
    title = "Australian Gas Production: Three-Year Forecasts",
    x = "Year",
    y = "Gas production (petajoules)"
  ) +
  guides(colour = guide_legend(title = "Model"))

Both models forecast continued quarterly seasonal fluctuations over the next three years. The damped and undamped forecast paths are very similar and nearly overlap on this plot. Damping therefore makes little visible difference over this forecast horizon. However, the plot alone does not establish whether damping improves forecast accuracy; this requires comparison against held-out observations.

Exercise 8.7: Does Damping Improve Forecast Accuracy?

gas_train <- gas_data |>
  slice_head(n = nrow(gas_data) - 12)

gas_test_models <- gas_train |>
  model(
    HoltWinters = ETS(
      Gas ~ error("M") + trend("A") + season("M")
    ),
    Damped = ETS(
      Gas ~ error("M") + trend("Ad") + season("M")
    )
  )

gas_test_fc <- gas_test_models |>
  forecast(h = 12)

gas_test_fc |>
  accuracy(gas_data) |>
  select(.model, .type, RMSE, MAE, MAPE) |>
  knitr::kable(digits = 4)
.model .type RMSE MAE MAPE
Damped Test 23.3370 21.0716 9.3567
HoltWinters Test 24.5081 21.9748 9.7585

The damped model performed better over the final 12 quarters. Its test RMSE was 23.3370, compared with 24.5081 for the undamped Holt-Winters model. It also had lower MAE (21.0716 versus 21.9748) and MAPE (9.3567% versus 9.7585%).

These results support using the damped model for this series, although the improvement is modest and applies to this particular three-year test period.

Exercise 8.8: Holt-Winters Models for Retail Turnover

Exercise 8.8(a): Examine the Seasonal Pattern

retail_data <- aus_retail |>
  filter(
    State == "New South Wales",
    Industry == "Takeaway food services"
  )

retail_data |>
  autoplot(Turnover) +
  labs(
    title = "New South Wales Takeaway Food Services",
    x = "Year",
    y = "Turnover (millions of Australian dollars)"
  )

New South Wales takeaway food turnover shows an overall upward trend, with some periods of decline. Seasonal fluctuations are relatively small early in the series and become larger as turnover increases.

Multiplicative seasonality is appropriate because the size of the seasonal fluctuations changes with the level of turnover. It models seasonal effects proportionally, whereas additive seasonality assumes fluctuations of approximately constant size.

Exercise 8.8(b): Fit Holt-Winters Models

retail_models <- retail_data |>
  model(
    HoltWinters = ETS(
      Turnover ~ error("M") + trend("A") + season("M")
    ),
    Damped = ETS(
      Turnover ~ error("M") + trend("Ad") + season("M")
    )
  )

retail_models
## # A mable: 1 x 4
## # Key:     State, Industry [1]
##   State           Industry                HoltWinters        Damped
##   <chr>           <chr>                       <model>       <model>
## 1 New South Wales Takeaway food services <ETS(M,A,M)> <ETS(M,Ad,M)>

Exercise 8.8(c): Compare One-Step Forecast RMSE

retail_training_accuracy <- retail_models |>
  accuracy()

retail_training_accuracy |>
  select(.model, .type, RMSE) |>
  knitr::kable(digits = 4)
.model .type RMSE
HoltWinters Training 9.7404
Damped Training 9.7338

The damped model has a one-step training RMSE of 9.7338 million Australian dollars, compared with 9.7404 for the undamped Holt-Winters model. I slightly prefer the damped method because it has the lower RMSE and reduces the continuation of the estimated trend at longer forecast horizons.

However, the difference in training RMSE is very small. This comparison alone does not establish that the damped model will forecast unseen observations more accurately.

Exercise 8.8(d): Check the Residuals

retail_models |>
  select(Damped) |>
  gg_tsresiduals()

### Exercise 8.8(d): Ljung–Box Test

retail_models |>
  select(Damped) |>
  augment() |>
  features(.innov, ljung_box, lag = 24, dof = 0) |>
  select(.model, lb_stat, lb_pvalue) |>
  knitr::kable(digits = 6)
.model lb_stat lb_pvalue
Damped 25.25887 0.391805

The innovation residuals fluctuate around zero without an obvious trend. Most ACF values lie within the dashed bounds, although a few spikes extend beyond them.

The Ljung–Box test at 24 lags gives a p-value of 0.391805. At the 5% significance level, we fail to reject the null hypothesis of no autocorrelation through these lags. Together, the plot and test suggest that the residuals are reasonably consistent with white noise, although this does not prove independence.

Exercise 8.8(e): Compare Test-Set Accuracy

retail_train <- retail_data |>
  filter(year(Month) <= 2010)

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

retail_test_models <- retail_train |>
  model(
    HoltWinters = ETS(
      Turnover ~ error("M") + trend("A") + season("M")
    ),
    Damped = ETS(
      Turnover ~ error("M") + trend("Ad") + season("M")
    ),
    SeasonalNaive = SNAIVE(Turnover)
  )

retail_test_fc <- retail_test_models |>
  forecast(new_data = retail_test)

retail_test_accuracy <- retail_test_fc |>
  accuracy(retail_data)

retail_test_accuracy |>
  select(.model, .type, RMSE, MAE, MAPE) |>
  knitr::kable(digits = 4)
.model .type RMSE MAE MAPE
Damped Test 69.0235 62.7774 14.7334
HoltWinters Test 70.3640 60.8845 15.2879
SeasonalNaive Test 96.8019 79.5469 16.2993

Using data through December 2010 for training, the damped model achieved the lowest test RMSE: 69.0235 million Australian dollars. The undamped Holt-Winters model had an RMSE of 70.3640, while seasonal naïve had an RMSE of 96.8019. Both Holt-Winters methods therefore improved on the seasonal naïve benchmark.

I prefer the damped model based on the exercise’s RMSE criterion. It also had the lowest MAPE, although the undamped model had a slightly lower MAE. The small RMSE difference between the two Holt-Winters models contrasts with their larger improvement over seasonal naïve.

Exercise 8.9: Box–Cox Transformation, STL and ETS

Estimate the Box–Cox Transformation Parameter

retail_lambda <- retail_train |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

retail_lambda
## [1] -0.1011142

Fit STL and ETS to the Transformed Training Data

retail_stl_fit <- retail_train |>
  model(
    STL_ETS = decomposition_model(
      STL(
        box_cox(Turnover, retail_lambda) ~
          season(window = "periodic")
      ),
      ETS(season_adjust ~ season("N"))
    )
  )

retail_stl_fit
## # A mable: 1 x 3
## # Key:     State, Industry [1]
##   State           Industry                                 STL_ETS
##   <chr>           <chr>                                    <model>
## 1 New South Wales Takeaway food services <STL decomposition model>

Compare Test-Set Forecast Accuracy

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

retail_stl_accuracy <- retail_stl_fc |>
  accuracy(retail_data)

bind_rows(
  retail_test_accuracy,
  retail_stl_accuracy
) |>
  select(.model, .type, RMSE, MAE, MAPE) |>
  arrange(RMSE) |>
  knitr::kable(digits = 4)
.model .type RMSE MAE MAPE
Damped Test 69.0235 62.7774 14.7334
HoltWinters Test 70.3640 60.8845 15.2879
STL_ETS Test 94.4427 90.7224 20.9947
SeasonalNaive Test 96.8019 79.5469 16.2993

The Box–Cox transformation parameter, estimated from the training data, was approximately −0.1011. After transforming the series, I applied STL with a fixed seasonal pattern and fitted a non-seasonal ETS model to the seasonally adjusted component.

On the same test period, this approach produced an RMSE of 94.4427 million Australian dollars, compared with 69.0235 for damped Holt-Winters. Its MAE and MAPE were also higher than those of both Holt-Winters models. Therefore, this Box–Cox/STL–ETS specification did not improve on the best previous forecasts.

Although its RMSE was slightly lower than seasonal naïve’s 96.8019, its MAE and MAPE were worse. Damped Holt-Winters remains my preferred method based on test RMSE.

Exercise 8.9: Compare Forecasts with Actual Turnover

bind_rows(
  retail_test_fc |> filter(.model == "Damped"),
  retail_stl_fc
) |>
  autoplot(
    retail_data |> filter(year(Month) >= 2005),
    level = NULL
  ) +
  labs(
    title = "Retail Turnover: Damped Holt-Winters versus STL–ETS",
    x = "Year",
    y = "Turnover (millions of Australian dollars)"
  ) +
  guides(colour = guide_legend(title = "Model"))

Both methods retain the seasonal pattern, but their trend forecasts differ. STL–ETS projects stronger growth and overestimates turnover for much of the test period. The damped Holt-Winters forecasts gradually flatten, underestimating turnover toward the end.

Neither method captures the entire test-period trajectory. However, damped Holt-Winters has the lower overall test RMSE, supporting its selection over this STL–ETS specification.