Introduction

The following RMD contains answers to exercises in chapter 8 of the Forecasting: Principles and Practice textbook for CUNY SPS DATA 624 Spring 2025 https://otexts.com/fpp3/expsmooth-exercises.html. This chapter focuses on exponential smoothing. The 8.8 Exercises answered here include 8.1, 8.5, 8.6, 8.7, 8.8, and 8.9.

Load Libraries

library(fpp3)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(tidyr)

Q1

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

  1. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of α and ℓ0, and generate forecasts for the next four months.

  2. Compute a 95% prediction interval for the first forecast using y ± 1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

# Save the data to pig_df
pig_df <- aus_livestock |>
  filter(Animal == "Pigs" & State == "Victoria")

# Visualize the data
pig_plot <- pig_df |>
  autoplot(Count) +
  labs(title = 'Counnt of Pigs Slaughtered in Victoria')
pig_plot

# Use ETS() to estimate the equivalent model
pig_model <- pig_df |>
  model(ses = ETS(Count ~ error('A') + trend('N') + season('N')))

# Get optimal values of alpha and l[0]
report(pig_model)
## 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
# Get forecast for the next 4 months
forecast(pig_model, h = 4)
## # 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>
# Visualize the forecast, with a filter on Month to see the forecast more clearly 
pig_model_plot <- pig_model |>
  forecast(h = 4) |>
  autoplot(filter(pig_df, Month >= yearmonth("2010 Jan"))) +
  labs(title = "Four Month Forecast Data of Pigs Slaughtered in Victoria")
pig_model_plot

Optimal values for alpha and l[0]:

alpha = 0.3221247 Interpretation of this alpha: The alpha can be between 0 and 1, so the alpha of 0.322 means the data has a relatively small influence on the forecast.

l[0] = 100,646.6 Interpretation of this l[0]: The series has an estimated starting average of 100,646.6.

The graph Four Month Forecast Data of Pigs Slaughtered in Victoria shows the forecast for 4 months in the future of the data. The ETS model does not account for trend, so the forecast looks like it may be slightly lower than the upward trend may suggest.

# Get the first forecast row
yHat <- pig_model |> 
  forecast(h = 4) |>
  pull(Count) |>
  head(1)

# Get the standard deviation of the residuals
sdev <- residuals(pig_model)$.resid |>
  sd()

# Calculate the lower and upper confidence intervals using y ± 1.96s 
lowCi <- yHat - 1.96 * sdev
highCi <- yHat + 1.96 * sdev
results <- c(lowCi, highCi)
results
## <distribution[2]>
## [1] N(76871, 8.7e+07)  N(113502, 8.7e+07)

The 95% confidence interval is [76,871, 113,502]. This range contains the actual number of pigs slaughtered an expected 95% of the time.

Q5

Prompt: 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.
# Get Luxembourg data
lux_df <- global_economy |>
  filter(Country == "Luxembourg")

# Plot using autoplot
lux_df |> 
  autoplot(Exports) +
  labs(y = "% of GDP", title = 'Luxembourg Exports')

The data here shows exports measured by the percentage of GDP. The overall trend is upward, especially increasing after 1980.

# Use an ETS(A,N,N) model to forecast the series
lux_model <- lux_df |> 
  model(ETS(Exports ~ error('A') + trend("N") + season('N')))

# Forecast the Luxembourg data to h=4
lux_fcst <- lux_model |> 
  forecast(h = 4)
lux_fcst
## # A fable: 4 x 5 [1Y]
## # Key:     Country, .model [1]
##   Country    .model                                                        Year
##   <fct>      <chr>                                                        <dbl>
## 1 Luxembourg "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))"  2018
## 2 Luxembourg "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))"  2019
## 3 Luxembourg "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))"  2020
## 4 Luxembourg "ETS(Exports ~ error(\"A\") + trend(\"N\") + season(\"N\"))"  2021
## # ℹ 2 more variables: Exports <dist>, .mean <dbl>
# Plot the forecast
lux_fcst |> 
  autoplot(lux_df) +
  labs(y = "% of GDP", title = 'Luxembourg Exports')

# Compute the RMSE values for the training data
accuracy(lux_model)
## # A tibble: 1 × 11
##   Country    .model        .type    ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <fct>      <chr>         <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Luxembourg "ETS(Exports… Trai…  2.44  7.54  5.57  1.48  4.45 0.983 0.991 0.116

The RMSE of the training data is 7.544208.

# Compare the results to those from an ETS(A,A,N) model
lux_model_comparison <- lux_df |>
  model(
    ANN = ETS(Exports ~ error('A') + trend('N') + season('N')),
    AAN = ETS(Exports ~ error('A') + trend('A') + season('N'))
  )

# Show the accuracies
accuracy(lux_model_comparison)
## # A tibble: 2 × 11
##   Country    .model .type       ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
##   <fct>      <chr>  <chr>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 Luxembourg ANN    Training  2.44  7.54  5.57 1.48   4.45 0.983 0.991 0.116 
## 2 Luxembourg AAN    Training  1.12  7.10  5.13 0.633  4.25 0.904 0.933 0.0252

