HW5
Gehad Gad
3/13/2022
Exponential smoothing
#Libraries required
#install.packages("forecast")
library(forecast)
library(fpp3)Exercise 8.1
Consider the 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.
#Examine the data
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 1800summary(aus_livestock$Animal)## Bulls, bullocks and steers Calves
## 4056 4464
## Cattle (excl. calves) Cows and heifers
## 4464 3546
## Lambs Pigs
## 3906 4464
## Sheep
## 4464summary(aus_livestock$State)## Australian Capital Territory New South Wales
## 3300 3810
## Northern Territory Queensland
## 3204 3810
## South Australia Tasmania
## 3810 3810
## Victoria Western Australia
## 3810 3810We need to filter the aus_livestock dataset for the animal to (Pigs), the state to (Victoria), and the months from January of 2005 to January of 2018.
Pigs <- aus_livestock %>%
filter(Animal == "Pigs" & State == "Victoria") %>%
filter(Month >= yearmonth(as.Date("2005-01-01")) & Month <= yearmonth(as.Date("2018-08-01")))autoplot(Pigs)+
labs(title = "Slaughter of Victorian “Pigs”")+
theme_replace()## Plot variable not specified, automatically selected `.vars = Count`
#Generate ETS model.
#ANN model is simple exponential smoothing with additive errors is given by the ets model Methods
Pigs$Count %>%
ets(model = "ANN")## ETS(A,N,N)
##
## Call:
## ets(y = ., model = "ANN")
##
## Smoothing parameters:
## alpha = 0.2904
##
## Initial states:
## l = 65501.6689
##
## sigma: 6956.338
##
## AIC AICc BIC
## 3742.316 3742.466 3751.615α = 0.2904 and ℓ = 65501.6689.
Generate forecasts for 4 months and plot the result:
#generate forecasts for 4 months
Pigs$Count %>%
ets(model = "ANN") %>%
forecast(h = 4) -> Pigs_forecast
#plot the result
Pigs_forecast %>% autoplot()
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.
Pigs_forecast## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 165 97932.55 89017.65 106847.5 84298.38 111566.7
## 166 97932.55 88649.29 107215.8 83735.03 112130.1
## 167 97932.55 88295.00 107570.1 83193.19 112671.9
## 168 97932.55 87953.28 107911.8 82670.58 113194.5We can see the Lo 95 and Hi 95 in the last two columns.
Lo 95 = 84298.38
Hi 95 = 111566.7
95% prediction interval computation
#Calculate lower 95%
s <- sd(Pigs_forecast$residuals)
Pigs_forecast$mean[1] - 1.96 * s ## [1] 84406.1#Calculate upper 95%
Pigs_forecast$mean[1] + 1.96 * s ## [1] 111459There is a small difference between the interval and the interval produced by R.
Exercise 8.5
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.
#global_economy$Country
#Filter for Exports for one country
Alb_Eco <- global_economy %>%
filter(Country == 'Albania')
#visualize the subset
Alb_Eco %>%
autoplot(Exports) +
labs(y = "Exports", title = "Albania Exports")## Warning: Removed 20 row(s) containing missing values (geom_path).
This is a yearly data where there appears to be a trend but not seasonality.
In the beginning of the series, Albania experienced a decline in Exports, then a small climb, then another fluctuating climb and decline, followed by a huge and continuous growth.
b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
#subset a training data - 2016, 2017 to be forecast
Alb_Eco_train <- Alb_Eco %>%
filter(Year <= year(as.Date("2015-01-01"))) #ETS (AAN) model - AIC 147.76
AlbaniaEST_holt <- Alb_Eco_train$Exports %>%
ets(model = "AAN") %>%
forecast(h = 2) ## Warning in ets(., model = "AAN"): Missing values encountered. Using longest
## contiguous portion of time series#calculate corresponding RMSE
sqrt(mean(AlbaniaEST_holt$residuals^2))## [1] 2.533468c. Compute the RMSE values for the training data.
AlbaniaETS_ses <- Alb_Eco_train$Exports %>%
ets(model = "ANN") %>%
forecast(h = 2) ## Warning in ets(., model = "ANN"): Missing values encountered. Using longest
## contiguous portion of time series#calculate corresponding RMSE
sqrt(mean(AlbaniaETS_ses$residuals^2))## [1] 2.533561d. 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.
The RMSE values are slightly smaller for the simple model than the holt model.
e. Compare the forecasts from both methods. Which do you think is best?
#ETS(A,A,N) Holt model:
AlbaniaEST_holt %>% autoplot()
#ETS(A,N,N) simple model:
AlbaniaETS_ses %>% autoplot()
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.
AlbaniaETS_ses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 57 27.26748 23.92646 30.60850 22.15783 32.37713
## 58 27.26748 22.54280 31.99216 20.04170 34.49326AlbaniaEST_holt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 57 27.28838 23.84466 30.73210 22.02167 32.55510
## 58 27.30928 22.43913 32.17944 19.86102 34.75754The last two columns of the first row is what we’re interested in:
Simple: Lo 95 = 22.15783, Hi 95 = 32.37713 Holt: Lo 95 = 22.02167, Hi 95 = 32.55510
Compare to the 95% prediction interval computation:
#95% prediction interval computation (simple)
s1 <- sd(AlbaniaETS_ses$residuals)
#lower 95%
AlbaniaETS_ses$mean[1] - 1.96 * s1 ## [1] 22.23458#upper 95%
AlbaniaETS_ses$mean[1] + 1.96 * s1 ## [1] 32.30038#95% prediction interval computation (holt)
s2 <- sd(AlbaniaEST_holt$residuals)
#lower 95%
AlbaniaEST_holt$mean[1] - 1.96 * s2 ## [1] 22.25488#upper 95%
AlbaniaEST_holt$mean[1] + 1.96 * s2## [1] 32.32189Simple and Holt model calculations were within 2 point of what was produced in R, the simple model was closer.
Exercise 8.6
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 to China
china_GDP <- global_economy %>%
filter(Country == 'China') %>%
mutate(GDP = GDP/1e9)
#GDP:Gross domestic product (in $USD February 2019).
china_GDP %>%
autoplot(GDP) +
labs(y="GDP in $USD", title="GDP in China")
There is a noticable upward trend to the Chinese economy.
fit <- china_GDP %>%
model(
SES = ETS(GDP ~ error("A") + trend("N") + season("N")),
Holt = ETS(GDP ~ error("A") + trend("A") + season("N")),
Damped = ETS(GDP ~ error("A") + trend("Ad") + season("N"))
)
# Forecast for 15 years
fc <- fit %>% forecast(h = 15)
fc %>%
autoplot(china_GDP, level=NULL) +
labs(y="Billions $USD", title="GDP: China") +
guides(colour = guide_legend(title = "Forecast"))
I used h = 15 (forecast for 15 years) and observed that the simple exponential smoothing model produced a horizontal line. The Holt approach produced a forecast that follows an exponential growth that clearly takes off after few years. The dampened method also starts to increase but not quite at the rate of the Holt method.
Exercise 8.7
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)
aus_gas <- aus_production %>%
select(Quarter, Gas)
fit <- aus_gas %>%
model(
mult_sea = ETS(Gas ~ error("M") + trend("A") + season("M")),
add_sea = ETS(Gas ~ error("A") + trend("A") + season("A")),
mult_sea_damp = ETS(Gas ~ error("M") + trend("Ad") + season("M"))
)# Forecast for 5 years (interval is quarters)
fc <- fit %>% forecast(h = 20)
fc %>%
autoplot(aus_gas, level=NULL) +
labs(y="Petajoules", title="Gas Production: Australia") +
guides(colour = guide_legend(title = "Forecast"))
Multiplicative seasonality is necessary because as the initial forecast plot indicates, the seasonal nature of the time series increases proportional to the level of the series itself.
There is a small difference is visible with the dampening applied in the multiplicative seasonality forecasts with and without dampening. Therefore, dampening does not improve the forecast.
Exercise 8.8
8.8 Recall your retail time series data (from Exercise 8 in Section 2.10).
set.seed(555555)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
autoplot(myseries)## Plot variable not specified, automatically selected `.vars = Turnover`
a. Why is multiplicative seasonality necessary for this series?
The retail time series above displays increasing proportional seasonality but the seasonal variation is not constant. As per the book, “multiplicative seasonality is necessary it is preferred when the seasonal variations are changing proportional to the level of the series.”
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
#create time series with Month, Turnover
timeseries <- ts(data=myseries$Turnover, frequency=12)
#apply HW multiplicative (without damping):
retail_hw <- ets(timeseries, model="MAM")
Fit_hw <- retail_hw %>% forecast(h=50)
Fit_hw %>% autoplot() 
#apply HW multiplicative (with damping):
retail_hw_damped <- ets(timeseries, model="MAM", damped = TRUE)
Fit_hw_damped <- retail_hw %>% forecast(h=50)
Fit_hw_damped %>% autoplot()
c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
sqrt(mean(Fit_hw$residuals^2))## [1] 0.09531695sqrt(mean(Fit_hw_damped$residuals^2))## [1] 0.09531695The damped model gives higher RMSE which is indicative of higher error and reduced accuracy. Therefore, the undamped model is better.
d. Check that the residuals from the best method look like white noise.
#check residuals of best method
checkresiduals(Fit_hw)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,M)
## Q* = 37.729, df = 8, p-value = 8.447e-06
##
## Model df: 16. Total lags used: 24e. 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?
#train model to the end of 2010
myseries_train <- myseries %>%
filter(year(Month) < 2011)
myseries_test <- myseries %>%
filter(year(Month) >= 2011)
#create timeseries
train <- ts(data=myseries_train$Turnover, frequency=12)
test <- ts(data=myseries_test$Turnover, frequency=12)
#create corresponding models
retail_hw_train <- ets(train, model="MAM")
retail_hw_test <- ets(test, model="MAM")#find test set RMSE
sqrt(mean(retail_hw_test$residuals^2))## [1] 0.07100762#produce forecasts for test data
fc_msd <- retail_hw_train %>%
forecast(new_data = anti_join(timeseries, timeseries_train))#observe whether model beat SNAIVE
retail_hw_train %>% accuracy()## ME RMSE MAE MPE MAPE MASE
## Training set 0.006112595 0.3117354 0.1771523 -0.4216475 6.791803 0.5195141
## ACF1
## Training set -0.0125651fc_msd %>% accuracy(timeseries)## ME RMSE MAE MPE MAPE MASE
## Training set 0.006112595 0.3117354 0.1771523 -0.4216475 6.791803 0.5195141
## Test set -1.687199753 1.8680247 1.6871998 -38.1427206 38.142721 4.9478555
## ACF1 Theil's U
## Training set -0.0125651 NA
## Test set 0.4494349 1.256166Exercise 8.9
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?
#STL decomposition and visualization on the data (no Box-Cox)
myseries_stl <- stl(train, s.window='periodic')
myseries_stl %>% autoplot() #visualization
#ETS with Box Cox
train %>%
ets(lambda = BoxCox.lambda(train)) %>%
forecast(h=10) -> ts_bc_ets
ts_bc_ets %>% autoplot() #visualization
ts_bc_ets %>% accuracy()## ME RMSE MAE MPE MAPE MASE
## Training set 0.04447833 0.307475 0.176615 1.374758 6.765866 0.5179382
## ACF1
## Training set -0.07853798#produce forecasts for test data
fc <- ts_bc_ets %>%
forecast(new_data = anti_join(timeseries, train))
fc %>% accuracy(timeseries)## ME RMSE MAE MPE MAPE MASE
## Training set 0.04447833 0.307475 0.176615 1.374758 6.765866 0.5179382
## Test set -0.80051479 0.898136 0.803767 -18.692926 18.746241 2.3571143
## ACF1 Theil's U
## Training set -0.07853798 NA
## Test set 0.17460590 1.144636After the STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data, the result shows that a Box Cox with ETS produced our best performing model.