Exponential Smoothing Exploration
Tai Chou-Kudu
2025-10-09
Exercise 8.1: Victorian Pig Slaughter
Part a: Simple Exponential Smoothing
pigs <- aus_livestock %>%
filter(Animal == "Pigs", State == "Victoria")
fit_pigs <- 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.3221247
##
## Initial states:
## l[0]
## 100646.6
##
## sigma^2: 87480760
##
## AIC AICc BIC
## 13737.10 13737.14 13750.07fc_pigs <- fit_pigs %>%
forecast(h = 4)
fc_pigs %>%
autoplot(pigs %>% filter(year(Month) >= 2010)) +
labs(title = "Victorian Pig Slaughter: 4-Month Forecast",
subtitle = "Simple Exponential Smoothing (ETS-ANN)",
y = "Thousands of Pigs", x = "Month") +
theme_minimal()
fc_pigsThe optimal smoothing parameter alpha = 0.32 shows the model is balancing recent observations with historical data pretty evenly. The initial level l_0 starts around 100.6 thousand pigs. The forecasts stay fairly stable around 95-100 thousand, which makes sense given how the series has been fluctuating lately. You can see the forecast intervals getting wider as we go further out, which reflects growing uncertainty.
Part b: Manual Prediction Interval
residuals_pigs <- augment(fit_pigs)
s <- sd(residuals_pigs$.resid, na.rm = TRUE)
first_forecast <- fc_pigs$.mean[1]
manual_lower <- first_forecast - 1.96 * s
manual_upper <- first_forecast + 1.96 * s
cat("Manual 95% Prediction Interval:\n")## Manual 95% Prediction Interval:cat(sprintf("Point Forecast: %.2f\n", first_forecast))## Point Forecast: 95186.56cat(sprintf("Standard Deviation (s): %.2f\n", s))## Standard Deviation (s): 9344.67cat(sprintf("Manual PI: [%.2f, %.2f]\n\n", manual_lower, manual_upper))## Manual PI: [76871.01, 113502.10]cat("R's 95% Prediction Interval:\n")## R's 95% Prediction Interval:print(hilo(fc_pigs, level = 95)[1,])## # 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 calculation using ŷ ± 1.96s gives us an interval around [76.9, 113.5] thousand pigs, while R’s interval is a bit narrower. This happens because R accounts for uncertainty in the smoothing parameters themselves, not just the residuals. Our simpler manual method only uses the residual standard deviation, which ends up being more conservative.
Exercise 8.5: Global Economy Exports
aus_exports <- global_economy %>%
filter(Country == "Australia") %>%
select(Year, Exports)Part a: Data Exploration
aus_exports %>%
autoplot(Exports) +
labs(title = "Australian Exports as Percentage of GDP",
y = "Exports (% of GDP)", x = "Year") +
theme_minimal()
The Australian exports data has some interesting patterns. There’s a general upward trend from the 1960s through the early 2010s, with exports climbing from around 12-14% of GDP to over 22%. The series is pretty volatile, especially in the earlier years. You can see some cyclical patterns that probably match up with commodity price cycles and global economic conditions. After 2010, things seem to flatten out or even decline a bit, which might reflect changes in trade patterns or commodity markets.
Parts b & c: ETS(A,N,N) Model
fit_ann <- aus_exports %>%
model(ETS_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.5756fc_ann <- fit_ann %>%
forecast(h = 10)
fc_ann %>%
autoplot(aus_exports) +
labs(title = "Australian Exports: Level-Only Forecast",
subtitle = "ETS(A,N,N) - No Trend Component",
y = "Exports (% of GDP)", x = "Year") +
theme_minimal()
accuracy(fit_ann)The ETS(A,N,N) model treats the data as moving around a slowly changing level without any trend. With an RMSE around 1.15 percentage points, it captures the general movements but isn’t perfect. The flat forecasts make sense since there’s no trend component.
Part d: Comparison with ETS(A,A,N)
fit_aan <- aus_exports %>%
model(ETS_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.6146accuracy_comparison <- bind_rows(
accuracy(fit_ann) %>% mutate(Model = "ETS(A,N,N)"),
accuracy(fit_aan) %>% mutate(Model = "ETS(A,A,N)")
)
accuracy_comparison %>%
select(Model, RMSE, MAE)cat("\nAICc Comparison:\n")##
## AICc Comparison:glance(fit_ann) %>% select(AICc)glance(fit_aan) %>% select(AICc)The ETS(A,A,N) model adds a trend component, which means one extra parameter. Looking at the comparison, the simpler model actually has a slightly lower AICc (258 vs 259), suggesting it fits marginally better when accounting for the extra parameter in the trended model. The improvement is pretty small though, and the RMSE values are less than 0.1 different in value, showing that while there’s been some historical growth, the trend isn’t super strong or persistent.
Part e: Forecast Comparison
fc_aan <- fit_aan %>%
forecast(h = 10)
bind_rows(
fc_ann %>% mutate(Model = "No Trend"),
fc_aan %>% mutate(Model = "With Trend")
) %>%
autoplot(aus_exports %>% filter(Year >= 1990), level = NULL) +
labs(title = "Comparing Level vs. Trend Forecasts",
subtitle = "Australian Exports: Which approach is more reasonable?",
y = "Exports (% of GDP)", x = "Year") +
theme_minimal()
The forecasts are quite different. The ETS(A,N,N) model shows stable exports around current levels, while ETS(A,A,N) keeps declining. If you look at recent data (post-2010), exports have been relatively flat with maybe some decline. I think the level-only model makes more sense here because the historical upward trend seems to have ended, and it doesn’t make sense to extrapolate a downward trend forever for a major exporting country like Australia.
Part f: Manual Prediction Intervals
resid_ann <- augment(fit_ann)
rmse_ann <- sqrt(mean(resid_ann$.resid^2, na.rm = TRUE))
first_fc_ann <- fc_ann$.mean[1]
cat("ETS(A,N,N) - First Forecast:\n")## ETS(A,N,N) - First Forecast:cat(sprintf("Manual 95%% PI: [%.3f, %.3f]\n",
first_fc_ann - 1.96 * rmse_ann,
first_fc_ann + 1.96 * rmse_ann))## Manual 95% PI: [18.359, 22.855]cat("R's Interval:\n")## R's Interval:print(hilo(fc_ann, level = 95)[1,])## # A tsibble: 1 x 5 [1Y]
## # Key: .model [1]
## .model Year Exports
## <chr> <dbl> <dist>
## 1 ETS_A… 2018 N(21, 1.4)
## # ℹ 2 more variables: .mean <dbl>, `95%` <hilo>resid_aan <- augment(fit_aan)
rmse_aan <- sqrt(mean(resid_aan$.resid^2, na.rm = TRUE))
first_fc_aan <- fc_aan$.mean[1]
cat("\nETS(A,A,N) - First Forecast:\n")##
## ETS(A,A,N) - First Forecast:cat(sprintf("Manual 95%% PI: [%.3f, %.3f]\n",
first_fc_aan - 1.96 * rmse_aan,
first_fc_aan + 1.96 * rmse_aan))## Manual 95% PI: [18.650, 23.027]cat("R's Interval:\n")## R's Interval:print(hilo(fc_aan, level = 95)[1,])## # A tsibble: 1 x 5 [1Y]
## # Key: .model [1]
## .model Year Exports
## <chr> <dbl> <dist>
## 1 ETS_A… 2018 N(21, 1.3)
## # ℹ 2 more variables: .mean <dbl>, `95%` <hilo>Again, the manual intervals are a bit wider than R’s. This pattern is consistent across both models. R’s intervals factor in the uncertainty around parameter estimates through formulas or simulation, while our manual method only uses the residual standard error. Both show pretty substantial uncertainty - the intervals span roughly 7-8 percentage points.
Exercise 8.6: Chinese GDP Forecasting
china_gdp <- global_economy %>%
filter(Country == "China") %>%
select(Year, GDP)
lambda_china <- china_gdp %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
fit_china <- china_gdp %>%
model(
Auto = ETS(GDP),
Damped = ETS(GDP ~ trend("Ad")),
BoxCox = ETS(box_cox(GDP, lambda_china)),
BoxCox_Damped = ETS(box_cox(GDP, lambda_china) ~ trend("Ad"))
)
glance(fit_china) %>%
select(.model, AIC, AICc, BIC)fc_china <- fit_china %>%
forecast(h = 30)
fc_china %>%
autoplot(china_gdp %>% filter(Year >= 1980), level = NULL) +
facet_wrap(~.model, scales = "free_y") +
labs(title = "Chinese GDP: Impact of Model Specification",
subtitle = "30-Year Forecasts Under Different Assumptions",
y = "GDP (USD)", x = "Year") +
theme_minimal()
This exercise really shows how much model choice matters for long-term forecasts. The automatic selection picks a multiplicative model with trend, giving us exponential growth that just keeps going. The damped trend version curves the growth down, avoiding those explosive long-term forecasts while still showing solid increases.
The Box-Cox transformation stabilizes the variance and makes growth more linear on the transformed scale, which translates to more moderate growth when you convert back. Combining Box-Cox with damping gives the most conservative forecasts, with growth eventually leveling off.
Given that China’s growth has been slowing down recently, the Box-Cox damped approach seems most realistic. The undamped models would have China’s economy at unrealistically high levels by 2050. This shows an important point: exponential smoothing with trend should usually be damped for long-term forecasts unless you have strong reason to believe exponential growth will continue forever. This can be true for an ETF like VGT, growing near 15% nowadays, but most likely not forever.
Exercise 8.7: Gas Production
gas_data <- aus_production %>%
select(Quarter, Gas)
gas_data %>%
autoplot(Gas) +
labs(title = "Australian Gas Production Over Time",
y = "Gas Production (Petajoules)", x = "Quarter") +
theme_minimal()
fit_gas <- gas_data %>%
model(
Multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
Mult_Damped = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)
report(fit_gas)glance(fit_gas) %>%
select(.model, AICc, BIC)fc_gas <- fit_gas %>%
forecast(h = 16)
fc_gas %>%
autoplot(gas_data %>% filter(year(Quarter) >= 1990), level = 80) +
labs(title = "Gas Production Forecasts: Damping Effect",
subtitle = "Comparing Regular vs. Damped Trend (4 Years Ahead)",
y = "Gas Production (Petajoules)", x = "Quarter") +
theme_minimal()
accuracy(fit_gas)Why multiplicative seasonality?
Looking at the plot, the seasonal ups and downs clearly get bigger as production increases. In the early years (pre-1990), the seasonal swings are pretty small, but by the 2000s, those same seasonal patterns create much larger differences. This proportional relationship between level and seasonal variation is what defines multiplicative seasonality. An additive model would underestimate seasonality when production is high and overestimate it when production is low.
Does damping improve forecasts?
The damped model has a slightly lower AICc, so there’s some improvement. More importantly, it gives more reasonable long-term forecasts that don’t assume the trend will keep going forever. Since gas production is limited by infrastructure, reserves, and demand, the damped approach makes more sense. The regular trend forecast shows production climbing aggressively, while the damped version lets it level off.
Exercise 8.8: Retail Time Series
set.seed(12345678)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
cat("Selected Series:\n")## Selected Series:print(unique(myseries[,c("State", "Industry")]))## # A tibble: 1 × 2
## State Industry
## <chr> <chr>
## 1 Northern Territory Clothing, footwear and personal accessory retailingmyseries %>%
autoplot(Turnover) +
labs(title = paste("Retail Turnover:", unique(myseries$Industry)),
subtitle = unique(myseries$State),
y = "Turnover ($Million AUD)", x = "Month") +
theme_minimal()
Part a: Why Multiplicative Seasonality?
You can clearly see that seasonal variation isn’t constant - it grows with the level of the series. In the early years, seasonal peaks and troughs are smaller, but as turnover increases, the seasonal swings get much larger in absolute terms. But if you look at percentage changes, they stay relatively stable. This is what defines multiplicative seasonality - the seasonal effect is a percentage of the current level rather than a fixed amount.
Part b: Holt-Winters Multiplicative Method
fit_retail <- myseries %>%
model(
HW_Mult = ETS(Turnover ~ error("M") + trend("A") + season("M")),
HW_Damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)
report(fit_retail)Both models capture the multiplicative seasonal pattern. The damped version adds flexibility by allowing the trend to weaken over time.
Part c: RMSE Comparison
accuracy(fit_retail) %>%
select(.model, RMSE, MAE) %>%
arrange(RMSE)glance(fit_retail) %>%
select(.model, AICc)The models perform very similarly. The regular multiplicative model has a slightly lower AICc (1776 vs 1783), but the difference is minimal and both models have nearly identical RMSE values. I prefer the damped model because retail trends rarely keep going forever - markets mature, competition changes, and what people want shifts over time.
Part d: Residual Diagnostics
best_model <- fit_retail %>% select(HW_Damped)
best_model %>% gg_tsresiduals()
augment(best_model) %>%
features(.innov, ljung_box, lag = 24, dof = 0)The residuals look reasonably close to white noise. Most lags in the ACF are within the confidence bounds, though there might be some minor leftover autocorrelation. The histogram looks roughly normal, maybe with slightly heavier tails. The Ljung-Box test helps us formally check if the residuals are white noise - if the p-value is above 0.05, the model has done a good job capturing the patterns in the data.
Part e: Test Set Performance
train <- myseries %>%
filter(year(Month) <= 2010)
test <- myseries %>%
filter(year(Month) > 2010)
fit_train <- train %>%
model(
HW_Mult = ETS(Turnover ~ error("M") + trend("A") + season("M")),
HW_Damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
SNAIVE = SNAIVE(Turnover)
)
fc_test <- fit_train %>%
forecast(h = nrow(test))
accuracy(fc_test, myseries) %>%
select(.model, RMSE, MAE, MAPE) %>%
arrange(RMSE)fc_test %>%
autoplot(myseries %>% filter(year(Month) >= 2005), level = NULL) +
labs(title = "Test Set Forecasts: Can We Beat Seasonal Naive?",
y = "Turnover ($Million AUD)", x = "Month") +
theme_minimal()
Looking at the test set results, we can see which method actually works better on new data. The seasonal naive method is a tough benchmark since it just uses the same month from last year. If our ETS models have lower RMSE on the test set, we’ve beaten this simple approach by capturing trend dynamics. The plot shows how well each method tracks the actual future data.
Exercise 8.9: STL Decomposition with ETS
lambda <- myseries %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
cat(sprintf("Guerrero's optimal lambda: %.4f\n", lambda))## Guerrero's optimal lambda: 0.0830fit_stl <- train %>%
model(
STL_ETS = decomposition_model(
STL(box_cox(Turnover, lambda) ~ season(window = 13)),
ETS(season_adjust ~ error("A") + trend("A") + season("N"))
)
)
report(fit_stl)## Series: Turnover
## Model: STL decomposition model
## Transformation: box_cox(Turnover, lambda)
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ETS(A,A,N)
## Smoothing parameters:
## alpha = 0.6146768
## beta = 0.0001000001
##
## Initial states:
## l[0] b[0]
## 1.008946 0.006661319
##
## sigma^2: 0.0053
##
## AIC AICc BIC
## 109.0046 109.2293 127.0519
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 0fit_all <- train %>%
model(
HW_Mult = ETS(Turnover ~ error("M") + trend("A") + season("M")),
HW_Damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
SNAIVE = SNAIVE(Turnover),
STL_ETS = decomposition_model(
STL(box_cox(Turnover, lambda) ~ season(window = 13)),
ETS(season_adjust ~ error("A") + trend("A") + season("N"))
)
)
fc_all <- fit_all %>%
forecast(h = nrow(test))
accuracy(fc_all, myseries) %>%
select(.model, RMSE, MAE, MAPE) %>%
arrange(RMSE)fc_all %>%
autoplot(myseries %>% filter(year(Month) >= 2008), level = NULL) +
labs(title = "Final Showdown: All Methods Compared",
subtitle = "Which approach best captures future retail dynamics?",
y = "Turnover ($Million AUD)", x = "Month") +
theme_minimal()
The STL+ETS approach works differently. By first decomposing the Box-Cox transformed series and then modeling only the seasonally adjusted part, we’re assuming the seasonal pattern is stable and can be estimated separately. This works well when seasonality is very regular and doesn’t interact much with the trend.
Looking at all methods on the test set, we can see which one actually performs best for this series. The STL approach sometimes wins because it handles the seasonal pattern more robustly, but it can fail if the seasonal pattern changes over time. The results depend on which retail series you picked, but usually the ETS multiplicative models do well because they let seasonality scale naturally with the level.
What’s interesting here is seeing how different approaches translate to actual forecast accuracy. The direct multiplicative ETS approach models everything together, while STL splits things up first. Neither is always better - it depends on your data.