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.knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(fpp3)
Consider the number of pigs slaughtered in Victoria, available in the
aus_livestock dataset.
pigs <- aus_livestock |>
filter(State == "Victoria", Animal == "Pigs")
pigs |> autoplot(Count) +
labs(title = "Pigs slaughtered in Victoria", y = "Count")
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 <- pigs |>
model(ses = 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.07
The optimal values are alpha = 0.3221 and l0 = 100,646.6. The small alpha means the level is updated slowly, so the forecasts rely on a long stretch of history rather than just the most recent months.
fc <- fit |> forecast(h = 4)
fc
## # 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
## # i 2 more variables: Count <dist>, .mean <dbl>
fc |> autoplot(pigs |> filter(year(Month) >= 2010)) +
labs(title = "Four-month forecasts from simple exponential smoothing", y = "Count")
All four forecasts equal 95,187, which is what simple exponential smoothing always does: the forecast is flat for every horizon.
s <- sd(augment(fit)$.resid)
s
## [1] 9344.666
y_hat <- fc$.mean[1]
c(lower = y_hat - 1.96 * s, upper = y_hat + 1.96 * s)
## lower upper
## 76871.01 113502.10
fc |> hilo(95) |> pull(`95%`) |> head(1)
## <hilo[1]>
## [1] [76854.79, 113518.3]95
The hand calculation gives [76,871, 113,502] and R
gives [76,855, 113,518]. The two are almost identical,
differing by only about 16 units at each end, which is less than 0.1% of
the forecast. The small gap comes from the fact that R uses the model’s
estimated variance, which divides by the number of observations, while
sd() divides by n - 1 and also subtracts the residual mean.
Either way the interval is very wide, roughly plus or minus 19% of the
forecast, which reflects how noisy this series is.
Data set global_economy contains the annual Exports from
many countries. Select one country to analyse. I use
Australia.
aus_exports <- global_economy |> filter(Country == "Australia")
aus_exports |> autoplot(Exports) +
labs(title = "Australian exports", y = "% of GDP")
Exports rose from about 13% of GDP in the 1960s to about 21% in 2017, so there is a clear long-run upward trend, but it is not steady. The series was flat through the 1960s, jumped in the 1970s, fell back in the late 1980s, and rose again from the mid-1990s. It peaked at 23% in 2009 and has drifted sideways since. The data is annual, so there is no seasonality, and the year-to-year movements look like irregular cycles tied to commodity prices and the exchange rate rather than a fixed pattern.
fit_ann <- aus_exports |>
model(ann = ETS(Exports ~ error("A") + trend("N") + season("N")))
report(fit_ann)
## Series: Exports
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.5659948
##
## Initial states:
## l[0]
## 12.98943
##
## sigma^2: 1.3621
##
## AIC AICc BIC
## 257.3943 257.8387 263.5756
fit_ann |> forecast(h = 10) |> autoplot(aus_exports) +
labs(title = "ETS(A,N,N) forecasts of Australian exports", y = "% of GDP")
accuracy(fit_ann) |> select(.model, RMSE, MAE, MASE)
## # A tibble: 1 x 4
## .model RMSE MAE MASE
## <chr> <dbl> <dbl> <dbl>
## 1 ann 1.15 0.914 0.928
The training RMSE is 1.147.
fit_both <- aus_exports |>
model(ann = ETS(Exports ~ error("A") + trend("N") + season("N")),
aan = ETS(Exports ~ error("A") + trend("A") + season("N")))
accuracy(fit_both) |> select(.model, RMSE, MAE, MASE)
## # A tibble: 2 x 4
## .model RMSE MAE MASE
## <chr> <dbl> <dbl> <dbl>
## 1 ann 1.15 0.914 0.928
## 2 aan 1.12 0.893 0.907
glance(fit_both) |> select(.model, sigma2, AIC, AICc, BIC)
## # A tibble: 2 x 5
## .model sigma2 AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ann 1.36 257. 258. 264.
## 2 aan 1.34 258. 259. 269.
tidy(fit_both)
## # A tibble: 6 x 4
## Country .model term estimate
## <fct> <chr> <chr> <dbl>
## 1 Australia ann alpha 0.566
## 2 Australia ann l[0] 13.0
## 3 Australia aan alpha 0.441
## 4 Australia aan beta 0.000100
## 5 Australia aan l[0] 12.7
## 6 Australia aan b[0] 0.137
The trended model fits very slightly better on the training data (RMSE 1.117 against 1.147), but it uses two extra parameters (beta and b0) to get there, and the estimated beta is tiny (0.0001), so the slope barely adapts and is essentially fixed at its initial value of 0.14 percentage points per year. The information criteria, which penalise the extra parameters, prefer the simpler model: AICc is 257.8 for ETS(A,N,N) against 259.5 for ETS(A,A,N). A 3% improvement in RMSE is not worth two extra parameters here, especially since the upward movement in this series comes in irregular bursts rather than as a constant trend.
fit_both |> forecast(h = 10) |> autoplot(aus_exports, level = NULL) +
labs(title = "ETS(A,N,N) and ETS(A,A,N) forecasts", y = "% of GDP")
ETS(A,N,N) forecasts a flat 20.6%, while ETS(A,A,N) adds a steady upward slope of about 0.14 percentage points per year, reaching 22.1% after ten years. The gap widens the further ahead you look, so the choice does matter at longer horizons. I prefer ETS(A,N,N). It has the better AICc, it is simpler, and extrapolating a trend ten years into the future is risky for a series driven by commodity cycles, which have gone up and down repeatedly over this period.
rmse <- accuracy(fit_both) |> select(.model, RMSE)
first_fc <- fit_both |> forecast(h = 1) |> as_tibble() |> select(.model, .mean)
left_join(first_fc, rmse, by = ".model") |>
mutate(lower = .mean - 1.96 * RMSE,
upper = .mean + 1.96 * RMSE)
## # A tibble: 2 x 5
## .model .mean RMSE lower upper
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ann 20.6 1.15 18.4 22.9
## 2 aan 20.8 1.12 18.6 23.0
fit_both |> forecast(h = 1) |> hilo(95) |> select(.model, Year, .mean, `95%`)
## # A tsibble: 2 x 4 [1Y]
## # Key: .model [2]
## .model Year .mean `95%`
## <chr> <dbl> <dbl> <hilo>
## 1 ann 2018 20.6 [18.31970, 22.89462]95
## 2 aan 2018 20.8 [18.57028, 23.10700]95
The hand-calculated intervals are [18.36, 22.86] for ETS(A,N,N) and [18.65, 23.03] for ETS(A,A,N). R gives [18.32, 22.89] and [18.57, 23.11]. The intervals match closely, with R’s being very slightly wider in both cases, because the RMSE divides the sum of squared errors by the number of observations while R’s variance estimate accounts for the parameters that were estimated.
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.
china <- global_economy |> filter(Country == "China")
china |> autoplot(GDP) +
labs(title = "Chinese GDP", y = "US$")
lambda <- china |> features(GDP, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.03446284
fit_china <- china |>
model(ets = ETS(GDP),
damped = ETS(GDP ~ trend("Ad")),
log = ETS(log(GDP)),
box_cox = ETS(box_cox(GDP, lambda)))
fit_china
## # A mable: 1 x 5
## # Key: Country [1]
## Country ets damped log box_cox
## <fct> <model> <model> <model> <model>
## 1 China <ETS(M,A,N)> <ETS(M,Ad,N)> <ETS(A,A,N)> <ETS(A,A,N)>
glance(fit_china) |> select(.model, AICc)
## # A tibble: 4 x 2
## .model AICc
## <chr> <dbl>
## 1 ets 3103.
## 2 damped 3107.
## 3 log -31.9
## 4 box_cox -137.
fit_china |>
forecast(h = 20) |>
autoplot(china, level = NULL) +
labs(title = "Twenty-year forecasts of Chinese GDP under four ETS options", y = "US$")
fit_china |>
forecast(h = 20) |>
as_tibble() |>
filter(Year %in% c(2027, 2037)) |>
mutate(forecast_trillions = round(.mean / 1e12, 1)) |>
select(.model, Year, forecast_trillions) |>
arrange(.model, Year)
## # A tibble: 8 x 3
## .model Year forecast_trillions
## <chr> <dbl> <dbl>
## 1 box_cox 2027 47.6
## 2 box_cox 2037 224.
## 3 damped 2027 18.7
## 4 damped 2037 23.9
## 5 ets 2027 19.7
## 6 ets 2037 27.3
## 7 log 2027 42.8
## 8 log 2037 172.
GDP was about $12.2 trillion in 2017, the last observation. The four options give wildly different answers for 2037:
The intuition: a transformation decides whether growth is treated as constant in dollars (no transformation) or constant in percent (log or Box-Cox), and damping decides whether growth is allowed to continue forever or is gradually choked off. The AICc values cannot be compared across transformations, since the models are fitted to different data. Judging by plausibility, the damped model is the most believable of the four. China’s growth rate has already slowed, and the untransformed exponential forecasts imply an economy many times the size of the world economy today.
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?
gas <- aus_production |> select(Gas)
gas |> autoplot(Gas) +
labs(title = "Australian gas production", y = "Petajoules")
fit_gas <- gas |>
model(auto = ETS(Gas),
damped = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
additive = ETS(Gas ~ error("A") + trend("A") + season("A")))
fit_gas
## # A mable: 1 x 3
## auto damped additive
## <model> <model> <model>
## 1 <ETS(M,A,M)> <ETS(M,Ad,M)> <ETS(A,A,A)>
report(fit_gas |> select(auto))
## 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
Why multiplicative seasonality is necessary: the seasonal swings grow as the level of the series grows. In the 1950s production was around 5 petajoules and the swings were barely visible, while by the 2000s production was around 220 and the peaks and troughs are about 30 petajoules apart. Additive seasonality assumes the seasonal effect is the same size in every year, which would make the early years far too variable and the later years far too smooth. The automatic selection confirms this by choosing ETS(M,A,M).
accuracy(fit_gas) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 3 x 5
## .model RMSE MAE MAPE MASE
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 auto 4.60 3.02 4.08 0.542
## 2 damped 4.59 3.03 4.10 0.544
## 3 additive 4.76 3.35 10.9 0.600
glance(fit_gas) |> select(.model, sigma2, AICc)
## # A tibble: 3 x 3
## .model sigma2 AICc
## <chr> <dbl> <dbl>
## 1 auto 0.00324 1682.
## 2 damped 0.00329 1685.
## 3 additive 23.6 1873.
The additive model is clearly worse on every measure (MAPE of 10.9% against 4.1%), which confirms the argument above.
fit_gas |>
select(auto, damped) |>
forecast(h = "3 years") |>
autoplot(gas |> filter(year(Quarter) >= 1990), level = NULL) +
labs(title = "ETS(M,A,M) and damped forecasts of gas production", y = "Petajoules")
Does damping improve the forecasts? Barely, and not enough to prefer it. The damped model has almost the same training RMSE (4.59 against 4.60) but a slightly worse AICc (1685 against 1682). A fairer test is time series cross-validation:
gas_cv <- gas |>
stretch_tsibble(.init = 60, .step = 1) |>
model(auto = ETS(Gas ~ error("M") + trend("A") + season("M")),
damped = ETS(Gas ~ error("M") + trend("Ad") + season("M"))) |>
forecast(h = 4) |>
accuracy(gas)
gas_cv |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
## .model RMSE MAE MAPE MASE
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 auto 7.60 5.69 4.75 1.02
## 2 damped 7.54 5.67 4.76 1.02
Over four-step-ahead forecasts the damped model is marginally better (RMSE 7.54 against 7.60), a difference of less than 1%. Damping is a sensible safeguard for long horizons, since it stops the trend growing without limit, but for this series it makes almost no practical difference.
Recall your retail time series data (from Exercise 7 in Section 2.10).
set.seed(624)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
distinct(as_tibble(myseries), State, Industry)
## # A tibble: 1 x 2
## State Industry
## <chr> <chr>
## 1 New South Wales Takeaway food services
myseries |> autoplot(Turnover) +
labs(title = "Retail turnover", y = "$ million")
The size of the seasonal swings grows with the level of the series. In the early 1980s turnover was around $80 million and the December peak stood only a few million above the rest of the year. By 2018 turnover was around $600 million and the December peak is tens of millions above the other months. The seasonal effect is roughly a constant percentage of the level, not a constant number of dollars, which is exactly what multiplicative seasonality models.
fit_hw <- myseries |>
model(hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
tidy(fit_hw) |> filter(term %in% c("alpha", "beta", "gamma", "phi"))
## # A tibble: 7 x 5
## State Industry .model term estimate
## <chr> <chr> <chr> <chr> <dbl>
## 1 New South Wales Takeaway food services hw alpha 0.715
## 2 New South Wales Takeaway food services hw beta 0.0309
## 3 New South Wales Takeaway food services hw gamma 0.106
## 4 New South Wales Takeaway food services hw_damped alpha 0.775
## 5 New South Wales Takeaway food services hw_damped beta 0.00752
## 6 New South Wales Takeaway food services hw_damped gamma 0.0968
## 7 New South Wales Takeaway food services hw_damped phi 0.980
fit_hw |> forecast(h = "3 years") |> autoplot(myseries, level = NULL) +
labs(title = "Holt-Winters' multiplicative forecasts, with and without damping",
y = "$ million")
The damping parameter is estimated at phi = 0.98, which is close to 1, so the damped forecasts bend only slightly below the undamped ones over three years.
accuracy(fit_hw) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 2 x 5
## .model RMSE MAE MAPE MASE
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 hw 9.74 7.13 3.08 0.308
## 2 hw_damped 9.73 7.10 3.03 0.307
glance(fit_hw) |> select(.model, sigma2, AICc)
## # A tibble: 2 x 3
## .model sigma2 AICc
## <chr> <dbl> <dbl>
## 1 hw 0.00174 4636.
## 2 hw_damped 0.00171 4627.
The one-step RMSE values are effectively tied: 9.74 for Holt-Winters and 9.73 for the damped version. The AICc separates them more clearly, favouring the damped model (4627 against 4636), so I prefer the damped version. It fits just as well, it is more cautious at longer horizons, and the information criterion supports it.
best <- myseries |>
model(ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
best |> gg_tsresiduals()
best |> augment() |> features(.innov, ljung_box, lag = 24)
## # A tibble: 1 x 5
## State Industry .model lb_stat lb_pvalue
## <chr> <chr> <chr> <dbl> <dbl>
## 1 New South Wales Takeaway food services "ETS(Turnover ~ erro~ 25.3 0.392
The residuals do look like white noise. Almost every ACF spike sits inside the significance bounds, and the Ljung-Box test gives a p-value of 0.39, so there is no evidence of leftover autocorrelation. The residuals are centred on zero with roughly constant variance, and the histogram is close to normal with slightly heavy tails. This is a big improvement on the seasonal naive model from Homework 3, whose residuals failed the same test badly.
train <- myseries |> filter(year(Month) < 2011)
fits <- train |>
model(snaive = SNAIVE(Turnover),
hw = ETS(Turnover ~ error("M") + trend("A") + season("M")),
hw_damped = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fc <- fits |> forecast(h = "8 years")
fc |> autoplot(myseries, level = NULL) +
labs(title = "Forecasts from the 2010 training set against the actual data",
y = "$ million")
accuracy(fc, myseries) |> select(.model, RMSE, MAE, MAPE, MASE)
## # A tibble: 3 x 5
## .model RMSE MAE MAPE MASE
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 hw 70.4 60.9 15.3 3.17
## 2 hw_damped 69.0 62.8 14.7 3.27
## 3 snaive 96.8 79.5 16.3 4.14
Yes, comfortably. The seasonal naive method has a test RMSE of 96.8, while Holt-Winters gets 70.4 and the damped version 69.0. That is an improvement of about 29%. The reason is visible in the plot: the seasonal naive forecasts simply repeat 2010 and fall further behind every year, while the Holt-Winters models keep following the upward trend.
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 <- train |> features(Turnover, features = guerrero) |> pull(lambda_guerrero)
lambda
## [1] -0.1011142
fits_stl <- train |>
model(
stl_ets = decomposition_model(
STL(box_cox(Turnover, lambda)),
ETS(season_adjust ~ season("N"))),
stl_ets_damped = decomposition_model(
STL(box_cox(Turnover, lambda)),
ETS(season_adjust ~ error("A") + trend("Ad") + season("N")))
)
fc_stl <- fits_stl |> forecast(h = "8 years")
fc_stl |> autoplot(myseries, level = NULL) +
labs(title = "STL plus ETS forecasts on the Box-Cox transformed series",
y = "$ million")
bind_rows(accuracy(fc, myseries), accuracy(fc_stl, myseries)) |>
select(.model, RMSE, MAE, MAPE, MASE) |>
arrange(RMSE)
## # A tibble: 5 x 5
## .model RMSE MAE MAPE MASE
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 hw_damped 69.0 62.8 14.7 3.27
## 2 hw 70.4 60.9 15.3 3.17
## 3 stl_ets_damped 74.8 66.6 14.6 3.47
## 4 snaive 96.8 79.5 16.3 4.14
## 5 stl_ets 110. 106. 24.1 5.52
bind_rows(accuracy(fits), accuracy(fits_stl)) |>
select(.model, RMSE, MASE) |>
arrange(RMSE)
## # A tibble: 5 x 3
## .model RMSE MASE
## <chr> <dbl> <dbl>
## 1 stl_ets 7.23 0.274
## 2 stl_ets_damped 7.30 0.276
## 3 hw 8.33 0.310
## 4 hw_damped 8.34 0.311
## 5 snaive 26.1 1
The Guerrero method gives lambda = -0.10, close to a log transformation, which again reflects the seasonal variation growing with the level.
On the training data the STL approach fits best of all (RMSE 7.23 against 8.33 for Holt-Winters), because the STL decomposition can follow a seasonal pattern that changes shape over time, which the Holt-Winters seasonal component cannot.
On the test set it is a different story. Plain STL plus ETS gives an RMSE of 109.9, which is worse even than the seasonal naive benchmark, while the damped version gives 74.8. Both are beaten by the Holt-Winters models at 70.4 and 69.0.
So the STL approach does not beat my best previous forecasts. The explanation is visible in the forecast plot: ETS applied to the seasonally adjusted series picks up the steep growth of the late 2000s and extends it in a straight line on the transformed scale, which back-transforms into forecasts that climb too quickly. Damping the trend fixes most of that problem but still does not match Holt-Winters. This is also a good reminder that the best in-sample fit is not the best forecaster: the STL model was first on the training data and last on the test data.