Imports

library(gridExtra)
library(ggplot2)
library(cowplot)
options(scipen=10000)
library(mlbench)
library(tidyverse)
library(corrplot)
library(AppliedPredictiveModeling)
library(caret)
library(DataExplorer)
library(kableExtra)
library(mice)
library(fpp3)

8.1

Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.

Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of
αlpha and ℓ0, and generate forecasts for the next four months.

pigs <- aus_livestock |> filter(Animal == "Pigs" & State == "Victoria")  |> dplyr::select(Count)
fit_model <- pigs |> model(ETS(Count~error("A")+trend("N")+season("N")))

pigs_fr <- fit_model |> forecast(h="4 months")
pigs_fr |> autoplot(pigs) + geom_line(aes(y=.fitted),data = augment(fit_model) ) + ggtitle("Number of Pigs slaughtered in Victoria")

Compute a 95% prediction interval for the first forecast using yt ± 1.96 “s” where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.

forecast_internvals = pigs_fr |> hilo(level = 95) |> mutate(lower = `95%`$lower,upper=`95%`$upper) |> mutate(ProducedBy="R")

prediction_intervals <- forecast_internvals[1,c('ProducedBy','lower','upper')] 
prediction_intervals[2, "lower"] <- forecast_internvals[1,]$.mean - 1.96*sd(augment(fit_model)$.resid)
prediction_intervals[2, "upper"] <-  forecast_internvals[1,]$.mean + 1.96*sd(augment(fit_model)$.resid)
prediction_intervals[2, "ProducedBy"] <- "Manual"

prediction_intervals |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
ProducedBy lower upper
R 76854.79 113518.3
Manual 76871.01 113502.1

The forecast intervals produced by R and your manual calculations are nearly identical, with only slight differences. The lower bounds differ by 16.22, and the upper bounds by 16.2. These small variations are likely due to rounding or differences in computational precision between R and manual methods. Despite these minimal discrepancies, the intervals are practically the same and would not significantly affect the interpretation of the forecast results.

8.5

Data set global_economy contains the annual Exports from many countries. Select one country to analyse.

Plot the Exports series and discuss the main features of the data.

usa_economy <- global_economy |> filter(Country == "United States" & Year <= 2016) |> dplyr::select(Exports)
usa_economy |> autoplot(Exports)

usa_economy_train <- usa_economy |> slice(1: (nrow(usa_economy)*0.7)  )
usa_economy_test <- usa_economy |> filter(Year > max(usa_economy_train$Year))

Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.

test_length <- nrow(usa_economy_test)
usa_model_fit <- usa_economy_train |> model(ANN = ETS(Exports~error("A")+trend("N")+season("N")))

usa_forecast <- usa_model_fit |> forecast(h=test_length)

usa_forecast |> autoplot() + autolayer(usa_economy_train,Exports)

Compute the RMSE values for the training data.

accuracy(usa_model_fit) |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
ANN Training 0.1414279 0.5318413 0.3771615 1.685702 4.710297 0.9744113 0.9871433 0.4272424

Compare the results to those from an ETS(A,A,N) model. Discuss the merits of the two forecasting methods for this data set.

usa_model_fit_trended <- usa_economy_train |> model(AAN = ETS(Exports~error("A")+trend("A")+season("N")))

usa_forecast_trended <- usa_model_fit_trended |> forecast(h=test_length)



usa_model_results <- accuracy(usa_model_fit)

usa_model_results <- rbind(usa_model_results, accuracy(usa_model_fit_trended))

usa_model_results |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
ANN Training 0.1414279 0.5318413 0.3771615 1.6857023 4.710297 0.9744113 0.9871433 0.4272424
AAN Training -0.0085090 0.5164241 0.3955093 -0.3757971 5.220122 1.0218136 0.9585278 0.4086637 
# usa_forecast_trended |> autoplot() + autolayer(usa_economy_train,Exports)

The ANN and the AAN models offer distinct advantages depending on the dataset’s characteristics. The AAN model, which includes a linear trend, has a lower RMSE (0.5164) and near-zero mean error (ME: -0.0085), indicating it handles data with trends better by producing less biased forecasts and fewer large errors. Its slightly lower autocorrelation in residuals (ACF1: 0.4087) also points to more independent errors. In contrast, the ANN model, which assumes no trend, has a lower MAPE (4.71%) and MASE (0.9744), making it marginally more accurate in terms of percentage-based and scaled errors. It also has a slightly lower MAE, suggesting it may be more suitable for stable datasets without significant trends. The choice between these models hinges on whether the data contains an underlying trend or remains stable over time.

