# Loading required packages
library(fpp3)DATA 624 Homework 5
8.1
Consider the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
pigs <- aus_livestock |>
filter(Animal == "Pigs", State == "Victoria")
pigs |> autoplot(Count)
The series is monthly with no clear seasonal pattern and a slowly changing level, so simple exponential smoothing is a reasonable starting point.
a. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of alpha and l0, and generate forecasts for the next four months.
fit <- pigs |>
model(ETS(Count ~ error("A") + trend("N") + season("N")))
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 Alpha is 0.3221, so each new month gets ~ 32% of the weight and the older history shares the other 68%, which means the level reacts at a moderate speed. The starting level l0 is 100,646.6 pigs and sigma squared is 87,480,760.
fc <- fit |> forecast(h = 4)
fc |> autoplot(pigs |> filter(year(Month) >= 2010))
Simple exponential smoothing has no trend or season, so all four forecasts sit at the same flat value of ~ 90000 pigs.
b. Compute a 95% prediction interval for the first forecast using y hat plus or minus 1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
s <- augment(fit) |> pull(.resid) |> sd(na.rm = TRUE)
yhat <- fc$.mean[1]
c(lower = yhat - 1.96 * s, upper = yhat + 1.96 * s) lower upper
76871.01 113502.10 fc |> hilo(level = 95) |> pull(`95%`) |> head(1)<hilo[1]>
[1] [76854.79, 113518.3]95My interval is (76,871, 113,502) and R gives (76,855, 113,518). R is wider by ~16 on each side, under 0.1% of the interval, because it uses the model’s own sigma (the square root of 87,480,760, ~ 9350) instead of the sample standard deviation of the residuals.
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.
aus <- global_economy |> filter(Country == "Australia")
aus |> autoplot(Exports)
Australian exports are yearly so there is no seasonality. The series moves up slowly over time with short ups and downs around that movement.
b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
fit_ann <- aus |>
model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
fit_ann |>
forecast(h = 5) |>
autoplot(aus)
The ANN forecast is a flat line at the last smoothed level.
c. Compute the RMSE values for the training data.
accuracy(fit_ann) |> select(.model, RMSE)# A tibble: 1 × 2
.model RMSE
<chr> <dbl>
1 ANN 1.15The training RMSE for ANN is 1.1468.
d. 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.
fit <- aus |>
model(
ANN = ETS(Exports ~ error("A") + trend("N") + season("N")),
AAN = ETS(Exports ~ error("A") + trend("A") + season("N"))
)
accuracy(fit) |> select(.model, RMSE)# A tibble: 2 × 2
.model RMSE
<chr> <dbl>
1 ANN 1.15
2 AAN 1.12glance(fit) |> select(.model, AICc)# A tibble: 2 × 2
.model AICc
<chr> <dbl>
1 ANN 258.
2 AAN 259.AAN fits the training data slightly better, with an RMSE of 1.1167 vs 1.1468 for ANN, a drop of 0.0301. But its AICc is 259.47 vs 257.84 for ANN, so the extra trend parameters cost more than they gain. ANN gives the better balance of fit and simplicity for this series.
e. Compare the forecasts from both methods. Which do you think is best?
fit |>
forecast(h = 5) |>
autoplot(aus, level = NULL)
The ANN line stays flat while the AAN line keeps the recent upward slope going, and the first forecast is 20.84 for AAN vs 20.61 for ANN. Because ANN has the lower AICc (257.84 versus 259.47), I think ANN is the best choice.
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 <- fit |> forecast(h = 5)
rmse <- accuracy(fit) |> select(.model, RMSE)
fc |>
as_tibble() |>
group_by(.model) |>
slice(1) |>
ungroup() |>
select(.model, .mean) |>
left_join(rmse, by = ".model") |>
mutate(lower = .mean - 1.96 * RMSE, upper = .mean + 1.96 * RMSE)# A tibble: 2 × 5
.model .mean RMSE lower upper
<chr> <dbl> <dbl> <dbl> <dbl>
1 AAN 20.8 1.12 18.6 23.0
2 ANN 20.6 1.15 18.4 22.9fc |>
hilo(level = 95) |>
as_tibble() |>
group_by(.model) |>
slice(1) |>
ungroup() |>
select(.model, `95%`)# A tibble: 2 × 2
.model `95%`
<chr> <hilo>
1 AAN [18.57028, 23.10700]95
2 ANN [18.31970, 22.89462]95The manual intervals are a little narrower than the ones from R. For ANN mine is (18.36, 22.85) and R gives (18.32, 22.89). For AAN mine is (18.65, 23.03) and R gives (18.57, 23.11). RMSE divides the squared errors by n, while Rs sigma divides by n minus the number of estimated parameters, which is 2 for ANN (alpha and l0) and 4 for AAN (alpha, beta, l0 and b0). That is why the gap is larger for AAN.
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.
china <- global_economy |> filter(Country == "China")
china |> autoplot(GDP)
China’s GDP grows faster and faster over time, so the growth looks exponential rather than linear. That is a good reason to try a Box Cox transformation.
lambda <- china |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda[1] -0.03446284The Guerrero lambda is ~ -0.034, very close to 0, so the Box Cox transformation is basically a log. On the log scale exponential growth becomes a straight line, which is the kind of trend ETS handles well.
fit <- china |>
model(
ETS = ETS(GDP),
Damped = ETS(GDP ~ trend("Ad")),
BoxCox = ETS(box_cox(GDP, lambda)),
BoxCoxDamped = ETS(box_cox(GDP, lambda) ~ trend("Ad"))
)
glance(fit) |> select(.model, AICc)# A tibble: 4 × 2
.model AICc
<chr> <dbl>
1 ETS 3103.
2 Damped 3107.
3 BoxCox -137.
4 BoxCoxDamped -133.AICc cannot be compared between the two groups, because the Box Cox models are measured on the transformed scale (~ -137 and -133) and the others on the original GDP scale (~ 3,103 and 3,107). Within each group the undamped model wins: ETS 3,103.2 vs Damped 3,106.7 (a gap of 3.5), and BoxCox -137.4 vs BoxCoxDamped -133.2 (a gap of 4.1). So damping does not improve the fit for this series.
fit |>
forecast(h = 30) |>
autoplot(china, level = NULL)
The damped trend models flatten out and end up below the undamped ones at long horizons. The Box Cox models grow in the transformed scale, so after back transforming they keep climbing faster and their intervals become uneven.
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?
aus_production |> autoplot(Gas)
The seasonal swings get larger as the level of gas production rises. A constant size seasonal pattern cannot capture that, so multiplicative seasonality is needed.
fit <- aus_production |>
model(
Auto = ETS(Gas),
HW = ETS(Gas ~ error("M") + trend("A") + season("M")),
HW_damped = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)
report(fit |> select(Auto))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 The automatic search chose ETS(M,A,M), which is multiplicative error, additive trend and multiplicative seasonality, matching what the plot shows. Alpha is 0.653, beta is 0.144 and gamma is 0.098, so the level reacts quickly while the trend and the seasonal pattern change slowly.
glance(fit) |> select(.model, AICc)# A tibble: 3 × 2
.model AICc
<chr> <dbl>
1 Auto 1682.
2 HW 1682.
3 HW_damped 1685.accuracy(fit) |> select(.model, RMSE, MAE, MAPE)# A tibble: 3 × 4
.model RMSE MAE MAPE
<chr> <dbl> <dbl> <dbl>
1 Auto 4.60 3.02 4.08
2 HW 4.60 3.02 4.08
3 HW_damped 4.59 3.03 4.10Auto and HW give identical results because the automatic search picked the same ETS(M,A,M). The damped version has a higher AICc (1685.09 vs 1681.79) nearly identical RMSE (4.5918 vs 4.5951), and a slightly higher MAE (3.031 vs 3.022) and MAPE (4.10 vs 4.08).
fit |>
forecast(h = "5 years") |>
autoplot(aus_production |> filter(year(Quarter) >= 2000), level = NULL)
Damping does not improve the forecasts here. Its RMSE is lower by only 0.0033 (4.5951 - 4.5918), while its MAE, MAPE and AICc all get worse. The damped forecast grows more slowly, but the extra phi parameter is not worth it, so I keep the undamped ETS(M,A,M).
8.8
Recall your retail time series data
set.seed(12345678)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
myseries |> autoplot(Turnover)
a. Why is multiplicative seasonality necessary for this series?
The seasonal swings grow together with the level of turnover, so the seasonal effect is a percentage of the level and not a fixed amount. This is the case where multiplicative seasonality fits. The fitted starting seasonal factors range from ~ 0.76 to 1.37 of the level, which shows the seasonal effect is a percentage of the level and not a fixed amount.
b. Apply Holt Winters multiplicative method to the data. Experiment with making the trend damped.
fit <- myseries |>
model(
HW = ETS(Turnover ~ error("M") + trend("A") + season("M")),
HW_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)
report(fit |> select(HW))Series: Turnover
Model: ETS(M,A,M)
Smoothing parameters:
alpha = 0.5262368
beta = 0.0001026197
gamma = 0.0001037781
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
2.347713 0.04437301 0.8398709 0.7617432 0.833525 1.373239 0.9998723 1.033418
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
1.032927 1.119514 1.140226 1.019343 0.9810424 0.8652791
sigma^2: 0.0045
AIC AICc BIC
1774.530 1776.274 1841.014 report(fit |> select(HW_damped))Series: Turnover
Model: ETS(M,Ad,M)
Smoothing parameters:
alpha = 0.5566535
beta = 0.0002378027
gamma = 0.1112567
phi = 0.9794538
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
2.527185 0.109247 0.8093026 0.746469 0.8360684 1.354604 1.003636 1.045657
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
1.034894 1.093245 1.160544 1.054001 0.9937055 0.8678735
sigma^2: 0.0046
AIC AICc BIC
1780.720 1782.674 1851.114 In the HW model alpha is 0.526 while beta (0.0001) and gamma (0.0001) are almost 0, so the trend and the seasonal pattern stay nearly fixed and only the level updates. In the damped model phi is 0.979, so the trend shrinks by ~2% at each step.
fit |>
forecast(h = "3 years") |>
autoplot(myseries |> filter(year(Month) >= 2010), level = NULL)
c. Compare the RMSE of the one step forecasts from the two methods. Which do you prefer?
accuracy(fit) |> select(.model, RMSE, MAE, MAPE)# A tibble: 2 × 4
.model RMSE MAE MAPE
<chr> <dbl> <dbl> <dbl>
1 HW 0.613 0.450 5.15
2 HW_damped 0.616 0.444 5.06glance(fit) |> select(.model, AICc)# A tibble: 2 × 2
.model AICc
<chr> <dbl>
1 HW 1776.
2 HW_damped 1783.The one step RMSE is almost the same, 0.6130 for HW and 0.6155 for HW_damped, a gap of 0.0025. HW also has the lower AICc (1776.27 vs 1782.67), while the damped model is marginally better on MAE (0.4444 vs 0.4497) and MAPE (5.06 vs 5.15). Since the question is about RMSE and AICc agrees, I prefer HW.
d. Check that the residuals from the best method look like white noise.
best <- accuracy(fit) |> slice_min(RMSE) |> pull(.model)
fit |> select(all_of(best)) |> gg_tsresiduals()
augment(fit) |>
filter(.model == best) |>
features(.innov, ljung_box, lag = 24)# A tibble: 1 × 5
State Industry .model lb_stat lb_pvalue
<chr> <chr> <chr> <dbl> <dbl>
1 Northern Territory Clothing, footwear and personal a… HW 34.5 0.0765The residuals are centered on 0 with no trend and a fairly steady spread, and the histogram is close to normal with a slight right tail. In the ACF almost every lag sits inside the blue bounds. Only lag 12 ( ~0.14) clearly crosses the bound, while lags 8 and 18 just touch it at ~ -0.11 and -0.10. With 24 lags we expect ~ 1 to cross by chance (24 × 0.05 = 1.2), so the residuals are close to white noise. The spike at lag 12 points to a small yearly pattern the model has not fully captured.
e. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naive approach from Exercise 7 in Section 5.11?
myseries_train <- myseries |>
filter(year(Month) < 2011)
fit_train <- myseries_train |>
model(
SNAIVE = SNAIVE(Turnover),
HW = ETS(Turnover ~ error("M") + trend("A") + season("M")),
HW_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)
fc_train <- fit_train |>
forecast(new_data = anti_join(myseries, myseries_train))Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc_train |>
accuracy(myseries) |>
select(.model, RMSE, MAE, MAPE)# A tibble: 3 × 4
.model RMSE MAE MAPE
<chr> <dbl> <dbl> <dbl>
1 HW 1.78 1.60 12.6
2 HW_damped 1.15 0.878 6.45
3 SNAIVE 1.55 1.24 9.06On the test set HW_damped is the clear winner with an RMSE of 1.151, vs 1.552 for seasonal naive and 1.783 for HW. HW_damped beats seasonal naive by 0.401 (1.552 - 1.151), ~ 26% lower, so yes it can be beaten. The undamped HW is worse than seasonal naive by 0.231 (1.783 - 1.552). Its additive trend keeps the full slope through the whole test period, while the damped trend flattens, which suggests the series grew less after 2010 than the undamped trend assumed. It also shows that the better model in sample (HW) was not the better forecaster.
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?
lambda <- myseries_train |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
lambda[1] 0.1971431Lambda is ~ 0.197, between a log (0) and a square root (0.5) transformation.
fit_stl <- myseries_train |>
model(
STL_ETS = decomposition_model(
STL(box_cox(Turnover, lambda) ~ season(window = "periodic"),
robust = TRUE),
ETS(season_adjust ~ season("N"))
)
)
fit_stl |> gg_tsresiduals()Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_point()`).Warning: Removed 12 rows containing non-finite outside the scale range
(`stat_bin()`).Warning: Removed 12 rows containing missing values or values outside the scale range
(`geom_rug()`).
fc_stl <- fit_stl |>
forecast(new_data = anti_join(myseries, myseries_train))Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc_stl |> autoplot(myseries |> filter(year(Month) >= 2008))
bind_rows(fc_train, fc_stl) |>
accuracy(myseries) |>
select(.model, RMSE, MAE, MAPE) |>
arrange(RMSE)# A tibble: 4 × 4
.model RMSE MAE MAPE
<chr> <dbl> <dbl> <dbl>
1 HW_damped 1.15 0.878 6.45
2 SNAIVE 1.55 1.24 9.06
3 HW 1.78 1.60 12.6
4 STL_ETS 4.24 3.81 29.4 The Box Cox step handles the growing seasonal swings, STL pulls the seasonal part out, and ETS forecasts only the seasonally adjusted data. Even so, STL with ETS does much worse than the earlier models, with a test RMSE of 4.236 vs 1.151 for HW_damped, 1.552 for seasonal naive and 1.783 for HW. That is ~ 3.7 times the error of the best model (4.236 /1.151 = 3.68), and the MAPE is 29.4% vs 6.4%. One likely reason is the back transformation: with lambda at 0.197 the original scale is the transformed value raised to about the power 1 / 0.197 = 5.1, so even a small trend in the ETS forecast gets magnified when it is converted back to turnover. The best forecasts are still the HW_damped ones.