Homework 5 - Exponential Smothing
Name: Charles Ugigabe.
Date: 10/15/23
library(fpp3)
library(tidyverse)Exercises 8.1
Consider the 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 alpha and level 0, and generate forecasts for the next four months.
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.
Solution 8.1
Solution 8.1A
The optimal values are: α = 0.3221247 ℓ = 100646.6 These are the forecasts for the next four months:
victoria.pigs <- aus_livestock %>%
filter(State == "Victoria", Animal == "Pigs")
fit <- victoria.pigs %>%
model(M1 = ETS(Count ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h = 4)fc %>%
autoplot(victoria.pigs)
fit %>% coef()## # A tibble: 2 × 5
## Animal State .model term estimate
## <fct> <fct> <chr> <chr> <dbl>
## 1 Pigs Victoria M1 alpha 0.322
## 2 Pigs Victoria M1 l[0] 100647.Solution 8.1B
The difference in intervals is due to the estimate value of the statistic. This interval is slightly smaller than the one produced by R.
- Calculated Prediction Interval: [76854.452, 113518.663]
- Predicted Prediction Interval:[76854.7889, 113518.32597]
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.07sd <- sqrt(87480760)
upper <- round(fc$.mean[1] + 1.96 * sd, 3)
lower <- round(fc$.mean[1] - 1.96 * sd, 3)
model.ci <- fc %>% hilo(95) %>% pull('95%') %>% head(1)paste0("Calculated Prediction Interval: [", lower, ", ", upper, "]")## [1] "Calculated Prediction Interval: [76854.452, 113518.663]"paste0("Predicted Prediction Interval:", model.ci)## [1] "Predicted Prediction Interval:[76854.7888896402, 113518.325972343]95"Exercise 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.
Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
Compute the RMSE values for the training data.
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.
Compare the forecasts from both methods. Which do you think is best?
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.
Solution 8.5
Solution 8.5A
The US exports have an inconsistent upward trend from 1960 to 2016. There were notable drops in exports within every decades, during the 1980s and late 2000s, reflecting a recession with those period
us.econ <- global_economy %>%
filter(Country == "United States")
us.econ <- us.econ[-58,]
autoplot(us.econ, Exports)
Solution 8.5B
fit <- us.econ %>%
model(PartB = ETS(Exports ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h = 4)
fc %>%
autoplot(us.econ) +
ggtitle("ETS(A, N, N) Forecast")
Solution 8.5C
The RMSE is 0.627 for the model.
accuracy(fit)## # A tibble: 1 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United States PartB Training 0.125 0.627 0.465 1.37 5.10 0.990 0.992 0.239Solution 8.5D
The RMSE for the ETS(A, A, N) is slightly lower than that of the ETS(A, N, N) model which suggest that the model forecast better. Generally, simple exponential smoothing (ETS(A, N, N)) is a better fitting model when the data has no clear trend or seasonal pattern.
fit2 <- us.econ %>%
model(PartD = ETS(Exports ~ error("A") + trend("A") + season("N")))
fc2 <- fit2 %>%
forecast(h = 4)
fc2 %>%
autoplot(us.econ) +
ggtitle("ETS(A, A, N) Forecast")
rbind(data.frame(accuracy(fit)), data.frame(accuracy(fit2)))## Country .model .type ME RMSE MAE MPE
## 1 United States PartB Training 0.12500385 0.6270672 0.4650723 1.37073229
## 2 United States PartD Training 0.01075622 0.6149566 0.4659441 -0.05696615
## MAPE MASE RMSSE ACF1
## 1 5.096552 0.9900993 0.9921329 0.2387874
## 2 5.188477 0.9919554 0.9729717 0.2375941Solution 8.5E
I believe that the ETS(A,A,N) method is better since it captures the increasing trend in the data and has a lower RMSE. Meanwhile, the ETS(A,N,N) method suggests that the exports will stay the same.
fit3 <- us.econ %>%
model(PartB = ETS(Exports ~ error("A") + trend("N") + season("N")),
PartD = ETS(Exports ~ error("A") + trend("A") + season("N")))
fc3 <- fit3 %>%
forecast(h = 4)
fc3 %>%
autoplot(us.econ) +
ggtitle("Part B and Part D Forecast Comparison")
Solution 8.5F
The calculated prediction intervals for both models are very similar to the one produced in R.
- Part B: Calculated Prediction Interval: [10.639, 13.142]
- Part B: Predicted Prediction Interval:[10.6395079134527, 13.1418592559602]
- Part D: Calculated Prediction Interval: [10.757, 13.257]
- Part D: Predicted Prediction Interval:[10.7566652294566, 13.2565616708691]
glance(fit3)## # A tibble: 2 × 10
## Country .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United States PartB 0.408 -88.6 183. 184. 189. 0.393 0.989 0.465
## 2 United States PartD 0.407 -87.5 185. 186. 195. 0.378 0.903 0.466sd <- sqrt(0.4075120)
upper <- round(fc$.mean[1] + 1.96 * sd, 3)
lower <- round(fc$.mean[1] - 1.96 * sd, 3)
model.ci <- fc %>% hilo(95) %>% pull('95%') %>% head(1)sd2 <- sqrt(0.4067128)
upper2 <- round(fc2$.mean[1] + 1.96 * sd2, 3)
lower2 <- round(fc2$.mean[1] - 1.96 * sd2, 3)
model.ci2 <- fc2 %>% hilo(95) %>% pull('95%') %>% head(1)paste0("Part B: Calculated Prediction Interval: [", lower, ", ", upper, "]")## [1] "Part B: Calculated Prediction Interval: [10.639, 13.142]"paste0("Part B: Predicted Prediction Interval:", model.ci)## [1] "Part B: Predicted Prediction Interval:[10.6395079134527, 13.1418592559602]95"paste0("Part D: Calculated Prediction Interval: [", lower2, ", ", upper2, "]")## [1] "Part D: Calculated Prediction Interval: [10.757, 13.257]"paste0("Part D: Predicted Prediction Interval:", model.ci2)## [1] "Part D: Predicted Prediction Interval:[10.7566652294566, 13.2565616708691]95"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.]
- Alpha effects the level of the forecast. Low alpha means more weight is applied to older observations. High alpha means that more weight is given to newer observations
- Beta effects the slope of the forecast. Low beta means more weight is applied to older observations. High beta means that more weight is given to newer observations
- Phi effects the severity of the damping effect. Low phi caues the damping effect to occur sooner. High phi pushes the damping effect further into the future.
- Boxcox seems to provide an exponential forecast that is much higher than the other methods.
Solution 8.6
china_economy <- global_economy %>%
filter(Country == "China") %>%
mutate(GDP_per_capita = GDP/Population)
china_economy %>%
autoplot(GDP_per_capita)
china_economy %>%
model(
`SES` = ETS(GDP_per_capita ~ error("A") + trend("N") + season("N")),
`Damped Holt's method (phi = 0.8)` = ETS(GDP_per_capita ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Damped Holt's method (phi = 0.85)` = ETS(GDP_per_capita ~ error("A") + trend("Ad", phi = 0.85) + season("N")),
`Damped Holt's method (phi = 0.9)` = ETS(GDP_per_capita ~ error("A") + trend("Ad", phi = 0.9) + season("N")),
`Damped Holt's method (phi = 0.98)` = ETS(GDP_per_capita ~ error("A") + trend("Ad", phi = 0.98) + season("N")),
) %>%
forecast(h = 15) %>%
autoplot(china_economy, level = NULL) +
labs(title = "China GDP Per Capita") +
guides(colour = guide_legend(title = "Forecast"))
As we can see here, as we decrease the damping parameter ϕ, the
forecasts begin to level out to a constant at an earlier time. By
forecasting 15 years into the future, we can see that we see significant
tapering at ϕ=0.8. This tapering becomes less evident as we increase
ϕ .
lambda <- china_economy %>%
features(GDP_per_capita, features = guerrero) |>
pull(lambda_guerrero)
china_economy %>%
model(
`SES` = ETS(box_cox(GDP_per_capita, lambda) ~ error("A") + trend("N") + season("N")),
`Damped Holt's method (phi = 0.8)` = ETS(box_cox(GDP_per_capita, lambda) ~ error("A") + trend("Ad", phi = 0.8) + season("N")),
`Damped Holt's method (phi = 0.85)` = ETS(box_cox(GDP_per_capita, lambda) ~ error("A") + trend("Ad", phi = 0.85) + season("N")),
`Damped Holt's method (phi = 0.9)` = ETS(box_cox(GDP_per_capita, lambda) ~ error("A") + trend("Ad", phi = 0.9) + season("N")),
`Damped Holt's method (phi = 0.98)` = ETS(box_cox(GDP_per_capita, lambda) ~ error("A") + trend("Ad", phi = 0.98) + season("N")),
) %>%
forecast(h = 15) %>%
autoplot(china_economy, level = NULL) +
labs(title = "China GDP Per Capita") +
guides(colour = guide_legend(title = "Forecast"))
Applying a Box-Cox Transformation has resulted in the following forecasts. As we increase phi, the slopes of the forecasts increase. With that being said, current UN projections suggest that the population of China will stagnate by 2029. Therefore, the forecasts from Box-Cox transforming the original data should definitely not be reported.
Exercise 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?
Solution 8.7
Multiplicative seasonality is necessary because the seasonality increases as over time. Additive would be preferred if the seasonality variability was fairly constant quarter after quarter. Both the Multiplicative and Damped Multiplicative models performed almost identical based on the plot and accuracy metrics. I would say that the damped effect does not improve the forecasts within the next few years.
aus_production <- aus_production
fit <- aus_production %>%
model(Multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
`Damped Multiplicative` = ETS(Gas ~ error("M") + trend("Ad") + season("M")))
fc <- fit %>%
forecast(h = 16)fc %>%
autoplot(aus_production)
data.frame(accuracy(fit))## .model .type ME RMSE MAE MPE
## 1 Multiplicative Training -0.114886465 4.595113 3.021727 0.1987614
## 2 Damped Multiplicative Training -0.004385562 4.591840 3.031478 0.3258961
## MAPE MASE RMSSE ACF1
## 1 4.082508 0.5420365 0.6059856 -0.01306290
## 2 4.104484 0.5437856 0.6055540 -0.02170345Exercise 8.8
Recall your retail time series data (from Exercise 8 in Section 2.10).
Why is multiplicative seasonality necessary for this series?
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?
Solution 8.8
Solution 8.8A
The variation in the seasonality changes over time, therefore a multiplicative seasonality component is preferred.
set.seed(15)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
autoplot(myseries, Turnover)
Solution 8.8B
The damped method appears to provide consistently lower projections than its non-damped counterpart model.
fit <- myseries %>%
model(`Holt Winters Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Damped Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fc <- fit %>%
forecast(h = 20)
fc %>%
autoplot(myseries)
Solution 8.8C
The Damped Multiplicative model has a lowest RMSE, therefore it fits the training data the best. I would prefer using the Damped Multiplicative model as it appears to be a better fit for the timeseries.
accuracy(fit)## # A tibble: 2 × 12
## State Indus…¹ .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 New So… Furnit… Holt … Trai… 0.759 11.3 7.96 0.0164 4.15 0.524 0.534 0.190
## 2 New So… Furnit… Dampe… Trai… 0.949 11.1 7.83 0.204 4.08 0.515 0.524 0.0261
## # … with abbreviated variable name ¹IndustrySolution 8.8D
The innovation residuals plot shows many spikes in the data, meaning that the residual variance fluctuates. Many lags are outside the blue interval, pointing toward autocorrelation. The histogram appears to be near normal wih a handful of outliers on each side of the mean. The statistical tests run below both reject the null hypothesis (pvalue < 0.05), indicating the possibility of non-zero autocorrelation within the first 24 lags. The residuals appear to not b white noise, indicating that there is most likely a more superior model to use to forecast this data.
fit <- myseries %>%
model(`Damped Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
fit %>%
gg_tsresiduals()
myseries %>% features(Turnover, box_pierce, lag = 24)## # A tibble: 1 × 4
## State Industry bp_stat bp_pv…¹
## <chr> <chr> <dbl> <dbl>
## 1 New South Wales Furniture, floor coverings, houseware and tex… 7963. 0
## # … with abbreviated variable name ¹bp_pvaluemyseries %>% features(Turnover, ljung_box, lag = 24)## # A tibble: 1 × 4
## State Industry lb_stat lb_pv…¹
## <chr> <chr> <dbl> <dbl>
## 1 New South Wales Furniture, floor coverings, houseware and tex… 8219. 0
## # … with abbreviated variable name ¹lb_pvalueSolution 8.8E
Both exponential smoothing approaches produced a better forecast than the naive approach. The Damped Multiplicative model was only a slightly better fit than the naive method, whereas the Holts-Winters Multiplicative model performed the best by far on the test data.
train.data <- myseries %>%
filter(year(Month) <= 2010)
test.data <- myseries %>%
filter(year(Month) > 2010)
fit <- train.data %>%
model(`Holt Winters Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")),
`Damped Multiplicative` = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
`Seasonal Naive` = SNAIVE(Turnover))
fc <- fit %>%
forecast(test.data)
fc %>%
autoplot(myseries, level = NULL)
fc %>% accuracy(test.data)## # A tibble: 3 × 12
## .model State Indus…¹ .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Damped Mu… New … Furnit… Test 57.5 71.4 58.5 14.4 14.7 NaN NaN 0.895
## 2 Holt Wint… New … Furnit… Test 12.4 24.3 19.8 2.78 5.21 NaN NaN 0.686
## 3 Seasonal … New … Furnit… Test 66.9 80.9 67.5 16.9 17.1 NaN NaN 0.902
## # … with abbreviated variable name ¹IndustryExercise 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?
Solution 8.9
we have to first transform the data using Box-Cox transformation. Then we can perform the STL decomposition before extracting the STL components. The STL decomposition provided the best RMSE out of all of the models. The both boxcox transformation models were more effective than the non-transformed model. A combination of model complexity and data manipulation is most likely to produce an optimal forecasting result. It should be noted that box-cox transformation effects the readability of the model results, making it less practical in certain situations.
lambda <- train.data %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
stl.model <- train.data %>%
model(`STL Boxcox` = STL(box_cox(Turnover, lambda)))
fit <- train.data %>%
model(`Simple ETS Boxcox` =ETS(box_cox(Turnover, lambda)),
`Complex ETS Boxcox` =ETS(box_cox(Turnover, lambda) ~ error("M") + trend("A") + season("M")),
`Non Transformed Holt Winters Multiplicative` = ETS(Turnover ~ error("M") + trend("A") + season("M")))stl.model %>%
components() %>%
autoplot()
fc <- fit %>%
forecast(test.data)
fc %>%
autoplot(myseries, level = NULL)
rbind(accuracy(stl.model), accuracy(fit))## # A tibble: 4 × 12
## State Indus…¹ .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 New South … Furnit… STL B… Trai… 0.103 6.88 4.92 -0.148 3.17 0.370 0.374
## 2 New South … Furnit… Simpl… Trai… 1.41 9.71 6.54 0.659 4.16 0.491 0.529
## 3 New South … Furnit… Compl… Trai… 0.253 10.3 7.13 -0.0143 4.56 0.536 0.563
## 4 New South … Furnit… Non T… Trai… -0.0511 10.2 7.02 -0.378 4.52 0.528 0.556
## # … with 1 more variable: ACF1 <dbl>, and abbreviated variable name ¹IndustryError received when trying to get test data forecasting for stl model.
fc <- stl.model %>%
forecast(test.data)## Error in `mutate()`:
## ℹ In argument: `STL Boxcox = (function (object, ...) ...`.
## Caused by error in `UseMethod()`:
## ! no applicable method for 'forecast' applied to an object of class "stl_decomposition"