Consider 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 α and ℓ0, and generate forecasts for the next four months.
# Let's visualize the aus_livestock dataset head(aus_livestock)
# A tsibble: 6 x 4 [1M]
# Key: Animal, State [1]
Month Animal State Count
<mth> <fct> <fct> <dbl>
1 1976 Jul Bulls, bullocks and steers Australian Capital Territory 2300
2 1976 Aug Bulls, bullocks and steers Australian Capital Territory 2100
3 1976 Sep Bulls, bullocks and steers Australian Capital Territory 2100
4 1976 Oct Bulls, bullocks and steers Australian Capital Territory 1900
5 1976 Nov Bulls, bullocks and steers Australian Capital Territory 2100
6 1976 Dec Bulls, bullocks and steers Australian Capital Territory 1800
# Let's get the number of pigs slaughtered in VictoriaVictoria_pigs <- aus_livestock |>filter(Animal=="Pigs", State=="Victoria")head(Victoria_pigs)
# A tsibble: 6 x 4 [1M]
# Key: Animal, State [1]
Month Animal State Count
<mth> <fct> <fct> <dbl>
1 1972 Jul Pigs Victoria 94200
2 1972 Aug Pigs Victoria 105700
3 1972 Sep Pigs Victoria 96500
4 1972 Oct Pigs Victoria 117100
5 1972 Nov Pigs Victoria 104600
6 1972 Dec Pigs Victoria 100500
# Let's plot the series to see pattern of the number of pigs slaughteredVictoria_pigs |>autoplot(Count) +labs(title="Number of Pigs slaughtered in Victoria")
# Let's Use the ETS() function to estimate the equivalent model for simple exponential smoothing and Find the optimal values of α and ℓ0Victoria_pigs_fit <- Victoria_pigs |>model(SES =ETS(Count ~error("A") +trend("N") +season("N")))report(Victoria_pigs_fit)
## Let's generate the forecast for the next four monthspigs_forecats <- Victoria_pigs_fit |>forecast(h=4)pigs_forecats
# A fable: 4 x 6 [1M]
# Key: Animal, State, .model [1]
Animal State .model Month
<fct> <fct> <chr> <mth>
1 Pigs Victoria SES 2019 Jan
2 Pigs Victoria SES 2019 Feb
3 Pigs Victoria SES 2019 Mar
4 Pigs Victoria SES 2019 Apr
# ℹ 2 more variables: Count <dist>, .mean <dbl>
Results
The optimal values of alpha (α) is 0.3221247 and the initial state ℓ0 is 100646.6.
Compute a 95% prediction interval for the first forecast using y_hat±1.96s where
s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
** First step:** Let’s manually compute a 95% prediction interval for the first forecast
# Let's calculate s is the standard deviation of the residuals.res <-residuals(Victoria_pigs_fit, type ="response")$.resids <-sd(res, na.rm =TRUE) s
[1] 9344.666
# Let's get the first forecast valuefirst_forecast <- pigs_forecats |>as_tibble() |>slice(1) |>pull(.mean)first_forecast
[1] 95186.56
# Let's manually calculate the lower and upper value of the 95% prediction interval for the first forecastInt_lower_value <- first_forecast -1.96*sInt_upper_value <- first_forecast +1.96*sInt_lower_value
[1] 76871.01
Int_upper_value
[1] 113502.1
# Let's compute the R's first forecast intervalfirst_forecast <- pigs_forecats |>hilo(level =95) |>head(1)R_first_forecast_Interval <- first_forecast[7]R_first_forecast_Interval
The manually computed interval for the first forecast is equivalent to [76871.01 ,113502.1] whereas the R’s 95% prediction interval is [76854.79, 113518.3].The lower value differ by only 16.22 and the upper value by 16.20. It was predictable that the interval would not be exactly the same because he manual formula is a simple approximation, while R’s interval accounts for the ETS model structure and it uses the exact critical value for a standard normal distribution to calculate the 1-step forecast variance rather than the sample standard deviation s of the in-sample residual series. In conclusion,both interval are close to each other.
EXERCISE 8.5
Data set global_economy contains the annual Exports from many countries. Let’s select one country to analyse the following. I will choose cameroon for my analysis.
a) Plot the Exports series and discuss the main features of the data.
# Dataset visualizationhead(global_economy)
# A tsibble: 6 x 9 [1Y]
# Key: Country [1]
Country Code Year GDP Growth CPI Imports Exports Population
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Afghanistan AFG 1960 537777811. NA NA 7.02 4.13 8996351
2 Afghanistan AFG 1961 548888896. NA NA 8.10 4.45 9166764
3 Afghanistan AFG 1962 546666678. NA NA 9.35 4.88 9345868
4 Afghanistan AFG 1963 751111191. NA NA 16.9 9.17 9533954
5 Afghanistan AFG 1964 800000044. NA NA 18.1 8.89 9731361
6 Afghanistan AFG 1965 1006666638. NA NA 21.4 11.3 9938414
# Let's filter data for CameroonCameroon_economy <- global_economy |>filter(Country=="Cameroon") |>select(Year, Exports) |>drop_na()Cameroon_economy
# Let's plot the Exports series and discuss the main features of the data.Cameroon_economy |>autoplot(Exports) +labs(title="Cameroon Exports as % of GDP from 1960 - 2017")
Interpretation
Cameroon exports (% of GDP) fluctuate between 15% and 33% with no persistent trend. The series does not show a seasonal cycle which can be expected for annual data. Moreover, Cameroon exports jump to around 33% in 1982 to 1985, then we observe a huge drop around 16 to 17% in 1987-1988 with a dip less than 15% in 1993 stands out as well. Since the mid-1990s the series shows a type of cyclical fluctuations in a narrower band around 20% to 29% and has been drifting down since 2012 and ending around 18% in 2017.
b) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
# Fit ETS modelfit <- Cameroon_economy |>model(ETS_ANN =ETS(Exports ~error("A") +trend("N") +season("N")))# Model summaryreport(fit)
# Let's Forecast for the next 5 years and plot the forecastsforecast_annual <- fit |>forecast(h=5)forecast_annual |>autoplot(Cameroon_economy) +labs(title="Cameroon exports: ETS(A,N,N) forecast")
Results
The values of alpha (α) is 0.713174 and the initial state ℓ0 is 22.97788 for the ETS ANN model. We know that a high value alpha (α) means the forecast depends in majority on the latest observations. The forecast exhibit a flat line around 19% for the next 5 years (2018-2022).
**c) Compute the RMSE values for the training data.
# Let's find the RMES in the accuracy measuresRMSE_ETS_ANN <-accuracy(fit) |>select(RMSE)RMSE_ETS_ANN
# A tibble: 1 × 1
RMSE
<dbl>
1 3.67
The RMSE values for the training data of the ETS(A,N,N) model is 3.667707.
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.)
# Let's fit the ETS models AANfit_AAN <- Cameroon_economy |>model(ETS_ANN =ETS(Exports ~error("A") +trend("A") +season("N")))report(fit_AAN)
# Let's fit both ETS models for a better comparisonfit_both_ETS <- Cameroon_economy |>model(ETS_ANN =ETS(Exports ~error("A") +trend("N") +season("N")),ETS_AAN =ETS(Exports ~error("A") +trend("A") +season("N")) )# Model summaryreport(fit_both_ETS)
Warning in report.mdl_df(fit_both_ETS): 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.
# Let's compare the training accuracy measureaccurary_measure <-accuracy(fit_both_ETS) accurary_measure
# 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 ETS_ANN Training -0.103 3.67 2.77 -2.18 12.5 0.998 0.963 0.0410
2 ETS_AAN Training 0.00421 3.67 2.79 -1.72 12.5 1.01 0.964 0.0468
Results
The ETS(A,A,N) model uses two extra parameters (β and the initial slope), and that hurts it on the information criteria (AIC, AICc and BIC). Indeed, the ETS(A,N,N) model informations criteria values are lower than their values in the ETS(A,A,N) model because we know lower is better on all three, so the simple model wins.
The values of alpha (α) and beta are respectively 0.7034786 and 0.0001000191 with both initial states ℓ0 and b(0) or slope are respectively 22.20102 and -0.062012 for the ETS AAN model. In the ETS(A,A,N) model has beta close to zero, the trend is nearly constant and their RMSE improves to 3.668642.
e) Compare the forecasts from both methods. Which do you think is best?
# Let's Forecast both models for the next 5 years and plot the forecastsforecast_both_models <- fit_both_ETS |>forecast(h=5)forecast_both_models |>autoplot (Cameroon_economy) +labs(title="Cameroon Exports: ETS(A,N,N) vs ETS(A,A,N)")
INterpretation
ANN exhibit a flat line at around 19%, whereas AAN produces a declining linear slope which follows the most recent observations. The forecast of the next five year for both models are within about 0.5% points of each other which is minimal. In conclusion, the ETS(A,N,N) model is better because it has fewer parameters and better AIC/AICc/BIC.
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.
# Let's pull out the first forecast for each model and training RMSEsRMSE_both <-accuracy(fit_both_ETS) |>select(.model,RMSE)forecast_1 <- forecast_both_models |>filter(Year==2018)# Let's manually calculate 95% prediction intervalsprediction_intervals <- forecast_1 %>%left_join(RMSE_both, by =".model") %>%mutate(Lower95 = .mean -1.96* RMSE,Upper95 = .mean +1.96* RMSE )prediction_intervals
# A fable: 2 x 7 [?]
# Key: .model [2]
.model Year Exports
<chr> <dbl> <dist>
1 ETS_A… 2018 N(19, 15)
2 ETS_A… 2018 N(19, 14)
# ℹ 4 more variables: .mean <dbl>, RMSE <dbl>, Lower95 <dbl>, Upper95 <dbl>
# Let's compute the R's first forecast intervalforecast_1 |>hilo(level=95) |>select(.model, Year, .mean, '95%')
The manual calculation uses the standard normal multiplier 1.96 whereas R computes intervals using the exact quantile (\(z_{0.975} \approx 1.959964\)). For first forecast the variance is just the standard deviation to second power (σ²), so the RMSE-based intervals match the model-based ones to the second decimal. In conclusion, for the first forecast, the manual calculation and those produced using R result will typically match but they would only differ for later horizons, where the model’s intervals widen.
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.]
# Let's filter the data of interest firstChina <- global_economy |>filter(Country=="China")head(China)
# A tsibble: 6 x 9 [1Y]
# Key: Country [1]
Country Code Year GDP Growth CPI Imports Exports Population
<fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 China CHN 1960 59716467625. NA NA 4.43 4.31 667070000
2 China CHN 1961 50056868958. -27.3 NA 3.49 3.87 660330000
3 China CHN 1962 47209359006. -5.58 NA 2.91 4.05 665770000
4 China CHN 1963 50706799903. 10.3 NA 2.86 4.01 682335000
5 China CHN 1964 59708343489. 18.2 NA 2.86 3.77 698355000
6 China CHN 1965 70436266147. 17.0 NA 3.19 3.64 715185000
# Let's determine the optimal Box-Cox transformation parameter using the Guerrero methodlambda <- China |>features(GDP, features = guerrero) |>pull(lambda_guerrero)lambda
[1] -0.03446284
# Let's fit a variety of ETS model specificationsfit_ETS <- China |>model(`Auto (ETS)`=ETS(GDP),`Holt's Linear`=ETS(GDP ~error("A") +trend("A") +season("N")),`Damped Linear`=ETS(GDP ~error("A") +trend("Ad") +season("N")),`Box-Cox Holt`=ETS(box_cox(GDP, lambda) ~error("A") +trend("A") +season("N")),`log`=ETS(log(GDP) ~error("A") +trend("A") +season("N")),`Damped Box-Cox`=ETS(box_cox(GDP, lambda) ~error("A") +trend("Ad") +season("N")) )# Let's generate forecasts with a large value of h=25 and plot the comparisonforecast_25h <- fit_ETS |>forecast(h=25)forecast_25h |>autoplot(China, level=NULL) +labs(title ="China GDPForecasts: ETS Model Variations for h=25",y ="GDP in USD",x ="Year")
Interpretation
The additive trend is extrapolating the estimated trend without damping. Also,damping trend reduces gradually the trend’s contribution, whereas damped Box-Cox combines variance stabilization with a slowing trend.
Intuition to take away
The trend component sets the shape of the forecast including flat, straight, damped or curved.
Damping controls how quickly the projected trend fades, whereas Box-Cox controls the scale on which the trend and errors are modeled.
EXERCISE 8.7
Find an ETS model for the Gas data from aus_production and forecast the next few years.
## Let's filter the data of interest and plot to visualize the gas proudction trend gas_data <- aus_production |>select(Quarter, Gas) |>filter(!is.na(Gas))gas_data |>autoplot(Gas)+labs(title="Quarterly Australian Gas production")
We observe an upward trend in the gas production and a seasonal fluctuations which become larger as the production increases.
Let’s find an ETS model and forecast the next few years.
c) Let’s forecast the next 4 years meaning (16 quarters)
forescast_gas <- fit_gas |>forecast(h =16)# Plot the forecastsforescast_gas |>autoplot(aus_production, level =NULL) +labs(title ="Australian Gas Production Forecasts",subtitle ="ETS models: ordinary versus damped trend",x ="Year",y ="Gas production" )
Result
Multiplicative seasonality is necessary because it assumes that seasonal effects are proportional to the underlying level of the series. In fact, Additive seasonality assumes a fixed-size swing, so it would be too big early on and too small later while A multiplicative form scales the pattern with the level.
Experiment with making the trend damped. Does it improve the forecasts?
For that, we need to compare the models In-sample accuracy and information criteria.
The AICc values for both models are very close. Often, the non dampedmodel has a slightly lower AICc because the data exhibits strong, persistent growth that does not yet show definitive signs of saturation over the sample period. Moreover, the non-damped often performs slightly better on a short period forecast whereas the damping is usually safer because isn’t supported by the likelihood, but it is a reasonable safeguard against the trend overshooting.
# A tsibble: 6 x 5 [1M]
# Key: State, Industry [1]
State Industry `Series ID` Month Turnover
<chr> <chr> <chr> <mth> <dbl>
1 Northern Territory Clothing, footwear and perso… A3349767W 1988 Apr 2.3
2 Northern Territory Clothing, footwear and perso… A3349767W 1988 May 2.9
3 Northern Territory Clothing, footwear and perso… A3349767W 1988 Jun 2.6
4 Northern Territory Clothing, footwear and perso… A3349767W 1988 Jul 2.8
5 Northern Territory Clothing, footwear and perso… A3349767W 1988 Aug 2.9
6 Northern Territory Clothing, footwear and perso… A3349767W 1988 Sep 3
Why is multiplicative seasonality necessary for this series?
# let's plot it firstmyseries |>autoplot(Turnover) +labs(title ="Turnover: NSW Cafes, Restaurants, and Catering Services",y ="Turnover ($M AUD)",x ="Month" )
We observe that the annual range (peak minus trough) was about less than 2 in the late 1980s, and around 10 in 2018, while the mean rose from about 3 to about 14. In that regard, Multiplicative Seasonality is necessary because the seasonal swing grows with the level meaning they are proportional to the trend level rather than constant in absolute dollars.
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
# Fit the modelsfit_Holt_Winters <- myseries |>model(`Holt-Winters (M,A,M)`=ETS(Turnover ~error("M") +trend("A") +season("M")),`Damped HW (M,Ad,M)`=ETS(Turnover ~error("M") +trend("Ad") +season("M")) )# Inspect estimated parameters and AICcreport(fit_Holt_Winters)
Warning in report.mdl_df(fit_Holt_Winters): 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.
The preferred method will be the one with the lower RMSE, provided its residuals are also reasonably well behaved. Therefore the the standard Holt-Winters is good.
Check that the residuals from the best method look like white noise.
After identifying the model with the lower RMSE, inspect its residuals using Residual diagnostics for ordinary Holt-Winters and Ljung-Box test for residual autocorrelation.
# Ljung-Box test for white noise# For monthly seasonal data, lag = 24 is recommended (dof = number of estimated states/parameters)best_fit |>augment() |>features(.innov, ljung_box, lag =24, dof =16)
Residuals are not centered or spread evenly over time, consistent with errors that scale with the level.
Most autocorrelation falls within the dashed line bound except for 3 lags.
The histogram
Now find the test set RMSE, while training the model to the end of 2010.
# Let's create the training and test splitstrain_data <- myseries |>filter(Month <=yearmonth("2010 Dec"))test_data <- myseries |>filter(Month >yearmonth("2010 Dec"))# Let's fit both models on training data firstfit_train <- train_data |>model(`Holt-Winters (M,A,M)`=ETS(Turnover ~error("M") +trend("A") +season("M")),`Damped HW (M,Ad,M)`=ETS(Turnover ~error("M") +trend("Ad") +season("M")) )# Let's Forecast over the test period (from 2011 Jan to end of series)test_length <-nrow(test_data)fc_test <- fit_train |>forecast(h = test_length)# Let's Evaluate out-of-sample test set accuracyfc_test |>accuracy(myseries) |>select(.model, RMSE, MAE, MAPE)
The test set RMSE evaluates how well the models generalize to unseen future observations over multiple forecast horizons. The model with the smallest test RMSE forecasts this particular series best over the five-year test period.
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?
# Let's get the same seriesset.seed(12345678)myseries <- aus_retail |>filter(`Series ID`==sample(aus_retail$`Series ID`, 1))# Let's train/test split (training to end of 2010)train_data <- myseries |>filter(Month <=yearmonth("2010 Dec"))test_data <- myseries |>filter(Month >yearmonth("2010 Dec"))# Let's Fit models: Previous benchmarks + Decomposition Modelfit_models <- train_data |>model(`Holt-Winters (M,A,M)`=ETS(Turnover ~error("M") +trend("A") +season("M")),`Damped HW (M,Ad,M)`=ETS(Turnover ~error("M") +trend("Ad") +season("M")),# STL decomposition on Box-Cox transformed data:`STL + ETS (Box-Cox)`=decomposition_model(STL(box_cox(Turnover, guerrero(Turnover)) ~season(window ="periodic")),ETS(season_adjust ~season("N"))) )#Let's Generate forecasts for the test horizonfc_test <- fit_models |>forecast(new_data = test_data)# Let's Evaluate out-of-sample accuracy on the test setfc_accuracy <- fc_test |>accuracy(myseries) |>select(.model, RMSE, MAE, MAPE, MASE)fc_accuracy
The Holt-Winter has the highest test RMSE, followed by Damped HW and finally STL + ETS (Box-Cox), and its flexible seasonal decomposition and Box–Cox transformation have helped. The Guerrero Box-Cox transformation linearizes the exponential drift and stabilizes the variation when seasonal swings become larger as turnover increases, allowing the non-seasonal ETS model on season_adjust to project a more disciplined, realistic trajectory.
In conclusion, STL + ETS is a useful alternative because it estimates seasonality separately and can accommodate evolving seasonal patterns when you use a non-periodic STL window.