The ANN model has a slightly higher RMSE, meaning the AAN model’s RMSE of 7.098 is lower and more accurate for this data.

# Compare the models
lux_model_comparison |>
  forecast(h = 4) |> 
  autoplot(lux_df, level = NULL) +
  labs(title = 'Luxembourg Annual Exports ANN Vs AAN Comparison')

Based on the trend of the graph, it looks like AAN is more accurate for predicting the model. The ANN prediction has a negative slope, which does not match the clear upward trend of the data.

# Get sdev of first forecast for AAN
lux_sdev <- lux_model_comparison |>
  select(Country, AAN) |>
  accuracy() |>
  transmute(Country, lux_sdev = RMSE)

# Graph the data
lux_model_comparison |>
  select(Country, AAN) |>
  forecast(h = 1) |>
  left_join(lux_sdev, by = 'Country') |>
  mutate(lowCi = Exports - 1.96 * lux_sdev,
         highCi = Exports + 1.96 * lux_sdev) |>
  select(Country, Exports, lowCi, highCi)
## # A fable: 1 x 5 [?]
## # Key:     Country [1]
##   Country       Exports
##   <fct>          <dist>
## 1 Luxembourg N(235, 54)
## # ℹ 3 more variables: lowCi <dist>, highCi <dist>, Year <dbl>

The confidence intervals for AAN calculated here show that the true value for export percentage should be in the (221, 249) range 95% of the time.

This is a much closer range than the ANN model, which expected (230, 289). The AAN model’s tighter range is preferable to the wider range of the ANN model here.

Q6

Prompt: 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.]

# Save chinese economic data
china_df <- global_economy |>
  filter(Country == "China")

# Visualize data
china_df |> autoplot(GDP) +
  labs(title = "Chinese GDP")

# Get optimized lambda for BoxCox
china_lambda <- china_df |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

# Try several ETS() options
china_ets_compare <- china_df |>
  model(
    Standard_ETS = ETS(GDP),
    BoxCox_ETS = ETS(box_cox(GDP, china_lambda)),
    Damped_ETS = ETS(GDP ~ trend('Ad', phi = 0.9)),
    Log_ETS = ETS(log(GDP))
  )

# Visualize using a large h
china_ets_compare |>
  forecast(h = 15) |>
  autoplot(china_df, level = NULL) +
  labs(title = 'Chinese GDP ETS Forecast Comparison')

The models vary based on their forecasting priorities. It looks like the BoxCox and Log ETS transformations are expecting a quick exponential increase in GDP, which is not common in an actual economy. The Damped and Standard ETS transformations looks much more feasable for the growth of a developed economy. Now, the data itself does show an exponential increase, as seen in the Chinese GDP visualization.

Q7

Prompt: 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?

# Get the gas data
gas_df <- aus_production |>
  select(Gas)

# Plot the data with autoplot
autoplot(gas_df, Gas)

# Get the gas data
gas_model <- gas_df |>
  model(fit = ETS(Gas))
report(gas_model)
## 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
# Visualize the model
gas_model |>
  forecast(h = 15) |>
  autoplot(filter(gas_df, Quarter >= yearquarter("1965 Q1")), level= NULL)

In this case, multiplicative seasonality is necessary due to the fact that the seasonality is increasing over time. I cropped the data to be 1965 and later to beter see the forecast data.

# Make the trend damped
gas_fit <- gas_df |>
  model(fit = ETS(Gas  ~ trend('Ad', phi = 0.9)))

# Visualize the data
gas_fit |>
  forecast(h = 15) %>%
  autoplot(filter(gas_df, Quarter >= yearquarter("1965 Q1")), level= NULL)

The forecast with the damped trajectory seems to have a similar seasonality variation with less of an upward trend. There still is a slightly upward trend, so I would say damping the forecast does make it more realistic.

Q8

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

a. Why is multiplicative seasonality necessary for this series?

b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

d. Check that the residuals from the best method look like white noise.

e. 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?

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

# Visualize the data
myseries |> 
  autoplot(Turnover)

Multiplicative seasonality is necessary for this series because the seasons change in variability over time. The variability of seasons seems to increase with the trend, until around 2011, where trend decreases and so does variability, until both decrease again starting in around 2014.

# Apply Holt-Winters’ multiplicative method to the data
myseries_model <- myseries |>
  model(
    'Holt Winters Multiplicative' = ETS(Turnover ~ error('M') + trend('A') + season('M')),
    'Holt Winters Damped' = ETS(Turnover ~ error('M') + trend('Ad') + season('M'))
  )

# Forecast the Holt-Winters
myseries_HW <- myseries_model |> forecast(h = 35)

# Visualize the data
myseries_HW |> autoplot(filter(myseries, Month >= yearmonth("2005 Jan")), level = NULL)

