Exponential Smoothing
Josef
10/9/2021
R Markdown
library(tsibbledata)
library(tsibble)##
## Attaching package: 'tsibble'## The following objects are masked from 'package:base':
##
## intersect, setdiff, unionlibrary(fable)## Loading required package: fabletoolslibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.1 ✓ dplyr 1.0.6
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()8.1. 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. b) Find the optimal values of α and ℓ and generate forecasts for the next four months. c) Compute a 95% prediction interval for the first forecast using y ± 1.96s where s is the standard deviation of the residuals. d) Compare your interval with the interval produced by R.
aus_livestock <- aus_livestock
glimpse(aus_livestock)## Rows: 29,364
## Columns: 4
## Key: Animal, State [54]
## $ Month <mth> 1976 Jul, 1976 Aug, 1976 Sep, 1976 Oct, 1976 Nov, 1976 Dec, 197…
## $ Animal <fct> "Bulls, bullocks and steers", "Bulls, bullocks and steers", "Bu…
## $ State <fct> Australian Capital Territory, Australian Capital Territory, Aus…
## $ Count <dbl> 2300, 2100, 2100, 1900, 2100, 1800, 1800, 1900, 2700, 2300, 250…unique(aus_livestock$Animal)## [1] Bulls, bullocks and steers Calves
## [3] Cattle (excl. calves) Cows and heifers
## [5] Lambs Pigs
## [7] Sheep
## 7 Levels: Bulls, bullocks and steers Calves ... Sheepfit <- aus_livestock %>%
filter(Animal == "Pigs") %>%
filter(State == "Victoria") %>%
model(ETS(Count ~ error("A") + trend("N") + season("N")))
fit## # A mable: 1 x 3
## # Key: Animal, State [1]
## Animal State `ETS(Count ~ error("A") + trend("N") + season("N"))`
## <fct> <fct> <model>
## 1 Pigs Victoria <ETS(A,N,N)>fc <- fit %>%
forecast(h = 4)
fit %>% report()## 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 report function tells us that alpha = 0.3579401 and gamma = 0.0001000139
fc %>%
filter(Animal == "Pigs") %>%
filter(State == "Victoria") %>%
autoplot(aus_livestock) +
geom_line(aes(y = .fitted), data = augment(fit)) +
labs(title = "Victoria", subtitle = "Pigs") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
glimpse(fc)## Rows: 4
## Columns: 6
## Key: Animal, State, .model [1]
## $ Animal <fct> "Pigs", "Pigs", "Pigs", "Pigs"
## $ State <fct> Victoria, Victoria, Victoria, Victoria
## $ .model <chr> "ETS(Count ~ error(\"A\") + trend(\"N\") + season(\"N\"))", "ET…
## $ Month <mth> 2019 Jan, 2019 Feb, 2019 Mar, 2019 Apr
## $ Count <dist> [<dist_normal[1]>], [<dist_normal[1]>], [<dist_normal[1]>], [<d…
## $ .mean <dbl> 95186.56, 95186.56, 95186.56, 95186.56fc %>%
filter(Animal == "Pigs") %>%
filter(State == "Victoria") %>%
autoplot(filter(aus_livestock, Month >= tsibble::yearmonth("2018 Jan"))) +
labs(title = "Victoria", subtitle = "Pigs") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
s <- augment(fit) %>% pull(.resid) %>%
sd()
s## [1] 9344.666yhat <- fc %>%
pull(Count) %>%
head(1)
yhat + c(-1, 1) * 1.96 * augment(fit) %>% pull(.resid) %>%
sd()## <distribution[2]>
## [1] N(76871, 8.7e+07) N(113502, 8.7e+07)76871## [1] 76871113502## [1] 113502fc %>%
mutate(interval = hilo(Count, 95)) %>%
pull(interval)## <hilo[4]>
## [1] [76854.79, 113518.3]95 [75927.17, 114445.9]95 [75042.22, 115330.9]95
## [4] [74194.54, 116178.6]9576854.79## [1] 76854.79113518.3## [1] 113518.38.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. b) Use an ETS(A,N,N) model to forecast the series, and plot the forecasts. c) Compute the RMSE values for the training data. 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.) e) Discuss the merits of the two forecasting methods for this data set. f) Compare the forecasts from both methods. Which do you think is best? g) Calculate a 95% prediction interval for the first forecast for each model, using the RMSE values and assuming normal errors. h) Compare your intervals with those produced using R.
Afghanistan <- global_economy %>%
filter(Country == "Afghanistan") %>%
filter(Year >= 2002)Afghanistan %>%
autoplot(Exports) +
labs(title = "Afghanistan", subtitle = "Exports") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
fit <- Afghanistan %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h = 5)
fc %>%
autoplot(Afghanistan) +
geom_line(aes(y = .fitted), color = 'red', data = augment(fit)) +
labs(title = "Afghanistan", subtitle = "Exports") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "red line = fitted value forecast")
fit %>% accuracy()## # A tibble: 1 x 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Afghani… "ETS(Exports … Trai… -2.39 5.64 4.78 -20.4 27.5 0.974 0.932 -0.263# RMSE = 3.982656fit <- Afghanistan %>%
model(
ANN_model = ETS(Exports ~ error("A") + trend("N") + season("N")),
AAN_model = ETS(Exports ~ error("A") + trend("A") + season("N")))
accuracy(fit) ## # A tibble: 2 x 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Afghanist… ANN_mod… Train… -2.39 5.64 4.78 -20.4 27.5 0.974 0.932 -0.263
## 2 Afghanist… AAN_mod… Train… -0.115 3.96 3.05 -6.01 17.1 0.622 0.653 -0.0298# 5.643453
# vs.
# 3.955708
# Discuss the merits of the two forecasting methods for this data set.# Compare the forecasts from both methods. Which do you think is best?
fit %>%
forecast(h = 4) %>%
autoplot(Afghanistan, level = NULL) +
labs(title = "Afghanistan", subtitle = "Exports and Forecast") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
# The AAN model looks like a better fit to me. It better captures the downward trend in the data. aan_rmse <- fit %>%
accuracy() %>%
filter(.model == "AAN_model") %>%
dplyr::select(Country, RMSE)
fit %>%
dplyr::select(Country, ANN_model) %>%
forecast(h = 1) %>%
left_join(aan_rmse, by = 'Country') %>%
mutate(lower_bound = Exports - 1.96 * RMSE) %>%
mutate(upper_bound = Exports + 1.96 * RMSE) %>%
dplyr::select(Country, Exports, lower_bound, upper_bound)## # A fable: 1 x 5 [?]
## # Key: Country [1]
## Country Exports lower_bound upper_bound Year
## <fct> <dist> <dist> <dist> <dbl>
## 1 Afghanistan N(6.4, 36) N(-1.4, 36) N(14, 36) 2018Forecast 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.]
Chinese_GDP <- global_economy %>%
filter(Country == "China")
fit <- Chinese_GDP %>%
model(ETS(GDP ~ error() + trend() + season()))
fc <- fit %>%
forecast(h = 5)
Chinese_GDP %>%
autoplot(GDP) 
fc %>%
autoplot(Chinese_GDP) +
geom_line(aes(y = .fitted), data = augment(fit)) +
labs(title = "China", subtitle = "GDP and Forecast") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
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?
fit <- aus_production %>%
model(ETS(Gas ~ error() + trend() + season()))
fc <- fit %>%
forecast(h = 5)
aus_production %>%
autoplot(Gas) 
fc %>%
autoplot(aus_production) +
geom_line(aes(y = .fitted), data = augment(fit)) +
labs(title = "Australia", subtitle = "Gas and Forecast") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
Recall your retail time series data (from Exercise 8 in Section 2.10). a) Why is multiplicative seasonality necessary for this series? b) Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer? Check that the residuals from the best method 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?
aus_retail <- aus_retail
class(aus_retail)## [1] "tbl_ts" "tbl_df" "tbl" "data.frame"nsw_cafes <- aus_retail %>% filter(Industry == "Cafes, restaurants and catering services") %>% filter(State == "New South Wales")
nsw_cafes %>% autoplot(Turnover) +
labs(title = "New South Wales", subtitle = "Cafes, restaurants and catering services") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
# Multiplicative seasonality is best when seasonality (the seasonality within each year) changes year-over-overfit <- nsw_cafes %>%
model(
'Holt-Winters Multiplicative' = ETS(Turnover ~ error("M")+ trend("A") + season("M")),
'Holt-Winters Damped' = ETS(Turnover ~ error("M")+ trend("Ad") + season("M")),
'Seasonal_Naive' = SNAIVE(Turnover)
)
fc <- fit %>% forecast(h = '12 years')
fc %>% autoplot(nsw_cafes, level = NULL) +
labs(title = "Australia NSW", subtitle = "Cafes, restaurants and catering services") +
theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
theme(plot.subtitle = element_text(hjust = 0.5)) +
labs(caption = "")
# The Holt-Winters Damped doesn't really do any better than a Seasonal Naive forecast. The Holt-Winters Multiplicative seems to do a lot better than both. 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?
nsw_lambda <- nsw_cafes %>%
features(Turnover, features = feasts::guerrero) %>%
pull(lambda_guerrero)
nsw_boxcox <- nsw_cafes %>%
mutate(bx = box_cox(Turnover, nsw_lambda))
fit <- nsw_boxcox %>%
model(
ETS_Boxcox = ETS(bx)
)
fit <- nsw_boxcox %>%
model(
'Holt-Winters Multiplicative' = ETS(Turnover ~ error("M")+ trend("A") + season("M")),
'Holt-Winters Damped' = ETS(Turnover ~ error("M")+ trend("Ad") + season("M")),
'Seasonal_Naive' = SNAIVE(Turnover)
)
accuracy(fit)## # A tibble: 3 x 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 New S… Cafes, r… Holt-… Trai… 1.16 16.7 11.9 0.174 4.19 0.337 0.350 0.180
## 2 New S… Cafes, r… Holt-… Trai… 1.35 16.4 11.7 0.326 4.17 0.329 0.343 0.0746
## 3 New S… Cafes, r… Seaso… Trai… 19.8 47.7 35.4 6.06 11.6 1 1 0.856# it looks like the holt winters damped model has the best fit of all the box-cox transformed models.