Homework 5
Catherine Dube
2026-03-08
Exercises: 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman.
Exercise 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 \(\ell\) 0 not, and generate forecasts for the next four months.
data(aus_livestock)
aus_livestock %>%
filter(Animal == "Pigs" , State =="Victoria") %>%
autoplot(Count)+
labs(title = "Count of Pigs Slaughtered in Victoria over Time (No Modeling or Forecasting)") 
pigs_df <- aus_livestock %>%
filter(Animal == "Pigs" , State =="Victoria")
pigs_fit <- pigs_df %>%
model(ETS(Count ~ error("A") + trend("N") + season("N")))
# specifically want a simple exponential smoothing model (N,N), I think I'm supposed to specify the error as additive
report(pigs_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.07pigs_fit %>%
forecast(h = 4) %>%
autoplot(pigs_df) +
labs(title = "Count of Pigs Slaughtered in Victoria over Time (N,N) Model") 
The optimal values of \(\alpha\) and \(\ell\) 0 were given in the report.
- \(\alpha\) = 0.32
- \(\ell\) 0 = 100646
b
Compute a 95% prediction interval for the first forecast using \(\hat{y} \pm 1.96s\) where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
# Getting the residuals from the model
residuals(pigs_fit)## # A tsibble: 558 x 5 [1M]
## # Key: Animal, State, .model [1]
## Animal State .model Month .resid
## <fct> <fct> <chr> <mth> <dbl>
## 1 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Jul -6447.
## 2 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Aug 7130.
## 3 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Sep -4367.
## 4 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Oct 17640.
## 5 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Nov -542.
## 6 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1972 Dec -4468.
## 7 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1973 Jan -8829.
## 8 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1973 Feb -6785.
## 9 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1973 Mar -5299.
## 10 Pigs Victoria "ETS(Count ~ error(\"A\") + trend(\"N\") + … 1973 Apr -18792.
## # ℹ 548 more rows# Getting the sd of the residuals
pigs_sd = sd(residuals(pigs_fit)$.resid)
# Getting y_hat
# this is the mean of the forecasted values
pigs_y_hat = (pigs_fit %>% forecast(h = 4))$.mean[1]
# Compute the prediction using the y_hat +/ 1.96*sd
pigs_prediction_lower = pigs_y_hat - 1.96*pigs_sd
pigs_prediction_upper = pigs_y_hat + 1.96*pigs_sd
pigs_compute_prediction_interval = rbind(pigs_prediction_lower, pigs_prediction_upper)
# Getting the prediction interval produced by R
pigs_fit %>%
forecast(h = 4) %>%
hilo(95) %>%
pull('95%') %>%
head(1)## <hilo[1]>
## [1] [76854.79, 113518.3]95The manually calculated prediction and the R calculated one are nearly the same
(768XX, 1135XX)
Exercise 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.
data(global_economy)
global_economy %>%
filter(Country == "Vietnam") %>%
autoplot(Exports) ## Warning: Removed 26 rows containing missing values or values outside the scale range
## (`geom_line()`).
I selected Vietnam to look at. It appears that exports from Vietnam were not tracked until the mid-late 1980s. There was a noticeable spike in exports in 1990, and since then the exports have been steadily increasing. There was also a noticeable dip in 2009 which recovered quickly.
b
Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
viet_exports <- global_economy %>%
filter(Country == "Vietnam") %>%
drop_na(Exports) ## I had to add this line because the missing data at the start was causing errors
viet_model <- viet_exports %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
viet_fc <- viet_model %>%
forecast(h = 6)
viet_fc %>%
autoplot(viet_exports)
c
Compute the RMSE values for the training data.
accuracy(viet_model)$RMSE## [1] 5.964381Rounding this RMSE up to 6, that means this model is off by about 6 units (percentage points of GDP).
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.
viet_model_AAN <- viet_exports %>%
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
viet_fc_AAN <- viet_model_AAN %>%
forecast(h = 6)
viet_fc_AAN %>%
autoplot(viet_exports)
accuracy(viet_model_AAN)$RMSE## [1] 5.120112The RMSE value of the AAN model is lower, suggesting it is a better fit. And speaking just visually speaking, the prediction looks a lot closer to what I would expect.
e
Compare the forecasts from both methods. Which do you think is best?
Including a trend component seems particularly important in order to have a clearer expectation of what is to come. I think the AAN model is a better fit.
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.
# Getting the sd of the residuals
viet_ANN_sd <- sd(residuals(viet_model)$.resid)
viet_AAN_sd <- sd(residuals(viet_model_AAN)$.resid)
# Getting y_hat
# this is the mean of the forecasted values
viet_ANN_y_hat <- viet_fc$.mean[1]
viet_AAN_y_hat <- viet_fc_AAN$.mean[1]
# Compute the prediction using the y_hat +/ 1.96*sd
viet_ANN_lower = viet_ANN_y_hat - 1.96*viet_ANN_sd
viet_ANN_upper = viet_ANN_y_hat + 1.96*viet_ANN_sd
viet_AAN_lower = viet_AAN_y_hat - 1.96*viet_AAN_sd
viet_AAN_upper = viet_AAN_y_hat + 1.96*viet_AAN_sd
viet_ANN_computed_interval = rbind(viet_ANN_lower, viet_ANN_upper)
viet_AAN_computed_interval = rbind(viet_AAN_lower, viet_AAN_upper)
viet_ANN_computed_interval## [,1]
## viet_ANN_lower 91.28988
## viet_ANN_upper 111.89540viet_AAN_computed_interval## [,1]
## viet_AAN_lower 93.55566
## viet_AAN_upper 113.94752# Getting the prediction interval produced by R
viet_ANN_interval_from_R <-viet_fc %>%
hilo(95) %>%
pull('95%') %>%
head(1)
viet_AAN_interval_from_R <- viet_fc_AAN %>%
hilo(95) %>%
pull('95%') %>%
head(1)
viet_ANN_interval_from_R## <hilo[1]>
## [1] [89.51929, 113.666]95viet_AAN_interval_from_R## <hilo[1]>
## [1] [93.02347, 114.4797]95These intervals seem a little less close together potentially than the ones from the previous homework problem, but that’s probably just because the magnitude of the numbers is smaller.
Exercise 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_df <- global_economy %>%
filter(Country == "China")
lambda <- china_df %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
china_fit <- global_economy %>%
filter(Country == "China") %>%
model(`Simple (A,N,N)` = ETS(GDP ~ error("A") + trend("N") + season("N")),
`Holt's method` = ETS(GDP ~ error("A") + trend("A") + season("N")),
`Damped (0.8) Holt's method` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Damped (0.9) Holt's method` = ETS(GDP ~ error("A") + trend("Ad", phi = 0.9) + season("N")),
`Box-Cox` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("A") + season("N")),
`Box-Cox Damped (0.8)` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Box-Cox Damped (0.9)` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad", phi = 0.9) + season("N")),
`Log` = ETS(log(GDP) ~ error("A") + trend("A") + season("N")),
`Log Damped (0.8)` = ETS(log(GDP) ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Log Damped (0.9)` = ETS(log(GDP) ~ error("A") + trend("Ad", phi = 0.9) + season("N"))
)
china_fc <- china_fit %>%
forecast(h = 20)
china_fc %>%
autoplot(china_df, level = NULL) +
labs(title="Forecast of China GDP With Multiple Methods") +
guides(colour = guide_legend(title = "Forecast Method"))
This is interesting to experiment with. The most obvious takeaway to me is how important damping is. In particular, Box-Cox Damped vs. Box-Cox is the difference between reality and something way beyond optimism. Aside from that, I also tried a couple phi values for dampening. The smaller (0.8) value pulls the prediction level much faster than the higher (0.9) value. I didn’t realize that both undamped Box-Cox and Log could be so unruly.
Exercise 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?
data(aus_production)
gas_fit <- aus_production %>%
select(Gas) %>%
model(`Simple (A,N,N)` = ETS(Gas ~ error("A") + trend("N") + season("N")),
`Additive` = ETS(Gas ~ error("A") + trend("A") + season("A")),
`Multiplicative` = ETS(Gas ~ error("M") + trend("A") + season("M")),
`Damped Additive` = ETS(Gas ~ error("M") + trend("Ad", phi = 0.85) + season("A")),
`Damped Multiplicative` = ETS(Gas ~ error("M") + trend("Ad", phi = 0.85) + season("M"))
)
gas_fit %>%
forecast(h=10) %>%
autoplot(aus_production, level = NULL)
report(gas_fit)## Warning in report.mdl_df(gas_fit): 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: 5 × 9
## .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Simple (A,N,N) 280. -1200. 2406. 2406. 2416. 277. 273. 12.6
## 2 Additive 23.6 -927. 1872. 1873. 1903. 22.7 29.7 3.35
## 3 Multiplicative 0.00324 -831. 1681. 1682. 1711. 21.1 32.2 0.0413
## 4 Damped Additive 0.0764 -1175. 2369. 2370. 2399. 29.7 37.3 0.0954
## 5 Damped Multiplicative 0.00347 -837. 1692. 1693. 1722. 20.9 32.5 0.0429Multiplicative is necessary here because the shape of the data is like a funnel widening. Because that amplitude is growing multiplicative seasonality is needed to account for its growth. Making the trend damped keeps the amplitude a smidge smaller, but it hard to tell whether or not it is better without checking the report.
Looking at the report, the damped multiplicative model did have a lower error than the multiplicative model.
Exercise 8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
myseries <- aus_retail %>%
filter(`Series ID` == "A3349476W")
myseries %>% autoplot()## Plot variable not specified, automatically selected `.vars = Turnover`
a
Why is multiplicative seasonality necessary for this series?
Multiplicative seasonality is necessary for this series because the data shows a trend, shows seasonality, and shows a growing funnel shape. An additive seasonal component would sum to approximately zero, which would not accurately represent this growth.
b
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
my_fit <- myseries %>%
model(`Simple (A,N,N)` = ETS(Turnover ~ error("A") + trend("N") + season("N")),
`Holt-Winters Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Holt-Winters Multiplicative & Damped (0.8)` = ETS(Turnover ~ error("M") + trend("Ad", phi = 0.8) + season("M")),
`Holt-Winters Multiplicative & Damped (0.9)` = ETS(Turnover ~ error("M") + trend("Ad", phi = 0.9) + season("M"))
)
my_fit %>%
forecast(h=45) %>%
autoplot(myseries, level = NULL)
c
report(my_fit)## Warning in report.mdl_df(my_fit): 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 × 11
## State Industry .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Queensla… Pharmac… Simpl… 2.14e+2 -2525. 5056. 5056. 5068. 213. 250. 9.43
## 2 Queensla… Pharmac… Holt-… 3.31e-3 -2041. 4115. 4117. 4185. 50.1 72.9 0.0433
## 3 Queensla… Pharmac… Holt-… 3.72e-3 -2063. 4161. 4162. 4230. 55.3 80.3 0.0463
## 4 Queensla… Pharmac… Holt-… 3.29e-3 -2037. 4108. 4110. 4178. 48.4 71.4 0.0428Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
It looks like the damped model at phi = 0.9 performed the best. Visually inspecting, all three non simple methods were actually quite similar in performance for this series of data though.
d
Check that the residuals from the best method look like white noise.
residuals(my_fit)## # A tsibble: 1,764 x 5 [1M]
## # Key: State, Industry, .model [4]
## State Industry .model Month .resid
## <chr> <chr> <chr> <mth> <dbl>
## 1 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Apr -0.270
## 2 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 May -0.261
## 3 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Jun 0.644
## 4 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Jul 1.98
## 5 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Aug 0.484
## 6 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Sep -1.51
## 7 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Oct -3.40
## 8 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Nov -2.13
## 9 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1982 Dec 3.33
## 10 Queensland Pharmaceutical, cosmetic and toiletry good… Simpl… 1983 Jan -4.11
## # ℹ 1,754 more rowsmyseries %>%
model(`Holt-Winters Multiplicative & Damped (0.9)` = ETS(Turnover ~ error("M") + trend("Ad", phi = 0.9) + season("M"))) %>%
gg_tsresiduals() +
ggtitle("Multiplicative Method and Damped 0.9")## 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.
The residual plot does appear to be white noise.However, the mean of the residuals is actually a bit left of 0, and there appears to potentially be some pattern in the acf plot.
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) < 2011)
fit_train <- myseries_train %>%
model(`Holt-Winters Multiplicative & Damped (0.9)` = ETS(Turnover ~ error("M") + trend("Ad", phi = 0.9) + season("M")),
snaive = SNAIVE(Turnover))
# Forecasts
fc <- fit_train %>%
forecast(new_data = anti_join(myseries, myseries_train))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc %>% autoplot(myseries, level = NULL)
accuracy(fit_train) %>%
select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Training Holt-Winters Multiplicative & Damped (0.9) 5.36
## 2 Training snaive 14.6fc %>% accuracy(myseries) %>%
select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test Holt-Winters Multiplicative & Damped (0.9) 77.6
## 2 Test snaive 44.3This is really interesting. Visually I think the Snaive model performs better, because it is predicting a higher value overall which is closer the the actual. However, the RMSE reveals that the Snaive model has a worse RMSE.
Exercise 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 <- myseries %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
# STL Breakout of the data
myseries %>%
model(STL(box_cox(Turnover,lambda_myseries) ~ season(window = "periodic"), robust = TRUE)) %>%
components() %>%
autoplot() +
ggtitle("STL of My Series (Retail Data)")
# the models
fit_train_with_STL <- myseries_train %>%
model(
#STL
`STL_Forecast` = decomposition_model(
STL(box_cox(Turnover, lambda_myseries) ~ season(window = "periodic"), robust = TRUE),
RW(season_adjust ~ drift()) # Forecasts the trend/remainder using a Random Walk with drift
),
`Holt-Winters Multiplicative & Damped (0.9)` = ETS(Turnover ~ error("M") + trend("Ad", phi = 0.9) + season("M")),
snaive = SNAIVE(Turnover)
)
# Forecasts
fc <- fit_train_with_STL %>%
forecast(new_data = anti_join(myseries, myseries_train))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc %>% autoplot(myseries, level = NULL)
accuracy(fit_train) %>%
select(.type, .model, RMSE)## # A tibble: 2 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Training Holt-Winters Multiplicative & Damped (0.9) 5.36
## 2 Training snaive 14.6fc %>% accuracy(myseries) %>%
select(.type, .model, RMSE)## # A tibble: 3 × 3
## .type .model RMSE
## <chr> <chr> <dbl>
## 1 Test Holt-Winters Multiplicative & Damped (0.9) 77.6
## 2 Test STL_Forecast 88.8
## 3 Test snaive 44.3# STL Residuals
fit_train_with_STL %>%
select(STL_Forecast) %>%
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()`).
The STL forecast continues overshooting the Turnover versus the other methods ( at least visually). It also has the highest RMSE value. Clearly the ETS method is a little more robust.