Again, the damped data trends upward less than the standard multiplicative method. In this case, it looks like the data will end up trending upward based on the overall trend of the graph so the damped method seems less effective here.

# Get RMSE accuracies
accuracy(myseries_model) |>
  select('.model', 'RMSE')
## # A tibble: 2 × 2
##   .model                       RMSE
##   <chr>                       <dbl>
## 1 Holt Winters Multiplicative  7.03
## 2 Holt Winters Damped          7.05

The RMSE is quite similar for both methods with a difference of around 0.02. Because this difference is so small, it is difficult to say that the multiplicative method is better at predicting than the corresponding damped method. I would consider these about the same as far as RMSE is concerned.

# Visualize the residuals of the Holt Winters data
myseries_model |>
  select('Holt Winters Multiplicative') |>
  gg_tsresiduals()

There are several outliers in the lag plot that would suggest there may be an issue with the model. Otherwise, the Innovation Residuals graph does appear mostly random, with increases in variability over time, espeicially in 2010. The residual histogram is lightly right skewed. These issues are signs of minor autocorrelation, but they may not be enough to flag the model as useless.

# Get data to the end of 2010.
myseries_train <- myseries |>
  filter(year(Month) <= 2010)

# Model the data
mst_model <- myseries_train |>
  model(
    'Holt Winters Multiplicative' = ETS(Turnover ~ error("M") + trend("A") + season("M")),
    'Holt Winters Damped' = ETS(Turnover ~ error('M') + trend('Ad') + season('M')),
    'Seasonal Naive' = SNAIVE(Turnover)
  )

# Get the comparison of the models forecasts
myseries_fcst <- mst_model |>
  forecast(new_data = anti_join(myseries, myseries_train, by = c("State", "Industry", "Series ID", "Month", "Turnover")))

# Visualize the data
myseries_fcst |>
  autoplot(filter(myseries, Month >= yearmonth("2005 Jan")), level = NULL)

# Get accuracy of the myseries training model
accuracy(mst_model) |>
  select(.type, .model, RMSE)
## # A tibble: 3 × 3
##   .type    .model                       RMSE
##   <chr>    <chr>                       <dbl>
## 1 Training Holt Winters Multiplicative  6.55
## 2 Training Holt Winters Damped          6.59
## 3 Training Seasonal Naive              11.6
# Get accuracy of the myseries forecast
myseries_fcst |>
  accuracy(myseries)  |>
  select(.type, .model, RMSE)
## # A tibble: 3 × 3
##   .type .model                       RMSE
##   <chr> <chr>                       <dbl>
## 1 Test  Holt Winters Damped          32.9
## 2 Test  Holt Winters Multiplicative  35.9
## 3 Test  Seasonal Naive               29.0

In the training data, the RMSE of Seasonal Naive is higher than that of the two Holt Winters models, suggesting the latter two models should be preferred.

In the testing of the models, the RMSE of the Seasonal Naive data is lower than that of the other models.

The differences between the training and testing results suggest that all three models perform relatively similarly in their predictive abilities.

Q9

Prompt: 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?

# Get the optimized lambda value for the BoxCox
myseries_lambda <- myseries_train |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

# Train the BoxCox
mst_boxcox <- myseries_train |>
  mutate(
    bc = box_cox(Turnover, myseries_lambda)
  )

# Model the STL and ETS BoxCoxes
mst_bc_model <- mst_boxcox |>
  model(
    'STL Box-Cox' = STL(bc ~ season(window = 'periodic'), robust = TRUE),
    'ETS Box-Cox' = ETS(bc)
  )

# Fit the model
mst_mp_best_fit <- mst_boxcox |>
  model(
    'Holt Winters Multiplicative' = ETS(Turnover ~ error('M') + trend('A') + season('M'))
  )

# Get STL and ETS BoxCox model accuracies
accuracy(mst_bc_model)
## # A tibble: 2 × 12
##   State  Industry .model .type      ME    RMSE     MAE     MPE  MAPE  MASE RMSSE
##   <chr>  <chr>    <chr>  <chr>   <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 Victo… Other s… STL B… Trai… 1.21e-4 0.00634 0.00504 4.79e-3 0.248 0.373 0.370
## 2 Victo… Other s… ETS B… Trai… 1.38e-5 0.00755 0.00593 1.97e-4 0.292 0.439 0.441
## # ℹ 1 more variable: ACF1 <dbl>
# Get Holt Winters model accuracy
accuracy(mst_mp_best_fit)
## # A tibble: 1 × 12
##   State    Industry     .model .type    ME  RMSE   MAE     MPE  MAPE  MASE RMSSE
##   <chr>    <chr>        <chr>  <chr> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 Victoria Other speci… Holt … Trai… 0.314  6.55  3.72 -0.0776  3.72 0.474 0.564
## # ℹ 1 more variable: ACF1 <dbl>

The RMSE values for the two Box-Cox types, STL and ETS are extremely low (<0.008). This makes them much more accurate than the Holt Winters model, which has a RMSE accuracy of 6.551.