Compare the forecasts from both methods. Which do you think is best?

usa_forecast_results <- accuracy(usa_forecast, usa_economy)
usa_forecast_results <- rbind(usa_forecast_results, accuracy(usa_forecast_trended, usa_economy))

usa_forecast_results  |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
ANN Test 0.9194272 1.822508 1.495330 6.275145 12.410716 3.863244 3.382733 0.8561814
AAN Test -0.6407437 1.197360 1.014339 -6.814088 9.587499 2.620584 2.222404 0.6711062 
plot1 <- usa_forecast |> autoplot() + autolayer(usa_economy,Exports) + ggtitle("USA Exports (A,N,N)")
plot2 <- usa_forecast_trended |> autoplot() + autolayer(usa_economy,Exports) + ggtitle("USA Exports (A,A,N)")

plot_grid(plot1,plot2)

The AAN model outperforms the ANN model across most key metrics. The AAN model has a lower RMSE (1.1974 vs. 1.8225), indicating fewer large errors, and a lower MAE (1.0143 vs. 1.4953), showing more accurate average forecasts. Additionally, its MAPE (9.59%) is better than the ANN model’s (12.41%), meaning the AAN model is more accurate percentage-wise. While the AAN model slightly underestimates (ME: -0.6407), it performs better overall and has lower autocorrelation in its errors. Thus, the AAN model is the better forecasting method for this dataset.

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.

usa_ann_internval = usa_forecast |> hilo(level = 95) |> mutate(lower = `95%`$lower,upper=`95%`$upper)
usa_aan_internval = usa_forecast_trended |> hilo(level = 95) |> mutate(lower = `95%`$lower,upper=`95%`$upper)

usa_interval_ann <- usa_ann_internval |> mutate(Model = "ANN",Producedby="R")
usa_interval_aan <- usa_aan_internval |> mutate(Model = "AAN",Producedby="R")

interval_com <- usa_interval_ann[1,c('Model','Producedby','lower','upper')] 
interval_com[2, ] <- usa_interval_aan[1,c('Model','Producedby','lower','upper')] 



aug_model <- augment(usa_model_fit)

interval_com[3, "lower"] <- usa_forecast[1,]$.mean - 1.96*sd(augment(usa_model_fit)$.resid) 
interval_com[3, "upper"] <-  usa_forecast[1,]$.mean + 1.96*sd(augment(usa_model_fit)$.resid) 
interval_com[3, "Producedby"] <- "Manual"
interval_com[3, "Model"] <- "ANN"

interval_com[4, "lower"] <- usa_forecast_trended[1,]$.mean - 1.96*sd(augment(usa_model_fit_trended)$.resid) 
interval_com[4, "upper"] <-  usa_forecast_trended[1,]$.mean + 1.96*sd(augment(usa_model_fit_trended)$.resid) 
interval_com[4, "Producedby"] <- "Manual"
interval_com[4, "Model"] <- "AAN"

interval_com |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
Model Producedby lower upper
ANN R 9.414667 11.55505
AAN R 9.580655 11.71755
ANN Manual 9.466846 11.50287
AAN Manual 9.623818 11.67439

The 95% prediction intervals calculated manually for both the ANN and AAN models are nearly identical to those produced by R, with only minor differences in the upper and lower bounds. For the ANN model, the manual interval is (9.466846, 11.50287), while R’s interval is (9.414667, 11.55505). Similarly, for the AAN model, the manual interval is (9.623818, 11.67439), and R’s interval is (9.580655, 11.71755). These small discrepancies are likely due to rounding or precision differences in the calculations, but overall, the results are consistent, demonstrating that the manual method aligns closely with R’s output.

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.

china_gdp <- global_economy |> filter(Country == "China") |> dplyr::select(GDP)
china_gdp |> autoplot(GDP) + ggtitle("China GDP")

china_model_fit <- china_gdp |> model(
  ANN = ETS(GDP~error("A")+trend("N")+season("N")),
  AAN = ETS(GDP~error("A")+trend("A")+season("N")),
  ADN = ETS(GDP~error("A")+trend("Ad")+season("N"))
  )

china_forecasts <- china_model_fit |> forecast(h=50)


china_forecasts |> autoplot() + autolayer(china_gdp,GDP) +facet_grid(~.model)

