library(tidyverse)
library(fpp3)
library(fable)Assignment 5: Exponential Smoothing
Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. Link to Chapter 8 exercises for reference.
8.1
Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
victoria_pigs <- aus_livestock |>
filter(Animal == 'Pigs' & State == 'Victoria')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.
To estimate the equivalent model for simple exponential smoothing, we will need to specify in the ETS function error("A") and trend("N") and season("N"). These are the inputs to produce a simple exponential smoothing model.
victoria_pigs_fit <- victoria_pigs |>
model(ETS(Count ~ error("A") + trend("N") + season("N")))
report(victoria_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.07 Based on the above, we can see that using simple exponential smoothing, the optimal value of α = 0.3221247. We can also see that the optimal value of ℓ0 = 100646.6.
Now we will generate forecasts for the next four months.
victoria_pigs_forecast <- victoria_pigs_fit |>
forecast(h = 4)
victoria_pigs_forecast |>
autoplot(victoria_pigs) +
labs(title = "Pigs Slaughtered in Victoria (1972-2018) with 4-Month Forecast
using Simple Exponential Smoothing (α = 0.3221247, ℓ0 = 100646.6)")
b. Compute a 95% prediction interval for the first forecast using y-hat ± 1.96 where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
First, we’ll identify the standard deviation of the residuals and multiply that by ±1.96 to get the upper and lower bounds of the 95% confidence interval. Then, we’ll compare that to the value produced in the forecast using R.
victoria_pigs_standarddeviation <- residuals(victoria_pigs_fit)$.resid |> sd()
victoria_pigs_standarddeviation[1] 9344.666victoria_pigs_yhat <- victoria_pigs_forecast$.mean[1]
victoria_pigs_yhat[1] 95186.56From the above we see that the standard deviation is 9344.666, and so we’ll multiply the standard deviation by 1.96 – we’ll subtract that value from 95186.56 to get the lower bound and add it to get the upper bound.
victoria_pigs_lower <-
victoria_pigs_yhat - (victoria_pigs_standarddeviation * 1.96)
victoria_pigs_upper <-
victoria_pigs_yhat + (victoria_pigs_standarddeviation * 1.96)
victoria_pigs_lower[1] 76871.01victoria_pigs_upper[1] 113502.1The computed value of the 95% prediction interval for the first forecast using y-hat ± 1.96 where s is the standard deviation of the residuals is from 76871.01 to 113502.1.
The value of the 95% confidence interval produced using R:
victoria_pigs_forecast |>
hilo(95) |>
pull('95%') |>
head(1)<hilo[1]>
[1] [76854.79, 113518.3]95We can see that the value of the 95% confidence interval for the first period forecast using R is from 76854.79 to 113518.3.
The interval produced using R (76854.79, 113518.3) is slightly broader (wider) but generally very similar and close to the computed interval (76871.01, 113502.1).
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.
Based on running ?global_economy, we know that Exports are the exports or goods and services as a percentage (%) of GDP. The global_economy dataset is a tsibble that details economic indicators featured by the World Bank from 1960 to 2017.
For this exercise, I selected Denmark and plotted its exports over time.
global_economy |>
filter(Country == "Denmark") |>
autoplot(Exports) +
labs(title = "Annual Exports from Denmark as a Percentage of GDP",
y = "Exports (% of GDP)")
Several key features emerge. The first visible feature is long-term growth. Denmark’s exports as a percentage of GDP show a clear upward trend from 1960 to 2017. Initially, exports hovered around 30% of GDP, but they steadily increased, reaching approximately 55% in the most recent years of the dataset. We also see variability (fluctuations) in the overall upward growth trend, including a period of decline in the mid-to-late 1970s, and increase in the late 1970s througr the early 1980s, and then another downturn in the mid-1980s. We can also see the impact of the global financial crisis (the Great Recession) around 2008, with a sharp drop in exports as a percentage of GDP, from nearly 55% to less than 50%. Exports have increased since that period of decline, with recent relative stability around 55% for the last 5-6 years of the dataset. It’s unclear if the visible fluctuations – periods of growth and decline, and overall trend – are a reflection of overall global economic trends, trade policies and agreements, increases in productivity or global demand for a certain export, or something else.
Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
denmark_exports_fit <- global_economy |>
filter(Country == "Denmark") |>
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
report(denmark_exports_fit)Series: Exports
Model: ETS(A,N,N)
Smoothing parameters:
alpha = 0.9998996
Initial states:
l[0]
32.31581
sigma^2: 4.1743
AIC AICc BIC
322.3492 322.7936 328.5305 denmark_exports_forecast <- denmark_exports_fit |>
forecast(h = 5)
denmark_exports_forecast |>
autoplot(global_economy) +
labs(title = "Annual Exports from Denmark as a Percentage of GDP
with a 5-year forecast using a ETS(A,N,N) Model",
y = "Exports (% of GDP)")
Compute the RMSE values for the training data.
We’ll compute the RMSE (root mean squared error) for the training data by calculating the residuals (difference between actual and fitted values) and then calculating the RMSE. We’ll use the accuracy() function from fable to do this.
denmark_exports_accuracy <- denmark_exports_fit |>
accuracy()
denmark_exports_accuracy$RMSE[1] 2.007573The RMSE for the ETS(A)
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.
denmark_exports_fit_AAN <- global_economy |>
filter(Country == "Denmark") |>
model(ETS(Exports ~ error("A") + trend("A") + season("N")))
report(denmark_exports_fit_AAN)Series: Exports
Model: ETS(A,A,N)
Smoothing parameters:
alpha = 0.9998998
beta = 0.0001000284
Initial states:
l[0] b[0]
32.46829 0.4192062
sigma^2: 4.1654
AIC AICc BIC
324.1161 325.2699 334.4183 denmark_exports_forecast_AAN <- denmark_exports_fit_AAN |>
forecast(h = 5)
denmark_exports_forecast_AAN |>
autoplot(global_economy) +
labs(title = "Annual Exports from Denmark as a Percentage of GDP
with a 5-year forecast using an ETS(A,A,N) Model",
y = "Exports (% of GDP)")
The previous ETS(A,N,N) model assumed the data followed a level pattern with no trend. The forecast is relatively stable. When we compare the above ETS(A,A,N) model we can see the notable difference in using a model that does include a trend component and assumes exports will continue with this growth trend at a steady rate.
denmark_exports_accuracy_AAN <- denmark_exports_fit_AAN |>
accuracy()
denmark_exports_accuracy_AAN$RMSE[1] 1.969295Compare the forecasts from both methods. Which do you think is best?
The ETS(A,A,N) model forecast appears to be a better representation of the existing trend that we observe in the training data, as it captures the trend component whereas ETS(A,N,N) does not capture any trend. The ETS(A,N,N) is a more conservative forecast in that it does not predict any future growth and is more based on the last observed level. The ETS(A,A,N) uses the additional parameter for the trend, which appears to improve accuracy in that it captures the visible growth trend over time. I think that the ETS(A,A,N) is better if we are comparing the two, and assuming that the growth trend will continue over time as it has over the last nearly 60 years.
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.
lower_denmark_exports_ANN <- denmark_exports_forecast$.mean[1] - (1.96 * denmark_exports_accuracy$RMSE)
upper_denmark_exports_ANN <- denmark_exports_forecast$.mean[1] + (1.96 * denmark_exports_accuracy$RMSE)
lower_denmark_exports_AAN <- denmark_exports_forecast_AAN$.mean[1] - (1.96 * denmark_exports_accuracy_AAN$RMSE)
upper_denmark_exports_AAN <- denmark_exports_forecast_AAN$.mean[1] + (1.96 * denmark_exports_accuracy_AAN$RMSE)lower_denmark_exports_ANN[1] 51.27509upper_denmark_exports_ANN[1] 59.14477lower_denmark_exports_AAN[1] 51.76921upper_denmark_exports_AAN[1] 59.48885The 95% prediction interval for the first forecast for each model using the RMSE values and assuming normal errors are: (51.275, 59.145) for the ETS(A,N,N) model (51.769, 59.489) for the ETS(A,A,N) model
The intervals are almost the same for the first forecast. The interval for the ETS(A,N,N) model is slightly wider than the ETS(A,A,N) model. The difference for the ETS(A,N,N) model is 7.86968 (59.14477 - 51.27509 = 7.86968) and the difference for the ETS(A,A,N) model is 7.71964 (59.48885 - 51.76921 = 7.71964).
Comparing each to the intervals produced using R:
ETS(A,N,N) model:
denmark_exports_forecast |>
hilo(95) |>
pull('95%') |>
head(1)<hilo[1]>
[1] [51.20551, 59.21435]95This is very similar though not exactly the same (as expected) to (51.275, 59.145).
ETS(A,A,N) model:
denmark_exports_forecast_AAN |>
hilo(95) |>
pull('95%') |>
head(1)<hilo[1]>
[1] [51.62888, 59.62918]95This is very similar though not exactly the same (as expected) to (51.769, 59.489)
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.]
global_economy |>
filter(Country == "China") |>
autoplot(GDP) +
labs(title = "Annual GDP from China",
y = "Gross domestic product (in $USD February 2019)")
Wow, GDP growth is extremely strong here since 2000 onward. It will be interesting to use and compare various models. Let’s get into it and see what different ETS models do.
We’ll start with an ETS(A,A,N) model (Holt’s method) because that seemed to represent the data better than a ETS(A,N,N) model in the previous exercise.
china_gdp_AAN_fit <- global_economy |>
filter(Country == "China") |>
model(ETS(GDP ~ error("A") + trend("A") + season("N")))
report(china_gdp_AAN_fit)Series: GDP
Model: ETS(A,A,N)
Smoothing parameters:
alpha = 0.9998964
beta = 0.5518569
Initial states:
l[0] b[0]
50284778074 3288256684
sigma^2: 3.87701e+22
AIC AICc BIC
3258.053 3259.207 3268.356 china_gdp_AAN_forecast <- china_gdp_AAN_fit |>
forecast(h = 5)
china_gdp_AAN_forecast |>
autoplot(global_economy) +
labs(title = "Annual GDP from China with a 5-year forecast
Using an ETS(A,N,N) Model",
y = "Gross domestic product (in $USD February 2019)")
This continued exponential growth is certainly following the trend, as expected with the ETS(A,A,N) model, but is this really the best forecast? Let’s compare it to some other options, including Damped Holt’s method. Let’s also look at the Box-Cox and Log methods because, just based on the way the data initially looks, this seems (on intuition alone) that it could be a useful comparison. And let’s increase the h to 20 years to get a broader sense of the ways these different models work.
First, if we’re using Box-Cox, let’s capture the optimal lambda using the guerrero method.
china_gdp_lambda <- global_economy |>
filter(Country == "China") |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
china_gdp_lambda[1] -0.03446284china_gdp_ETScompared_fit <- global_economy |>
filter(Country == "China") |>
model(`Holt’s method (Additive, Additive, No Seasonality)` = ETS(
GDP ~ error("A") + trend("A") + season("N")),
`Damped Holt’s Method (Additive, Damped)` = ETS(
GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Box-Cox Transformation` = ETS(box_cox(
GDP,china_gdp_lambda) ~ error("A") + trend("Ad") + season("N")),
`Damped Box-Cox Transformation` = ETS(box_cox(
GDP,china_gdp_lambda)~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Log Transformation` = ETS(log(
GDP) ~ error("A") + trend("A") + season("N")),
`Damped Log Transformation` = ETS(log(
GDP) ~ error("A") + trend("Ad", phi = 0.8) + season("N"))
)
china_gdp_ETScompared_forecast <- china_gdp_ETScompared_fit |>
forecast(h = 20)
china_gdp_ETScompared_forecast |>
autoplot(global_economy, level = NULL) +
labs(title="China's GDP: ETS models compared") +
guides(colour = guide_legend(title = "Forecast Model"))
china_gdp_ETScompared_accuracy <- china_gdp_ETScompared_fit |>
accuracy()
china_gdp_ETScompared_accuracy# A tibble: 6 × 11
Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 China Holt’… Trai… 2.36e10 1.90e11 9.59e10 1.41 7.62 0.442 0.453 0.00905
2 China Dampe… Trai… 5.94e10 2.07e11 9.83e10 2.55 8.03 0.454 0.494 -0.0283
3 China Box-C… Trai… 6.35e 9 1.96e11 1.02e11 1.77 7.26 0.472 0.468 0.135
4 China Dampe… Trai… 4.21e10 1.94e11 9.43e10 2.61 7.23 0.435 0.463 -0.0964
5 China Log T… Trai… -3.48e10 2.86e11 1.25e11 0.744 7.21 0.577 0.683 0.654
6 China Dampe… Trai… 4.23e10 1.95e11 9.43e10 2.60 7.23 0.435 0.464 -0.0987 This is fascinating and so informative to see the different models over h = 20 years.
I selected phi=0.8 for the damped models as it is a commonly used starting point in time series forecasting. This allows the trend to persist while gradually decaying over time, which aligns with the historical GDP data. Given China’s rapid growth up to the dataset’s endpoint, an undamped model would assume continued exponential expansion, which is unlikely in the presence of known/real economic cycles, policy changes, and global downturns. The damped trend prevents over-extrapolation while still capturing the momentum of past growth.
The model with the lowest RMSE is Holt’s method (A,A,N), which means it fits the training data more closely - it closely follows past GDP growth. While this means it fits the data trends in the dataset most accurately of the various models in comparison, this doesn’t necessarily mean it’s the best forecast, especially because it means that we expet everything to continue exactly as it’s been. We know that exponential growth rarely continues unchecked in real life - for instance, the global COVID-19 pandemic and the economic ramifications globally, or the 2008 global recession, or many other factors. So, while Holt’s has the lowest RMSE out of all those in this model comparison, for other reasons we may prefer to use a different model like Damped Holt’s, which is more cautious and often better for long-term forecasting, making it more reasonable choice here. More on this below.
Model comparison: Holt’s Method (AAN): forecast behavior shows continued strong trend growth Damped Holt’s Method: forecast behavior is more conservative, growth trend slows over time Box-Cox transformation: Adjusts for variance instability, slope still steep with growth trend Damped Box-Cox: More “reasonable” than Box-Cox alone, but still shows notable growth trend Log transformation: Transforms GDP into a more stable pattern but as a result appears to exaggerate future growth forecastby convering the exponential growth into a linear trend Damped Log transformation: More convervative in growth forecast than Log transformation, appears more realistic as a result
We see that the Damped Box-Cox, the Damped Log Transformation, and Holt’s method are all very close in proximity in the above graph. Let’s take the Log and Box-Cox out of the plot to compare the others more closely.
The below tells us a lot about how important context is and how important it is to compare results at different levels. If we were looking at a 3-to-5-year forecast, all these 4 models would look roughly the same. If shorter-term (3-5-year) forecasting is the goal, Holt’s method might be reasonable since trends often persist in the short term. If long-term forecasting (20+ years) is needed, Damped Holt’s would be better: the Damped Holt’s method balances capturing the historical trend without over-extrapolating unrealistic exponential growth.
china_gdp_ETScompared_fit_fewer <- global_economy |>
filter(Country == "China") |>
model(`Holt’s method (Additive, Additive, No Seasonality)` = ETS(
GDP ~ error("A") + trend("A") + season("N")),
`Damped Holt’s Method (Additive, Damped)` = ETS(
GDP ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Damped Box-Cox Transformation` = ETS(box_cox(
GDP,china_gdp_lambda)~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Damped Log Transformation` = ETS(log(
GDP) ~ error("A") + trend("Ad", phi = 0.8) + season("N"))
)
china_gdp_ETScompared_forecast_fewer <- china_gdp_ETScompared_fit_fewer |>
forecast(h = 20)
china_gdp_ETScompared_forecast_fewer |>
autoplot(global_economy, level = NULL) +
labs(title="China's GDP: ETS models compared") +
guides(colour = guide_legend(title = "Forecast Model"))
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)
Looking at the data prior to any modeling or forecasting, we are reminded that the Gas data in aus_production shows increasing variance in the seasonal trend. There is a visible seasonality and its trend increases over time (it is not a stable seasonality where the same peaks and troughs recur over time – the variance increases in the seasonality). For this reason a multiplicative seasonality of the model will be necessary. Let’s see which model is determined to be optimal:
aus_production |>
model(ETS(Gas)) |>
report()Series: Gas
Model: ETS(M,A,M)
Smoothing parameters:
alpha = 0.6528545
beta = 0.1441675
gamma = 0.09784922
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
sigma^2: 0.0032
AIC AICc BIC
1680.929 1681.794 1711.389 The model chosen here is the ETS(M,A,M) model – selected because it has the smallest AICc value, whcih as we learned in the book means it’s generally the model with the best forecast.
Now, knowing the seasonality is quarterly, in order to see the next few years of data we’ll set h = 12 (3 years of 4 seasonal cycles).
aus_production |>
model(ETS(Gas)) |>
forecast(h = 12) |>
autoplot(aus_production)
That looks pretty good to me! Now, per the question prompt, let’s experiment with including an additive trend into the model: we’ll go from the ETS(M,A,M) above to the ETS(M,Ad,M) model:
aus_production |>
model(ETS(Gas ~ error("M") + trend("Ad") + season("M"))) |>
report()Series: Gas
Model: ETS(M,Ad,M)
Smoothing parameters:
alpha = 0.6489044
beta = 0.1551275
gamma = 0.09369372
phi = 0.98
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
5.858941 0.09944006 0.9281912 1.177903 1.07678 0.8171255
sigma^2: 0.0033
AIC AICc BIC
1684.028 1685.091 1717.873 We can see that the AICc is just slightly larger in the ETS(M,Ad,M) model (AICc = 1685.091) compared with the ETS(M,A,M) model (AICc = 1681.794). Since AICc favors models with a better trade-off between fit and complexity, the non-damped model is marginally preferred, but the difference is small.
aus_production |>
model(
`ETS(M,A,M)` = ETS(Gas ~ error("M") + trend("A") + season("M")),
`ETS(M,Ad,M) (Damped)` = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
) |>
forecast(h = 12) |>
autoplot(aus_production) +
labs(title = "Comparison of ETS(M,A,M) vs. ETS(M,Ad,M) Forecasts")
aus_production_models <- aus_production |>
model(
`ETS(M,A,M)` = ETS(Gas ~ error("M") + trend("A") + season("M")),
`ETS(M,Ad,M) (Damped Trend)` = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)
accuracy(aus_production_models)# 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(M,A,M) Trai… -0.115 4.60 3.02 0.199 4.08 0.542 0.606 -0.0131
2 ETS(M,Ad,M) (Dampe… Trai… -0.00439 4.59 3.03 0.326 4.10 0.544 0.606 -0.0217The forecasts look almost identical, meaning the damping adjustment has little visible effect over the 12 quarters in our forecast. The forecasts look so similar such as to appear the same in the plot above.
Interestingly the damped trend model has a lower RMSE - ever so slightly. The RMSE difference is less than .01: RMSE = 4.595113 for the ETS(M,A,M) model and RMSE = 4.591840 for the ETS(M,Ad,M) (Damped Trend) model.
Damping did not significantly change the forecasts, likely because the trend in gas production is strong and persistent. The RMSE and AICc differences are too small to justify preferring one model over the other. If we were increasing our forecast into the more-distant future (more than a few years) it would make sense to prefer the model that uses damping, which would prevent over-extrapolation of trends. But as noted above, and as visualized in the plots above and the accompanying data about the models, the ETS(M,A,M) model and ETS(M,Ad,M) (Damped Trend) model appear nearly the same in the 12-quarter forecast above.
8.8
Recall your retail time series data (from Exercise 7 in Section 2.10).
set.seed(5889)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))Autoplot of myseries for reference:
autoplot(myseries)Plot variable not specified, automatically selected `.vars = Turnover`
Why is multiplicative seasonality necessary for this series?
The multiplicative seasonality is necessary for this series because, similar to the Gas data from aus_production, the seasonal fluctuations increase in magnitude as the trend rises. This suggests that the variation in seasonality scales with the level of the series. In other words, the variation in seasonality does not remain constant over time. This is why a multiplicative seasonality component is necessary.
The necessity of the multiplicative seasonality is confirmed by the fable package: the ETS(M,A,M) model (which represents the Holt-Winters Multiplicative method) is identified as the optimal choice.
myseries |>
model(ETS(Turnover)) |>
report()Series: Turnover
Model: ETS(M,A,M)
Smoothing parameters:
alpha = 0.5375567
beta = 0.0001448645
gamma = 0.0001523882
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
8.841997 0.09803851 0.9814989 0.9005267 0.8926347 1.409731 1.058104 0.9952105
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
0.9734195 0.992143 0.9772829 0.9259911 0.9665959 0.9268618
sigma^2: 0.0023
AIC AICc BIC
2905.523 2906.970 2975.037 Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
The ETS(M,A,M) model is Holt-Winters’ multiplicative method (shorthand A,M for additive, multiplicative), which was selected as the optimal choice. We’ll compare it to Damped Holt-Winters’ multiplicative method. I abbreviated Holt-Winters as HW in the below.
myseries_models <- myseries |>
model(
`HW Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Damped HW Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
)myseries_models |>
forecast(h = 24) |>
autoplot(myseries) +
labs(title = "Comparison of Holt-Winters' Multiplicative vs Damped Holt-Winters' Multiplicative Models")
myseries |>
model(
`HW Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Damped HW Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M"))
) |>
report()Warning in report.mdl_df(model(myseries, `HW Multiplicative` = ETS(Turnover ~ :
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: 2 × 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 Australi… Other r… HW Mu… 0.00234 -1436. 2906. 2907. 2975. 1.91 2.53 0.0373
2 Australi… Other r… Dampe… 0.00241 -1440. 2916. 2918. 2990. 1.95 2.60 0.0376The AICc = 2906.970 for HW Multiplicative. For Damped HW Multiplicative, the AICc = 2917.895. The HW Multiplicative model has the lower AICc (preferred), which makes sense at it was selected as the optimal model. This indicates that adding a damped trend doesn’t improve the model performance.
accuracy(myseries_models)# 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 Australian … Other r… HW Mu… Trai… 0.0119 1.38 1.04 -0.195 3.74 0.452 0.469
2 Australian … Other r… Dampe… Trai… 0.117 1.39 1.04 0.269 3.73 0.453 0.473
# ℹ 1 more variable: ACF1 <dbl>For a 24 month forcast (h = 24) when comparing these models, the RMSE = 1.383601 for HW Multiplicative and for Damped HW Multiplicative, the RMSE = 1.394638. The numbers are very similar, so the accuracy of the forecasts are not vastly different. That said, the lower RMSE indicates the forecasts are more accurate. So, once again, the HW Multiplicative is preferred.
Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
First we’ll extract the fitted values for one-step forecasts and then compute the RMSE.
accuracy(myseries_models) |>
select(.model, RMSE)# A tibble: 2 × 2
.model RMSE
<chr> <dbl>
1 HW Multiplicative 1.38
2 Damped HW Multiplicative 1.39As noted previously, the non-damped model is preferred.
Check that the residuals from the best method look like white noise.
myseries |>
model(ETS(Turnover ~ error("M") + trend("A") + season("M"))) |>
gg_tsresiduals() +
ggtitle("Holt-Winters Multiplcative")
There are several instances of lags spiking (bottom left) - this is not what we’d expect white noise to look like. The histogram (bottom right) of residuals does resemble what we’d expect for white noise. The innovation residuals (top) do appear to fluctuate somewhat randomly. We’d need to do further testing to totally confirm what’s going on with the bottom left plot which indicates the residuals from the best method do not look like white noise.
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?
Since this model doesn’t pass the residual diagnostic tests, it’s not clear that we should move forward and use it to make predictions on test data. However, since that’s what the prompt is asking for, we’ll go ahead.
myseries_train <- myseries |>
filter(year(Month) <= 2010)
myseries_test <- anti_join(myseries, myseries_train)Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`myseries_train
myseries_train_fit <- myseries_train |>
model(ETS(Turnover ~ error("M") + trend("A") + season("M")))
myseries_train_forecast <- myseries_train_fit |>
forecast(new_data = myseries_test)
myseries_train_forecast |> autoplot(myseries) +
labs(title = "Holt-Winters Multiplicative Model: ETS(M,A,M)
Training Set (Up to End of 2010) and Forecast Compared to Test Data")
accuracy(myseries_train_forecast, data = myseries_test) |>
select(RMSE, MAPE)# A tibble: 1 × 2
RMSE MAPE
<dbl> <dbl>
1 6.42 14.2The results of the Seasonal Naive model used in Exercise 7 of chapter 5 had a RMSE = 7.832651, so this (RMSE = 6.41) is definitely an improvement.
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?
myseries_lambda <- myseries |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
myseries_lambda[1] 0.08542275myseries$Turnover_BoxCox <- box_cox(myseries$Turnover, myseries_lambda)myseries_STL <- myseries |>
model(STL(Turnover_BoxCox ~season(window = "periodic")))
myseries_components <- components(myseries_STL)
myseries$Turnover_BoxCox_STL <- myseries_components$trend + myseries_components$season_adjustmyseries_stl_train <- myseries |>
filter(year(Month) <= 2010)
myseries_stl_test <- anti_join(myseries, myseries_stl_train)Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover,
Turnover_BoxCox, Turnover_BoxCox_STL)`myseries_stl_fit <- myseries_stl_train |>
model(ETS(Turnover_BoxCox_STL ~ error("M") + trend('A') + season("M")))
myseries_stl_forecast <- myseries_stl_fit |>
forecast(new_data = myseries_stl_test)accuracy(myseries_stl_forecast, data = myseries_stl_test) |>
select(RMSE, MAPE)# A tibble: 1 × 2
RMSE MAPE
<dbl> <dbl>
1 0.368 3.76The above exercise revealed an RMSE of 6.419511, so the current RMSE of 0.3682462 represents a significant improvement. However, since the data has been transformed, RMSE may not be the best metric for direct comparison.
Looking at MAPE, the result we saw in the final section of 8.8 (above) had a MAPE of 14.21092, whereas the current MAPE is 3.761589, which represents a notable reduction (a significant improvement). This supports the conclusion of improved overall performance using an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data.