DATA 624 - Assignment 5
Gullit Navarrete
2025-10-03
Objective:
Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. Please submit both the link to your Rpubs and the .pdf file with your run code
8.1
Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
a. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of alpha and ell0, and generate forecasts for the next four months.
library(fpp3)## Warning: package 'fpp3' was built under R version 4.4.3## Registered S3 method overwritten by 'tsibble':
## method from
## as_tibble.grouped_df dplyr## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──## ✔ tibble 3.2.1 ✔ tsibble 1.1.6
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.1 ✔ feasts 0.4.2
## ✔ lubridate 1.9.3 ✔ fable 0.4.1
## ✔ ggplot2 3.5.1## Warning: package 'tsibble' was built under R version 4.4.3## Warning: package 'tsibbledata' was built under R version 4.4.3## Warning: package 'feasts' was built under R version 4.4.3## Warning: package 'fabletools' was built under R version 4.4.3## Warning: package 'fable' was built under R version 4.4.3## ── 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()pigs_vic <- aus_livestock |>
filter(State == "Victoria", Animal=="Pigs") |>
select(Month, Count)
# Single Exponential Smoothing
fit <- pigs_vic |>
model(SES = ETS(Count ~ error("A") + trend("N") + season("N")))
# Optimal alpha and initial level (ℓ0)
report(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# Forecast next 4 months
fc <- fit |> forecast(h = 4)
autoplot(pigs_vic, Count) +
autolayer(fc, level = NULL) +
labs(title = "Victoria pigs: SES", x = NULL, y = "Count")
Compute a 95% prediction interval for the first forecast using y hat+-1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
s <- residuals(fit) |>
as_tibble() |>
summarise(s = sd(.resid, na.rm = TRUE)) |>
pull(s)
# Forecast point
yhat1 <- fc |> as_tibble() |> slice(1) |> pull(.mean)
lower_manual <- yhat1 - 1.96 * s
upper_manual <- yhat1 + 1.96 * s
fc1_r <- fc |>
hilo(95) |>
unpack_hilo(`95%`) |>
as_tibble() |>
slice(1)
lower_r <- fc1_r$`95%_lower`
upper_r <- fc1_r$`95%_upper`
tibble(
method = c("y hat+-1.96s", "R"),
lower = c(lower_manual, lower_r),
mean = c(yhat1, yhat1),
upper = c(upper_manual, upper_r)
)The two intervals are very similar as R’s interval is only a little wider with 16-17 units on each side than the manual y hat+-1.96s interval.
8.5
Data set global_economy contains the annual Exports from many countries. Select one country to analyse.
a. Plot the Exports series and discuss the main features of the data.
cntry <- "Australia"
exp_ts <- global_economy |>
filter(Country == cntry) |>
select(Year, Exports)
autoplot(exp_ts, Exports) +
labs(title = paste(cntry, "Exports"),
y = "Exports (GDP %)", x = NULL)
b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# Use the same 'exp_ts' from part (a)
fit_ann <- exp_ts |>
model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
fc_ann <- fit_ann |> forecast(h = 10)
autoplot(exp_ts, Exports) +
autolayer(fc_ann) +
labs(title = paste(cntry, "Exports ETS(A,N,N)"),
y = "Exports (% of GDP)", x = NULL)
c. Compute the RMSE values for the training data.
# RMSE on the training data
acc_train <- fit_ann |>
accuracy() |>
dplyr::select(.model, RMSE)
acc_traind. 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.
aus_exports <- global_economy |>
filter(Country == "Australia") |>
select(Year, Exports)
# Fit ETS with Holt's linear
fit_ann <- aus_exports |>
model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
fit_aann <- aus_exports |>
model(AAN = ETS(Exports ~ error("A") + trend("A") + season("N")))
# RMSE on training data
bind_rows(
fit_ann |> accuracy() |> dplyr::select(.model, RMSE),
fit_aann |> accuracy() |> dplyr::select(.model, RMSE)
)e. Compare the forecasts from both methods. Which do you think is best?
When comparing the forecasts, teh no-trend model ETS(A,N,N) just keeps the forecast flat at the last level, while the trend model ETS(A,A,N) continues the upward climb. Because the data clearly trends up while its error is a bit lower (1.12 vs 1.15) the trend model is the better pick.
f. 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.
fc_ann <- fit_ann |> forecast(h = 1)
fc_aann <- fit_aann |> forecast(h = 1)
# Training RMSE
s_ann <- fit_ann |> accuracy() |> dplyr::pull(RMSE)
s_aann <- fit_aann |> accuracy() |> dplyr::pull(RMSE)
# Point forecasts
yhat_ann <- fc_ann |> as_tibble() |> dplyr::slice(1) |> dplyr::pull(.mean)
yhat_aann <- fc_aann |> as_tibble() |> dplyr::slice(1) |> dplyr::pull(.mean)
# 95% PIs using +-1.96*s
man_ann <- c(lower = yhat_ann - 1.96*s_ann, upper = yhat_ann + 1.96*s_ann)
man_aann <- c(lower = yhat_aann - 1.96*s_aann, upper = yhat_aann + 1.96*s_aann)
# R's 95% PIs
r_ann <- fc_ann |>
hilo(95) |>
as_tibble() |>
dplyr::slice(1) |>
unpack_hilo(`95%`) |>
dplyr::transmute(lower = `95%_lower`, upper = `95%_upper`)
r_aann <- fc_aann |>
hilo(95) |>
as_tibble() |>
dplyr::slice(1) |>
unpack_hilo(`95%`) |>
dplyr::transmute(lower = `95%_lower`, upper = `95%_upper`)
# Final comparison
tibble::tibble(
model = c("ETS(A,N,N) manual", "ETS(A,N,N) R", "ETS(A,A,N) manual", "ETS(A,A,N) R"),
lower = c(man_ann["lower"], r_ann$lower, man_aann["lower"], r_aann$lower),
mean = c(yhat_ann, yhat_ann, yhat_aann, yhat_aann),
upper = c(man_ann["upper"], r_ann$upper, man_aann["upper"], r_aann$upper)
)Both are the same result as in the forecast is the same and the ranges are nearly identical. R’s intervals are just a tiny bit wider because it uses a fuller uncertainty calculation than our quick +-1.96 RMSE trick.
8.6
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.]
china <- global_economy |>
filter(Country == "China") |>
transmute(Year, GDP_tril = GDP/1e12)
# First test plot
autoplot(china, GDP_tril) + labs(title = "China GDP", y = "USD (trillions)")
# Experiment with multiple ETS options
# let ETS choose
# addtrend: additive (undamped) trend, no seasonality
# damped: additive damped trend (Ad)
# boxcox_damped: damped trend fitted on a Box–Cox scale
fits <- china |>
model(
auto = ETS(GDP_tril),
addtrend = ETS(GDP_tril ~ trend("A")),
damped = ETS(GDP_tril ~ trend("Ad")),
boxcox_damped = ETS(box_cox(GDP_tril, 0) ~ trend("Ad")) # log transform (λ=0)
)
# 30 year forecast to see significant differences
fc <- forecast(fits, h = 30)
autoplot(fc) +
facet_wrap(vars(.model), scales = "fixed") +
labs(title = "ETS forecasts for China GDP", y = "USD (trillions)")
8.7
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?
gas <- aus_production |> select(Quarter, Gas) |> filter(!is.na(Gas))
autoplot(gas, Gas) + labs(title = "Australian gas production")
# Two ETS models: additive trend with multiplicative seasonality,
# with and without damping of the trend
fit <- gas |>
model(
ANN_M = ETS(Gas ~ error("A") + trend("A") + season("M")),
AdNN_M = ETS(Gas ~ error("A") + trend("Ad") + season("M")) # damped trend
)
# Forecast next 4 years
fc <- fit |> forecast(h = "4 years")
autoplot(fc, gas) + labs(title = "ETS forecasts for Gas")
# Compare fit criteria
glance(fit) |> dplyr::select(.model, AICc, BIC)# In-sample error metrics
accuracy(fit) |> dplyr::select(.model, RMSE, MASE)dplyr::left_join(
glance(fit) |> dplyr::select(.model, AICc, BIC),
accuracy(fit) |> dplyr::select(.model, RMSE, MASE),
by = ".model"
)Multiplicative seasonality is needed because the seasonal swings get bigger as gas production rises, so therefore the size of the ups and downs scales with the level an additive model would miss that. Adding a damped trend would only flattens the longrun forecasts and doesn’t really improve accuracy here, so the undamped model is the best choice based on AICc/BIC and RMSE.
8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
a. Why is multiplicative seasonality necessary for this series?
The seasonal ups and downs get bigger as the series grows, so the effect scales with the level. Multiplicative seasonality captures that such as the peaks and dips grow with the trend, while additive keeps a fixed size instead.
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
myseries <- aus_retail |>
filter(State == "Northern Territory",
Industry == "Clothing, footwear and personal accessory retailing") |>
select(Month, Turnover)fits <- myseries |>
model(
HW_M = ETS(Turnover ~ error("A") + trend("A") + season("M")),
HW_M_damped = ETS(Turnover ~ error("A") + trend("Ad") + season("M"))
)
fc <- forecast(fits, h = "36 months")
autoplot(fc, myseries) +
facet_wrap(vars(.model), ncol = 1) +
labs(title = "Holt Winters multiplicative forecasts",
y = "Turnover")
c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
acc <- fit |> accuracy() |> dplyr::select(.model, RMSE)
accBased on the results, I prefer ANN_M as it has the lower one-step RMSE 4.19 compared to the 4.22 for the damped version, so it’s more accurate.
d. Check that the residuals from the best method look like white noise.
best_fit <- myseries |> model(HW_M_damped = ETS(Turnover ~ error("A") + trend("Ad") + season("M")))
best_fit |> gg_tsresiduals()## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
augment(best_fit) |> # Ljung–Box test for remaining autocorrelation
features(.innov, ljung_box, lag = 24, dof = 3)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?
myseries_train <- myseries |> filter(year(Month) <= 2010)
myseries_test <- anti_join(myseries, myseries_train)## Joining with `by = join_by(Month, Turnover)`fits_train <- myseries_train |>
model(
HW_M = ETS(Turnover ~ error("A") + trend("A") + season("M")),
HW_M_damped = ETS(Turnover ~ error("A") + trend("Ad") + season("M")),
SNAIVE = SNAIVE(Turnover ~ lag("year"))
)
fc_test <- forecast(fits_train, new_data = myseries_test)
acc_test <- accuracy(fc_test, myseries) |>
dplyr::select(.model, RMSE) |>
arrange(RMSE)
acc_testBased on the results, the damped Holt–Winters multiplicative model has lower RMSE (1.362746) than SNAIVE (1.551981), so it beats the seasonal naïve benchmark, while the undamped HW_M does not (1.771684).
8.9
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?
# train/tset
yseries_train <- myseries |> filter(year(Month) < 2011)
myseries_test <- anti_join(myseries, myseries_train)## Joining with `by = join_by(Month, Turnover)`lambda <- myseries_train |>
features(Turnover, guerrero) |>
dplyr::pull(lambda_guerrero)
# STL on Box–Cox data
fit_stl <- myseries_train |>
model(
STL_ETS = decomposition_model( # ETS on the seasonally-adjusted series
STL(box_cox(Turnover, lambda) ~ season(window = "periodic")),
ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
)
)
# test period and test RMSE
fc_stl <- fit_stl |> forecast(new_data = myseries_test)
acc_stl <- fc_stl |> accuracy(myseries) |> dplyr::select(.model, RMSE)
fit_hw <- myseries_train |>
model(HW_M_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fc_hw <- fit_hw |> forecast(new_data = myseries_test)
acc_hw <- fc_hw |> accuracy(myseries) |> dplyr::select(.model, RMSE)
dplyr::bind_rows(acc_hw, acc_stl)Based on these results, the performance of STL_ETS was not as accurate with a RMSE of 1.212853, compared to my previous best HW_M_damped’s RMSE of 1.151260.