DATA624_Assignment5
2025-10-04
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 L0, and generate forecasts for the next four months.
’
fit <- aus_livestock %>%
filter(State == "Victoria") %>%
filter(Animal == "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.07fc <- fit %>%
forecast(h = 5)
fc %>%
autoplot(aus_livestock)
Alpha = 0.3221247 and the initial L0 is 100646.6
b ) Compute a 95% prediction interval for the first forecast using ^y ± 1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
resid <- augment(fit) %>% select(.resid)
s <- sd(resid$.resid, na.rm = TRUE)
yhat <- fc$.mean[1]
manual_lower <- yhat - 1.96 * s
manual_upper <- yhat + 1.96 * s
fc_intervals <- fc %>%
hilo(level = 95) %>%
mutate(
r_lower = `95%`$lower,
r_upper = `95%`$upper
)
comparison <- fc_intervals %>%
slice(1) %>%
select(.mean, r_lower, r_upper) %>%
mutate(
manual_lower = manual_lower,
manual_upper = manual_upper
) %>%
select(.mean, manual_lower, manual_upper, r_lower, r_upper)
comparison## # A tsibble: 1 x 6 [1M]
## .mean manual_lower manual_upper r_lower r_upper Month
## <dbl> <dbl> <dbl> <dbl> <dbl> <mth>
## 1 95187. 76871. 113502. 76855. 113518. 2019 JanThe confidence interval ranges are almost identical, with the R estimated range being slightly wider at both ends
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.
global_economy %>%
filter(Country == "India") %>%
autoplot(Exports)
The exports for India show a generally upward trend starting at 5
percent at the beginning of 1960s with a small dip during the 1980s
followed by exponential growth all the way up to 25 percent in the 2010s
followed by a sharp fall rise and another fall
b) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
india_exports <- global_economy %>%
filter(Country == "India") %>%
select(Year, Exports)
fit <- india_exports %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h = 5)
fc %>%
autoplot(india_exports)
c) Compute the RMSE values for the training data.
accuracy(fit)## # A tibble: 1 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "ETS(Exports ~ error(… Trai… 0.251 1.19 0.798 2.05 7.25 0.983 0.991 -0.0331d) 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_trend <- india_exports %>%
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
accuracy(fit_trend)## # A tibble: 1 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "ETS(Exports ~ erro… Trai… -0.0501 1.18 0.736 0.239 6.98 0.906 0.978 -0.0294The RMSE is lower for the model with the trend component. Simple Exponential smoothing is more stable for short-term forecasting more centered around around the recent level. While ETS(A,A,N) adds one parameter for the trend and results in a lower RMSE and therefore a better fit, and is better suited fro a longer term forecast especially when we can see a clear trend in the data .
e) Compare the forecasts from both methods. Which do you think is best?
fc_trend <- fit_trend %>%
forecast(h = 5)
fc_trend %>%
autoplot(india_exports)
Here the trend related smoothing equation helps the model actually understand the recently decreasing characteristics of the exports and bring that into the forecasting, while the simple exponential smoothing models forecasts followed a more stable line across the next 5 years forecasts.
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.
yhat_simple<- fc$.mean[1]
yhat_trend <- fc_trend$.mean[1]
rmse_simple <- accuracy(fit) %>% pull(RMSE)
rmse_trend <- accuracy(fit_trend) %>% pull(RMSE)
manual_simple <- c(
lower = yhat_simple - 1.96 * rmse_simple,
upper = yhat_simple + 1.96 * rmse_simple
)
manual_trend <- c(
lower = yhat_trend - 1.96 * rmse_trend,
upper = yhat_trend + 1.96 * rmse_trend
)
r_simple <- fc %>%
slice(1) %>%
hilo(level = 95) %>%
mutate(
r_lower = `95%`$lower,
r_upper = `95%`$upper
)
r_trend <- fc_trend %>%
slice(1) %>%
hilo(level = 95) %>%
mutate(
r_lower = `95%`$lower,
r_upper = `95%`$upper
)
tibble(
Model = c("ETS(A,N,N)", "ETS(A,A,N)"),
Forecast = c(yhat_simple, yhat_trend),
Manual_Lower = c(manual_simple["lower"], manual_trend["lower"]),
Manual_Upper = c(manual_simple["upper"], manual_trend["upper"]),
R_Lower = c(r_simple$r_lower, r_trend$r_lower),
R_Upper = c(r_simple$r_upper, r_trend$r_upper)
)## # A tibble: 2 × 6
## Model Forecast Manual_Lower Manual_Upper R_Lower R_Upper
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS(A,N,N) 19.0 16.7 21.4 16.7 21.4
## 2 ETS(A,A,N) 18.2 15.9 20.5 15.8 20.68.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.
china_gdp <- global_economy %>%
filter(Country == "China") %>%
select(Year, GDP)
lambda <- china_gdp %>% features(GDP, features = guerrero) %>% pull(lambda_guerrero)
china_fit <- china_gdp %>%
model(
`ANN` = ETS(GDP ~ error("A") + trend("N") + season("N")),
`AAN ` = ETS(GDP ~ error("A") + trend("A") + season("N")),
`AMN ` = ETS(GDP ~ error("A") + trend("A") + season("N")),
`AAdN` = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
`Box-Cox` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("A") + season("N")),
`Damped Box-Cox` = ETS(box_cox(GDP,lambda) ~ error("A") + trend("Ad") + season("N"))
)
china_fc <- china_fit %>% forecast(h = 15)
china_fc %>%
autoplot(china_gdp, level = NULL) 
accuracy(china_fit)## # A tibble: 6 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 "ANN" Trai… 2.10e11 4.16e11 2.13e11 8.14 11.0 0.983 0.991 0.789
## 2 "AAN " Trai… 2.36e10 1.90e11 9.59e10 1.41 7.62 0.442 0.453 0.00905
## 3 "AMN " Trai… 2.36e10 1.90e11 9.59e10 1.41 7.62 0.442 0.453 0.00905
## 4 "AAdN" Trai… 2.95e10 1.90e11 9.49e10 1.62 7.62 0.438 0.454 -0.00187
## 5 "Box-Cox" Trai… -2.75e10 2.88e11 1.25e11 0.607 7.17 0.578 0.688 0.665
## 6 "Damped Box-C… Trai… 6.35e 9 1.96e11 1.02e11 1.77 7.26 0.472 0.468 0.135We see that using holts method to incorporate trend based smoothing into the forecasting helps the model a lot more when compared to the simple exponential smoothing method of ANN, with AAN and AMN giving similar RMSE values. Damped versions AAdN provide a little more realistic version when it comes to producing better long term forecasts if the growth is expected to even out .
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) 
Here we see that clear seasonal fluctuations along with a clear trend
pattern that shows seasonal variations that increase over time. This
points towards the necessity of a multiplicative seasonality
pattern.
Gas_fit <- aus_production %>%
model(
`ETS(A,A,A)` = ETS(Gas ~ error("A") + trend("A") + season("A")), # additive trend & additive seasonality
`ETS(A,A,M)` = ETS(Gas ~ error("A") + trend("A") + season("M")), # additive trend & multiplicative seasonality
`Damped ETS(A,Ad,M)` = ETS(Gas ~ error("A") + trend("Ad") + season("M")) # additive damped trend & multiplicative seasonality
)
Gas_fc <- Gas_fit %>% forecast(h = "8 years")
Gas_fc %>%
autoplot(aus_production, level = NULL) 
We see that the multiplicative forecast allows for increasing
fluctuations in seasonality within the predictions, while the additive
forecast retains the same variability. The damped forecasts look to be a
slightly better fit when considering long term forecasts as it reflects
how the trend gradually slows in different quarters and is not always at
a consistent rate of increase
8.8 Recall your retail time series data (from Exercise 7 in Section 2.10).
a) Why is multiplicative seasonality necessary for this series?
set.seed(90483484)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries %>%
autoplot(Turnover)
Similar to the previous example here the seasonal variance seems to
increase as the level of the series grows so we need to use
multiplicative seasonality that better scales with the level and is
better at capturing the seasonal peaks and lows that widen over
time.
b) Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
fit_myseries <- myseries %>%
model(`HW_Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`HW_Damped_Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fit_myseries %>%
forecast(h=50) %>%
autoplot(myseries, level = NULL) +
facet_wrap(~ .model)
c) Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
accuracy(fit_myseries)## # A tibble: 2 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 South Austr… Superma… HW_Mu… Trai… 0.203 7.63 5.85 -0.0762 2.23 0.334 0.350
## 2 South Austr… Superma… HW_Da… Trai… 1.06 7.87 6.09 0.279 2.28 0.348 0.361
## # ℹ 1 more variable: ACF1 <dbl>Here for one-step forecasts the Multiplicative model seems to have a lower RSME than the Damped multiplicative model showing us that damping maybe unnecessary. The damping reduces the trend and affects short term forecast accuracy
d) Check that the residuals from the best method look like white noise.
fit_myseries %>%
select(HW_Multiplicative) %>%
gg_tsresiduals()
Residuals seems almost normally distributed with not skewing and not
specific patterns. While they seem to oscillate around 0 there seems to
be a slightly higher variations in the earlier periods. The acf shows
peaks and bottoms random but crossing the significance intervals meaning
that the residuals are not fully uncorrelated and it doesn’t fully
satisfy the conditions for it being white noise
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_myseries_train <- myseries_train %>%
model(`HW_Multiplicative` = ETS(Turnover ~ error("M") + trend("M") + season("M")),
`HW_Damped_Multiplicative` = ETS(Turnover ~ error("M") + trend("Md") + season("M")),
snaive = SNAIVE(Turnover))
fc_myseries_train <- fit_myseries_train %>%
forecast(new_data = anti_join(myseries, myseries_train))
fc_myseries_train %>% autoplot(myseries, level = NULL)
fc_myseries_train %>% accuracy(myseries) %>% select(.model, RMSE)## # A tibble: 3 × 2
## .model RMSE
## <chr> <dbl>
## 1 HW_Damped_Multiplicative 109.
## 2 HW_Multiplicative 68.5
## 3 snaive 124.The HW_Multiplicative method was able to beat the Snaive here.
8.9 or 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_retail <- myseries %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
stl_fit <- myseries %>%
model(STL(box_cox(Turnover, lambda_retail) ~ season(window = "periodic"), robust = TRUE))
Decomp <- components(stl_fit)
myseries_adj <- myseries %>%
mutate(Turnover_sa = Decomp$season_adjust)
fit <- myseries_train %>%
model(ETS(Turnover ~ error("M") + trend("A") + season("N")))
fc <- fit %>%
forecast(new_data = anti_join(myseries_adj, myseries_train))
accuracy_summary <- bind_rows(
fit %>%
accuracy() %>%
select(.model, .type, RMSE) %>%
mutate(Data = "Training"),
fc %>%
accuracy(myseries_adj) %>%
select(.model, .type, RMSE) %>%
mutate(Data = "Test")
)
accuracy_summary## # A tibble: 2 × 4
## .model .type RMSE Data
## <chr> <chr> <dbl> <chr>
## 1 "ETS(Turnover ~ error(\"M\") + trend(\"A\") + season(\"N\")… Trai… 15.2 Trai…
## 2 "ETS(Turnover ~ error(\"M\") + trend(\"A\") + season(\"N\")… Test 31.9 Test