The AAN model you describe forecasts a steady linear trend, but its prediction intervals provide important insights into the model’s level of uncertainty. While the upper bound of the interval increases linearly, the lower bound follows a logarithmic growth pattern, eventually curving downward and dipping below the last observed value. This suggests that, although the model predicts continuous growth, it remains cautious about the possibility of declines. The downward curvature of the lower bound implies that the model sees potential for sharp downturns or underperformance, likely due to volatility or past sudden drops in the data. Ultimately, while the forecast projects growth, the model incorporates significant downside risk, reflecting a conservative approach to future uncertainty.

In contrast, the damped AAdN model produces forecasts with a curving line, indicating that the trend is expected to slow over time, following a pattern of logarithmic growth. This behavior reflects the model’s assumption that, while the trend will continue, it will taper off gradually rather than maintain a constant rate of growth. The prediction intervals mirror this pattern, showing that the model accounts for uncertainty, but it anticipates that uncertainty will grow more gradually as time progresses. Compared to the non-damped AAN model, this approach is more conservative, suggesting a more stable long-term outlook with less volatility. The damped model is particularly useful when there is an expectation that the trend will slow down, offering a more balanced and cautious forecast that accounts for future stabilization.

The ANN model, on the other hand, produces a forecast with a horizontal line, meaning it detects no trend and expects future values to remain constant over time. This indicates that the model assumes stability in the data without anticipating any significant upward or downward movements. The prediction intervals in this model grow linearly, reflecting increasing uncertainty as the forecast horizon extends, even though the forecasted values remain flat. This behavior suggests that while the model predicts no change in the future, it still acknowledges growing uncertainty the further out it forecasts. Overall, the ETS(ANN) model is best suited for datasets where no clear trend exists, and future values are expected to remain stable over time.

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_gas <- aus_production |> dplyr::select(Gas)
aus_gas_train <- aus_gas |> slice(1: (nrow(aus_gas)*0.7)  )
aus_gas_train |> autoplot(Gas) + ggtitle("Australia Gas Training Data")

aus_gas_fit <- aus_gas_train |> model(
  # AAA = ETS(Gas~error("A")+trend("A")+season("A")),
  AAM = ETS(Gas~error("A")+trend("A")+season("M")),

  ADM = ETS(Gas~error("A")+trend("Ad")+season("M"))

  )
forecastRange <- nrow(aus_gas) - nrow(aus_gas_train) 
aus_gas_forecasts <- aus_gas_fit |> forecast(h= forecastRange)
aus_gas_forecasts |> autoplot() + autolayer(aus_gas,Gas)+facet_grid(.model~.)+ ggtitle("Australia Gas")

accuracy(aus_gas_forecasts, aus_gas) |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
AAM Test 8.813859 11.98350 9.519342 4.444384 4.836461 1.940582 1.693807 0.6598361
ADM Test 25.757802 31.18278 25.789631 12.780896 12.799275 5.257390 4.407529 0.9092521

In this experiment, the ETS model uses an additive trend and multiplicative seasonality. Multiplicative seasonality is necessary because the amplitude of seasonal fluctuations increases or decreases in proportion to the level of the trend. As the trend grows or declines, the seasonal variations become more pronounced or diminish accordingly. In addition to the multiplicative seasonality, the ETS model includes a dampened trend, which curves downward over time. This means that while the gas production initially shows an increasing trend, the model expects the growth rate to slow down gradually.

On the other hand, the comparison between the AAM and ADM models shows that the AAM model performs better overall, with lower error metrics across the board. It has a lower RMSE (11.18) and MAE (8.83) compared to the ADM model’s higher RMSE (31.00) and MAE (25.68), indicating smaller forecast errors. The AAM model also has lower MAPE (4.50%) compared to ADM’s 12.75%, suggesting more accurate percentage-based predictions. Additionally, AAM shows less autocorrelation in residuals (ACF1: 0.637 vs. ADM’s 0.908), indicating more reliable error independence. Overall, AAM offers better predictive performance.

Therefore, multiplicative seasonality is essential to model the dynamic relationship between the trend and seasonal variations effectively.

8.8

Recall your retail time series data (from Exercise 7 in Section 2.10).

set.seed(12345678)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1)) |> dplyr::select(Turnover)
myseries |> autoplot(Turnover)+labs(title = "Autralian Retail Turnover", x="Month")

Why is multiplicative seasonality necessary for this series?

Multiplicative seasonality is necessary because the data exhibits strong seasonal patterns, with the amplitude of seasonal fluctuations gradually increasing over time. This variation suggests that the seasonality grows in proportion to the trend, making multiplicative seasonality more appropriate for capturing the changing seasonal effects.

Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.

myseries_model_fit <- myseries |> model(
  MAM = ETS(Turnover~error("M")+trend("A")+season("M")),
  ADN = ETS(Turnover~error("M")+trend("Ad")+season("M"))
  )

