Consider 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 α and ℓ0, and generate forecasts for the next four months.
# Filter the dataset for Pigs in Victoriavictoria_pigs <- aus_livestock |>filter(Animal =="Pigs", State =="Victoria")# Plot the datavictoria_pigs |>autoplot(Count) +labs(y ="Count", title ="Pigs in Victoria")
# Fit the dataset using Simple Exponential Smoothing with ANNvictoria_pigs_fit <- victoria_pigs |>model(ANN =ETS(Count ~error("A") +trend("N") +season("N")))# Find the optimal values of alpha and l0report(victoria_pigs_fit)
# Generate forecasts for the next four monthsvictoria_pigs_fc <- victoria_pigs_fit |>forecast(h =4)print(victoria_pigs_fc)
# A fable: 4 x 6 [1M]
# Key: Animal, State, .model [1]
Animal State .model Month
<fct> <fct> <chr> <mth>
1 Pigs Victoria ANN 2019 Jan
2 Pigs Victoria ANN 2019 Feb
3 Pigs Victoria ANN 2019 Mar
4 Pigs Victoria ANN 2019 Apr
# ℹ 2 more variables: Count <dist>, .mean <dbl>
# Plot the forecast for the next four monthsvictoria_pigs_fc |>autoplot(victoria_pigs) +geom_line(aes(y = .fitted), col="#D55E00",data =augment(victoria_pigs_fit)) +labs(y ="Count", title ="Pigs Slaughtered in Victoria")
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.
# Calculate standard deviation of the residuals (s)victoria_pigs_residuals_data <-augment(victoria_pigs_fit)victoria_pigs_s <-sd(victoria_pigs_residuals_data$.resid, na.rm =TRUE)# Get the point forecast for the first month (y_hat)victoria_pigs_y_hat <- victoria_pigs_fc$.mean[1]# Compute the 95% prediction intervalvictoria_pigs_manual_lower <- victoria_pigs_y_hat - (1.96* victoria_pigs_s)victoria_pigs_manual_upper <- victoria_pigs_y_hat + (1.96* victoria_pigs_s)cat("95% PI: [", victoria_pigs_manual_lower, ", ", victoria_pigs_manual_upper, "]\n")
95% PI: [ 76871.01 , 113502.1 ]
# Extract the interval produced by R for comparisonvictoria_pigs_r_interval <- victoria_pigs_fc |>slice(1) |>hilo(level =95)print(victoria_pigs_r_interval)
# A tsibble: 1 x 7 [1M]
# Key: Animal, State, .model [1]
Animal State .model Month
<fct> <fct> <chr> <mth>
1 Pigs Victoria ANN 2019 Jan
# ℹ 3 more variables: Count <dist>, .mean <dbl>, `95%` <hilo>
The 95% prediction interval from the manual calculation is very similar to the 95% prediction interval generated by R. There is a very small difference between the values where the R-generated interval is slightly increased by about 15 or less than a hundredth of a percent for both the start and end of the interval.
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.
# Filter the dataset for Exports from Chinachina_exports <- global_economy |>filter(Country =="China") |>select(Exports)# Plot the datachina_exports |>autoplot(Exports) +labs(y ="Exports of goods and services (% of GDP)", title ="Exports from China")
Looking at the plot of exports from China, we can see that it had a strong positive trend for most of the plot and then reversed into a negative trend around 2005/2006 and onwards.
b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# Fit the dataset using Simple Exponential Smoothingchina_exports_fit <- china_exports |>model(ANN =ETS(Exports ~error("A") +trend("N") +season("N")))# Find the optimal values of alpha and l0report(china_exports_fit)
# Generate forecasts for the next twelve monthschina_exports_fc <- china_exports_fit |>forecast(h =12)print(china_exports_fc)
# A fable: 12 x 4 [1Y]
# Key: .model [1]
.model Year
<chr> <dbl>
1 ANN 2018
2 ANN 2019
3 ANN 2020
4 ANN 2021
5 ANN 2022
6 ANN 2023
7 ANN 2024
8 ANN 2025
9 ANN 2026
10 ANN 2027
11 ANN 2028
12 ANN 2029
# ℹ 2 more variables: Exports <dist>, .mean <dbl>
# Plot the forecast for the next twelve monthschina_exports_fc |>autoplot(china_exports) +geom_line(aes(y = .fitted), col="#D55E00",data =augment(china_exports_fit)) +labs(y ="Exports of goods and services (% of GDP)", title ="Exports from China")
c. Compute the RMSE values for the training data.
# Calculate the RMSE value for the training dataaccuracy(china_exports_fit) |>select(.model, RMSE)
# A tibble: 1 × 2
.model RMSE
<chr> <dbl>
1 ANN 1.90
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 the dataset using AANchina_exports_fit2 <- china_exports |>model(AAN =ETS(Exports ~error("A") +trend("A") +season("N")))# Find the optimal values of alpha and l0 for AANchina_exports_fit2 |>report()
# Generate forecasts for the next twelve monthschina_exports_fc2 <- china_exports_fit2 |>forecast(h =12)print(china_exports_fc2)
# A fable: 12 x 4 [1Y]
# Key: .model [1]
.model Year
<chr> <dbl>
1 AAN 2018
2 AAN 2019
3 AAN 2020
4 AAN 2021
5 AAN 2022
6 AAN 2023
7 AAN 2024
8 AAN 2025
9 AAN 2026
10 AAN 2027
11 AAN 2028
12 AAN 2029
# ℹ 2 more variables: Exports <dist>, .mean <dbl>
# Plot the forecast for the next twelve monthschina_exports_fc2 |>autoplot(china_exports) +geom_line(aes(y = .fitted), col="#D55E00",data =augment(china_exports_fit2)) +labs(y ="Exports of goods and services (% of GDP)", title ="Exports from China")
The plot for the ANN model shows a more symmetrical forecast with a straight line and even area above and below the line. The plot for AAN is downward trending because the dataset was trending downwards towards the end so the forecast also is downward trending with more emphasis on lower values.
The advantage of the AAN model is being able to take into account the most recent trend because it takes into account an extra variable. While this could be beneficial, it could also be detrimental if that recent downward trend was a one-off event like the Covid lock-down. The ANN model gives a more neutral approach where it doesn’t factor in the latest trend so I believe that’s a safer forecast.
e. Compare the forecasts from both methods. Which do you think is best?
If I had to choose between the two forecast approaches, I would choose the ANN model just because it gives a balanced view to the forecast rather than an biased approached. If another black swan event like Covid happened again, then it would make the forecasts from the AAN model irrelevant. But even just a rare recent dip could also affect the accuracy of the AAN model as well.
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.
# Calculate RMSE values from accuracy metrics for the ANN modelchina_exports_metrics <-accuracy(china_exports_fit)rmse_ann <- china_exports_metrics |>pull(RMSE)# Slice the first forecast pointchina_exports_fc_first <- china_exports_fc |>slice(1)yhat_ann <- china_exports_fc_first |>pull(.mean)# Calculate 95% intervalscat("Manual ANN 95% PI: [", yhat_ann -1.96* rmse_ann, ", ", yhat_ann +1.96* rmse_ann, "]\n")
# A tsibble: 1 x 5 [1Y]
# Key: .model [1]
.model Year Exports
<chr> <dbl> <dist>
1 AAN 2018 N(19, 3.9)
# ℹ 2 more variables: .mean <dbl>, `95%` <hilo>
For the AAN model, the manually calculated prediction interval was also a bit smaller than the interval produced by R.
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.]
# Filter the dataset for China's GDPchina_gdp <- global_economy |>filter(Country =="China") |>select(GDP)# Plot the datachina_gdp |>autoplot(GDP) +labs(y ="GDP ($USD)", title ="China's GDP")
# Calculate the optimal lambda parameter for Box-Cox using the Guerrero methodlambda <- china_gdp |>features(GDP, features = guerrero) |>pull(lambda_guerrero)# compare standard ETS, Damped Trend Method, and Box-Cox Transformationchina_gdp_fit <- china_gdp |>model(`Standard ETS`=ETS(GDP),`Damped Trend`=ETS(GDP ~trend("Ad")),`Box-Cox Transformation`=ETS(box_cox(GDP, lambda)) )# Generate forecasts for the next twenty yearschina_gdp_fc <- china_gdp_fit |>forecast(h =20)# Plot the forecastschina_gdp_fc |>autoplot(china_gdp, level =NULL) +labs(y ="GDP ($USD)", title ="China's GDP Forecasts")
The comparison plot shows a comparison of various different ETS methods. In the plot, the Box-Cox Transformation using an optimal lambda parameter from the Guerrero method is the one with the biggest difference. Since China’s GDP was exponentially increasing, the Box-Cox transformation takes it to another level and continues that trend. Meanwhile the Standard ETS takes the trend and continues indefinitely with the positive trend. The Damped Trend method takes a more conservative approach and initially trends upwards but slows down or dampens as it goes. I think the Damped Trend method makes the most sense in this case because it would be wishful thinking to believe that the GDP will keep increasing at the current rate into the future or sustain the exponential growth shown with the Box-Cox Transformation.
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?
# Plot the raw data to observe the seasonalityaus_production |>autoplot(Gas) +labs(y ="Petajoules", title ="Australian Gas Production")
# Fit an automatic ETS model and a damped trend modelaus_production_fit <- aus_production |>model(`Standard ETS`=ETS(Gas),`Damped Trend`=ETS(Gas ~trend("Ad") +season("M")) )
# Forecast the next five yearsaus_production_fc <- aus_production_fit |>forecast(h ="5 years")
# Plot the forecastsaus_production_fc |>autoplot(aus_production) +labs(title ="Forecasts of Australian Gas Production", y ="Petajoules")
# Check model accuracy and information criteriaaccuracy(aus_production_fit)
# A tibble: 2 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Standard ETS Training -0.115 4.60 3.02 0.199 4.08 0.542 0.606 -0.0131
2 Damped Trend Training -0.00439 4.59 3.03 0.326 4.10 0.544 0.606 -0.0217
The multiplicative seasonality is necessary for this dataset because the seasonal peaks and troughs get bigger as the data trends positive. This means that an additive seasonality would not be appropriate since usually that’s for a dataset where the peaks and troughs are more uniform.
In order to test if the Damped Trend model improved forecasts, I compared the RMSE and AICc values for both the Standard ETS and the Damped Trend model. Based on the AICc values, the Standard ETS model is the better choice since it has a lower value, 1681.794 vs 1685.091 for the Damped Trend model. Based on the RMSE values, they were very similar with less than a 0.004 difference. I think in this case, the simpler model, the Standard ETS model is the better choice.
Exercise 8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
set.seed(8585)myseries <- aus_retail |>filter(`Series ID`==sample(aus_retail$`Series ID`,1))myseries |>autoplot(Turnover) +labs(y ="Turnover", title ="Other retailing n.e.c. in New South Wales")
a. Why is multiplicative seasonality necessary for this series?
Multiplicative seasonality is necessary for this series because the height of the peaks and troughs get larger as the dataset trends positively. Additive seasonality would not be appropriate for this dataset because that is used for seasonality that is more even.
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
# Apply Holt-Winters' multiplicative methodmyseries_fit <- myseries |>model(`HW Multiplicative`=ETS(Turnover ~error("M") +trend("A") +season("M")) )# Forecast the next five yearsmyseries_fc <- myseries_fit |>forecast(h ="5 years")# Plot the forecastsmyseries_fc |>autoplot(aus_retail) +labs(y ="Turnover", title ="Other retailing n.e.c. in New South Wales")
# Apply Holt-Winters' Damped multiplicative methodmyseries_fit_damped <- myseries |>model(`HW Damped`=ETS(Turnover ~error("M") +trend("Ad") +season("M")) )# Forecast the next five yearsmyseries_fc_damped <- myseries_fit_damped |>forecast(h ="5 years")# Plot the forecastsmyseries_fc_damped |>autoplot(aus_retail) +labs(y ="Turnover", title ="HW Damped Other retailing n.e.c. in New South Wales")
These plots show both the Holt-Winters’ Multiplicative method as well as the Holt-Winters’ Damped method. The HW Damped has a slightly more conservative trend since the goal of the Damped method is to eventually become a flat line.
c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
# Calculate RMSE for the HW Multiplicativeaccuracy(myseries_fit) |>select(.model, RMSE)
I prefer the HW Multiplicative method because it has a lower RMSE score which means that it fits the historical data better. With that said, I can see the benefit of the HW Damped method because it eventually flattens out instead of increasing indefinitely which can lead to unrealistic forecasts long term.
d. Check that the residuals from the best method look like white noise.
The residuals look like white noise for the standard HW Multiplicative method. The acf plot shows only 1 spike at the very beginning that goes past the dashed blue line. The histogram sounds a relatively normal distribution around 0. The residuals plot shows random scatter around 0 with no clear pattern.
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?
# Create the training model to the end of 2010myseries_training <- myseries |>filter_index(. ~"2010 Dec")# Train all three models on the training datatraining_fit <- myseries_training |>model(`HW Multiplicative`=ETS(Turnover ~error("M") +trend("A") +season("M")),`HW Damped`=ETS(Turnover ~error("M") +trend("Ad") +season("M")),`Seasonal Naive`=SNAIVE(Turnover) )# Forecast over the test periodtest_data <-anti_join(myseries, myseries_training, by =c("State", "Industry", "Month"))training_fc <- training_fit |>forecast(new_data = test_data)# Compare the test set RMSEaccuracy(training_fc, myseries) |>select(.model, RMSE) |>arrange(RMSE)
The standard HW Multiplicative model has the lowest RMSE score. It beats out both the SNAIVE model and the HW Damped model.
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?
# Find the optimal Box-Cox lambda parameter on the training datalambda <- myseries_training |>features(Turnover, features = guerrero) |>pull(lambda_guerrero)# Fit the STL decomposition model alongside the previous best HW modelmyseries_fit_comparison <- myseries_training |>model(`HW Multiplicative`=ETS(Turnover ~error("M") +trend("A") +season("M")),`STL + ETS`=decomposition_model(STL(box_cox(Turnover, lambda)),ETS(season_adjust ~season("N")) ) )# Forecast using the test data datesmyseries_fc_comparison <- myseries_fit_comparison |>forecast(new_data = test_data)# Compare the test set RMSEaccuracy(myseries_fc_comparison, myseries) |>select(.model, RMSE) |>arrange(RMSE)
The RMSE score of the STL + ETS model was lower at 187.5048 than the previous best standard HW Multiplicative model’s RMSE score of 193.9849. This means that the STL + ETS model fits the training data better.