myseries_forecasts <- myseries_model_fit |> forecast(h=30)


myseries_forecasts |> autoplot() + autolayer(myseries,Turnover) +facet_grid(~.model)

Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

accuracy(myseries_model_fit) |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
MAM Training -0.0128385 0.6130408 0.4497149 -0.4689010 5.147100 0.5134257 0.5290064 -0.0094057
ADN Training 0.0397571 0.6155067 0.4444157 -0.0722921 5.061188 0.5073758 0.5311343 -0.0436615

While both models perform similarly in terms of RMSE, the ADN model has a slight edge in MAE, MAPE, and other metrics, making it a slightly better choice overall for one-step forecasts.

Check that the residuals from the best method look like white noise.

myseries_model_fit |> dplyr::select(ADN) |> gg_tsresiduals()

The residuals are normally distributed and resemble white noise, indicating no obvious patterns in the errors. Additionally, the autocorrelation values remain within acceptable thresholds, suggesting that there is no significant correlation between residuals over time.

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?

turnover_train <-  myseries |> filter(year(Month) <= 2010)
turnover_test <-  myseries |> filter(year(Month) > 2010)
turnover_model_fit <- turnover_train |> model(
  MAM = ETS(Turnover~error("M")+trend("A")+season("M")),
  Snaive= SNAIVE(Turnover)
  )

turnover_forecasts <- turnover_model_fit |> forecast(h= nrow(turnover_test))
turnover_forecasts |> autoplot() + autolayer(myseries,Turnover)  +facet_grid(~.model)

# myseries |> autoplot(Turnover)+labs(title = "Autralian Retail Turnover", x="Month")
turnover_forecast_snaive_man <- accuracy(turnover_forecasts, myseries)
turnover_forecast_snaive_man |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
MAM Test -1.5388814 1.782947 1.597517 -12.292535 12.642079 1.746761 1.468981 0.4951790
Snaive Test 0.8364583 1.551981 1.240625 5.940129 9.064092 1.356528 1.278687 0.6011751

The comparison between the MAM and seasonal naïve (Snaive) models shows that the Snaive model performs better, with a lower RMSE (1.5520 vs. 1.7829) and more accurate MAE (1.2406 vs. 1.5975) and MAPE (9.06% vs. 12.64%). This indicates that the Snaive model captures seasonal patterns more effectively, providing more precise forecasts. Despite the MAM model’s complexity, it cannot outperform the simplicity of the Snaive model in this case, making the Snaive model the better choice for forecasting.

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 <- turnover_train |> features(Turnover,features = "guerrero") |> pull(lambda_guerrero)

turnover_stl_fit <- turnover_train |> model(STL(box_cox(Turnover,lambda) ~ trend() + season(window="periodic") ))
turnover_stl_dec <- turnover_stl_fit |> components() |> dplyr::select(season_adjust)

turnover_ets_fit <- turnover_stl_dec |> model(ETS(season_adjust))


turnover_ets_boxcox_forecast <-  turnover_ets_fit |> forecast(h= 10)

turnover_decomposed_test <-  myseries |> model(STL(box_cox(Turnover,lambda) ~ trend() + season(window="periodic") )) |> components() |> dplyr::select(season_adjust)

Comparison Results

model_result <- accuracy(turnover_ets_boxcox_forecast,turnover_decomposed_test) |>
  rbind(turnover_forecast_snaive_man[2,])

model_result  |> kable() %>% kable_styling() %>%  kable_classic(full_width = F)
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
ETS(season_adjust) Test -0.1738302 0.193366 0.1738302 -5.509745 5.509745 0.9187345 0.7700024 -0.0544203
Snaive Test 0.8364583 1.551981 1.2406250 5.940129 9.064092 1.3565275 1.2786869 0.6011751

When comparing the ETS(season_adjust) model to the Snaive model on relevant scale-independent metrics, ETS(season_adjust) consistently outperforms despite the seasonal component being transformed using the Box-Cox method. The ETS model shows a significantly lower MASE (0.92 vs. 1.36) and RMSSE (0.77 vs. 1.28), indicating smaller errors relative to the scale and baseline. Additionally, the ACF1 for ETS (-0.05) suggests almost no autocorrelation in the residuals, pointing to well-behaved, independent errors, whereas the Snaive model has a high ACF1 (0.60), indicating residual patterns and less reliable forecasts. Overall, ETS(season_adjust) offers superior accuracy and reliability compared to the Snaive model, making it the better choice even with the transformation applied.