DSCI-FXP-SP19- Assignment #5 Exponential Smoothing

Question 1

Data set books contains the daily sales of paperback and hardcover books at the same store. The task is to forecast the next four days’ sales for paperback and hardcover books.

Plot the series and discuss the main features of the data.

p1 <- autoplot(books[,1]) + xlab("Day") + ylab("Sales") + ggtitle("Sales of Paperback Books")

p2 <- autoplot(books[,2]) + xlab("Day") + ylab("Sales") + ggtitle("Sales of Hardcover Books")

p3 <- ggAcf(books[,1], main = NULL)

p4 <- ggAcf(books[,2], main = NULL)

grid.arrange(p1, p2, p3, p4, ncol = 2)

According to the time plots, the sales of paperback and hardcover books both have an increasing trend. However, no seasonal patter is observed in both time series. The correlograms show the same observation.

Use the ses() function to forecast each series, and plot the forecasts.

# Forecast for paperback books
fcast_paper_ses <- ses(books[,1], h = 4)
fcast_paper_ses

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 31       207.1097 162.4882 251.7311 138.8670 275.3523
## 32       207.1097 161.8589 252.3604 137.9046 276.3147
## 33       207.1097 161.2382 252.9811 136.9554 277.2639
## 34       207.1097 160.6259 253.5935 136.0188 278.2005

# Forecast for hardcover books
fcast_hardcover_ses <- ses(books[,2], h = 4)
fcast_hardcover_ses

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 31       239.5601 197.2026 281.9176 174.7799 304.3403
## 32       239.5601 194.9788 284.1414 171.3788 307.7414
## 33       239.5601 192.8607 286.2595 168.1396 310.9806
## 34       239.5601 190.8347 288.2855 165.0410 314.0792

p1 <- autoplot(fcast_paper_ses) + autolayer(fitted(fcast_paper_ses), series = "Fitted") + xlab("Day") + ylab("Paperback Books Sales")

p2 <- autoplot(fcast_hardcover_ses) + autolayer(fitted(fcast_hardcover_ses), series = "Fitted") + xlab("Day") + ylab("Hardcover Books Sales")

grid.arrange(p1, p2, ncol = 2)

Compute the RMSE values for the training data in each case.


writeLines("RMSE for paperback book sales forecast")
## RMSE for paperback book sales forecast
round(accuracy(fcast_paper_ses)[2], 2)
## [1] 33.64

writeLines("RMSE for hardcover book sales forecast")
## RMSE for hardcover book sales forecast
round(accuracy(fcast_hardcover_ses)[2], 2)
## [1] 31.93

The RMSE for the paperback book sales forecast is 33.64, and the RMSE for the hardcover book sales forecast is 31.93.

Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.

# Forecast for paperback books
fcast_paper_holt <- holt(books[,1], h = 4)
fcast_paper_holt

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 31       209.4668 166.6035 252.3301 143.9130 275.0205
## 32       210.7177 167.8544 253.5811 145.1640 276.2715
## 33       211.9687 169.1054 254.8320 146.4149 277.5225
## 34       213.2197 170.3564 256.0830 147.6659 278.7735

# Forecast for hardcover books
fcast_hardcover_holt <- holt(books[,2], h = 4)
fcast_hardcover_holt

##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 31       250.1739 212.7390 287.6087 192.9222 307.4256
## 32       253.4765 216.0416 290.9113 196.2248 310.7282
## 33       256.7791 219.3442 294.2140 199.5274 314.0308
## 34       260.0817 222.6468 297.5166 202.8300 317.3334

p1 <- autoplot(fcast_paper_holt) + autolayer(fitted(fcast_paper_holt), series = "Fitted") + xlab("Day") + ylab("Paperback Books Sales")

p2 <- autoplot(fcast_hardcover_holt) + autolayer(fitted(fcast_hardcover_holt), series = "Fitted") + xlab("Day") + ylab("Hardcover Books Sales")

grid.arrange(p1, p2, ncol = 2)

Compare the RMSE measures of Holt’s method for the two series to those of simple exponential smoothing in the previous question. (Remember that Holt’s method is using one more parameter than SES.) Discuss the merits of the two forecasting methods for these data sets.


writeLines("RMSE for paperback book sales forecast")
## RMSE for paperback book sales forecast
round(accuracy(fcast_paper_holt)[2], 2)
## [1] 31.14

writeLines("RMSE for hardcover book sales forecast")
## RMSE for hardcover book sales forecast
round(accuracy(fcast_hardcover_holt)[2], 2)
## [1] 27.19

The RMSE for the paperback book sales forecast is 31.14, and the RMSE for the hardcover book sales forecast is 27.19. The RMSEs of Holt’s method for the two series are smaller than those of the simple exponential smoothing.

When the time series has a clear trend, it is better to use Holt’s method thatn the simple exponential smoothing. Becasue the SES method cannot capture the trend feature. On the contrary, if no trend is observed over time, the simple exponential smoothing may give a more reasonable forecast.

Compare the forecasts for the two series using both methods. Which do you think is best?

Based on the RMSE, Holt’s method is best for these two series.

Calculate a 95% prediction interval for the first forecast for each series, using the RMSE values and assuming normal errors. Compare your intervals with those produced using ses and holt.


writeLines("Output: 95% Prediction Interval for SES forecast for paperback book sales")
## Output: 95% Prediction Interval for SES forecast for paperback book sales
round(fcast_paper_ses$lower[1, "95%"],2)
##    95% 
## 138.87
round(fcast_paper_ses$upper[1, "95%"],2)
##    95% 
## 275.35

writeLines("Calculated: 95% Prediction Interval for SES forecast for paperback book sales")
## Calculated: 95% Prediction Interval for SES forecast for paperback book sales
round((fcast_paper_ses$mean[1] - 1.96*accuracy(fcast_paper_ses)[2]),2)
## [1] 141.18
round((fcast_paper_ses$mean[1] + 1.96*accuracy(fcast_paper_ses)[2]),2)
## [1] 273.04

writeLines("Output: 95% Prediction Interval for SES forecast for hardcover book sales")
## Output: 95% Prediction Interval for SES forecast for hardcover book sales
round(fcast_hardcover_ses$lower[1, "95%"],2)
##    95% 
## 174.78
round(fcast_hardcover_ses$upper[1, "95%"],2)
##    95% 
## 304.34

writeLines("Calculated: 95% Prediction Interval for SES forecast for hardcover book sales")
## Calculated: 95% Prediction Interval for SES forecast for hardcover book sales
round((fcast_hardcover_ses$mean[1] - 1.96*accuracy(fcast_hardcover_ses)[2]),2)
## [1] 176.98
round((fcast_hardcover_ses$mean[1] + 1.96*accuracy(fcast_hardcover_ses)[2]),2)
## [1] 302.14

writeLines("Output: 95% Prediction Interval for Holt's forecast for paperback book sales")
## Output: 95% Prediction Interval for Holt's forecast for paperback book sales
round(fcast_paper_holt$lower[1, "95%"],2)
##    95% 
## 143.91
round(fcast_paper_holt$upper[1, "95%"],2)
##    95% 
## 275.02

writeLines("Calculated: 95% Prediction Interval for Holt's forecast for paperback book sales")
## Calculated: 95% Prediction Interval for Holt's forecast for paperback book sales
round((fcast_paper_holt$mean[1] - 1.96*accuracy(fcast_paper_holt)[2]),2)
## [1] 148.44
round((fcast_paper_holt$mean[1] + 1.96*accuracy(fcast_paper_holt)[2]),2)
## [1] 270.5

writeLines("Output: 95% Prediction Interval for Holt's forecast for hardcover book sales")
## Output: 95% Prediction Interval for Holt's forecast for hardcover book sales
round(fcast_hardcover_holt$lower[1, "95%"],2)
##    95% 
## 192.92
round(fcast_hardcover_holt$upper[1, "95%"],2)
##    95% 
## 307.43

writeLines("Calculated: 95% Prediction Interval for Holt's forecast for paperback book sales")
## Calculated: 95% Prediction Interval for Holt's forecast for paperback book sales
round((fcast_hardcover_holt$mean[1] - 1.96*accuracy(fcast_hardcover_holt)[2]),2)
## [1] 196.87
round((fcast_hardcover_holt$mean[1] + 1.96*accuracy(fcast_hardcover_holt)[2]),2)
## [1] 303.47

Forecast for the paperback book sales:

Forecast Output 95% PI: SES - [138.87, 275.35], Holt’s - [143.91, 275.02]
Calculated 95% PI: SES - [141.18, 273.04], Holt’s - [148.44, 270.5]

Forecast for the hardcover book sales:

Forecast Output 95% PI: SES - [174.78, 304.34], Holt’s - [192.92, 307.43]
Calculated 95% PI: SES - [176.98, 302.14], Holt’s - [196.87, 303.47]

The calculated prediction intervals are smaller than the forecasting outputs.

Question 2

For this exercise use data set visitors, the monthly Australian short-term overseas visitors data, May 1985-April 2005.

Make a time plot of your data and describe the main features of the series.

autoplot(visitors) + xlab("Month") + ylab("Number of Visitors") + ggtitle("Australian Short-term Overseas Visitors (May 1985 - Apr 2005)")

p1 <- ggsubseriesplot(visitors, main = NULL) + xlab("Month") + ylab("Number of Visitors")

p2 <- ggAcf(window(visitors), main = NULL)

grid.arrange(p1, p2, ncol = 2, top = "Australian Short-term Overseas Visitors (May 1985 - Apr 2005)")

According to the time plot, the time series has an upward trend and strong seasonality. The subseries plot and correlogram confirm the existence of trend and seasonality.The seasonality increases over time.

Split your data into a training set and a test set comprising the last two years of available data. Forecast the test set using Holt-Winters’ multiplicative method.

#Split data into a training set and a test set
visitors_test <- window(visitors, start = c(2004, 1))
visitors_train <- window(visitors, end = c(2003, 12))

autoplot(visitors_train) +
  autolayer(holt(visitors_train, type = "multiplicative"), series = "Holt's method", PI = FALSE) + 
  xlab("Month") + ylab("Number of Visitors") +
  ggtitle("Holt's Forecast for Australian Short-term Overseas Visitors")

Why is multiplicative seasonality necessary here?

Becasue the seasonal variation increases over time.

Forecast the two-year test set using each of the following methods:

(an ETS model; an additive ETS model applied to a Box-Cox transformed series; a seasonal naïve method; an STL decomposition applied to the Box-Cox transformed data followed by an ETS model applied to the seasonally adjusted (transformed) data. )

##ETS Model
fit_ets <- ets(visitors_train)

fc1 <- forecast(fit_ets, h = length(visitors_test))
p1 <- autoplot(fc1, xlab = NULL, ylab = NULL) + ggtitle("ETS Model") + autolayer(visitors_test) -> p1

##ETS Additive model
fit_ets_add <- ets(visitors_train, lambda = BoxCox.lambda(visitors_train), additive.only = TRUE)

fc2 <- forecast(fit_ets_add, h = length(visitors_test))
p2 <- autoplot(fc2, xlab = NULL, ylab = NULL) + ggtitle("ETS Additive Model with Box-Cox Transformation") + autolayer(visitors_test)

##Seasonal Naïve Method
fc3 <- fit_naive_season <- snaive(visitors_train, h = length(visitors_test))

p3 <- autoplot(fc3, xlab = NULL, ylab = NULL) + ggtitle("Seasonal Naïve Model") + autolayer(visitors_test)

##ETS model with STL decomposition and Box-Cox transformation
fc4 <- fit_ets_stl <- stlf(visitors_train, t.window = 13, s.window = "periodic", h = length(visitors_test), robust = TRUE, method = "ets", lambda = BoxCox.lambda(visitors_train))

p4 <- autoplot(fc4, xlab = NULL, ylab = NULL) + ggtitle("ETS Model with STL Decomposition and Box-Cox Transformation") + autolayer(visitors_test)

grid.arrange(p1, p2, p3, p4, ncol = 1)

According to the plots above, the ETS models work better than the Sesonal Naïve Model. It makes sense, because the ETS models have more parameters that can capture more features of the time series.

Which method gives the best forecasts? Does it pass the residual tests?


writeLines("RMSE for the ETS Model:")
## RMSE for the ETS Model:
accuracy(fc1, visitors_test)[, "RMSE"]
## Training set     Test set 
##     14.95912     22.90797

writeLines("RMSE for the ETS Additive model:")
## RMSE for the ETS Additive model:
accuracy(fc2, visitors_test)[, "RMSE"]
## Training set     Test set 
##     15.51547     20.22801

writeLines("RMSE for the Seasonal Naïve Method:")
## RMSE for the Seasonal Naïve Method:
accuracy(fc3, visitors_test)[, "RMSE"]
## Training set     Test set 
##     31.24792     59.02680

writeLines("RMSE for the ETS model with STL decomposition and Box-Cox transformation:")
## RMSE for the ETS model with STL decomposition and Box-Cox transformation:
accuracy(fc4, visitors_test)[, "RMSE"]
## Training set     Test set 
##     15.39434     24.35689

According to the RMSE values for the training data, the ETS model is the best. However, the ETS Additive model has the smallest RMSE for the testing data. Since it is reasonable to evaluate a model using testing data, the ETS Additive model is the best for this time series.

summary(fit_ets_add)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = visitors_train, additive.only = TRUE, lambda = BoxCox.lambda(visitors_train)) 
## 
##   Box-Cox transformation: lambda= 0.2709 
## 
##   Smoothing parameters:
##     alpha = 0.6565 
##     beta  = 1e-04 
##     gamma = 0.1586 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 8.7195 
##     b = 0.0978 
##     s = -0.2696 0.1704 0.2909 -0.1056 1.3525 0.5126
##            0.1142 -0.6103 -0.2138 0.0478 -0.5689 -0.7202
## 
##   sigma:  0.2435
## 
##      AIC     AICc      BIC 
## 597.6209 600.9575 659.0306 
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1.168056 15.51547 11.22989 0.1124176 4.046645 0.4299918
##                     ACF1
## Training set -0.02809955

checkresiduals(fit_ets_add)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 15.475, df = 7, p-value = 0.03037
## 
## Model df: 17.   Total lags used: 24

According to the plots above, the residuals pass the tests as follows:

The residuals are around zero, which means the mean is zero.
The variance is constant over time. No pattern is observed in the correlogram.
The residuals follow a normal distribution.
Compare the same four methods using time series cross-validation with the tsCV()function instead of using a training and test set. Do you come to the same conclusions?

##ETS model
fets <- function(x, h){
  forecast(ets(x), h = h)
}
e1 <- tsCV(visitors, fets, h = 1)

##ETS additive model with Box-Cox transformation
fets_add <- function(x, h){
  forecast(ets(x, lambda = BoxCox.lambda(x), additive.only = TRUE), h = h)
}
e2 <- tsCV(visitors, fets_add, h = 1)

##Seasonal Naïve Method
e3 <- tsCV(visitors, snaive, h = 1)

##ETS model with STL decomposition and Box-Cox transformation
fets_stl <- function(x, h){
  stlf(x, h, t.window = 13, s.window = "periodic", robust = TRUE, method = "ets", lambda = BoxCox.lambda(x))
}
e4 <- tsCV(visitors, fets_stl, h = 1)

#Compare MSE:
writeLines("MSE for the ETS Model:")
## MSE for the ETS Model:
mean(e1^2, na.rm = TRUE)
## [1] 343.3552

writeLines("MSE for the ETS Additive model:")
## MSE for the ETS Additive model:
mean(e2^2, na.rm = TRUE)
## [1] 355.8652

writeLines("MSE for the Seasonal Naïve Method:")
## MSE for the Seasonal Naïve Method:
mean(e3^2, na.rm = TRUE)
## [1] 1074.971

writeLines("MSE for the ETS model with STL decomposition and Box-Cox transformation:")
## MSE for the ETS model with STL decomposition and Box-Cox transformation:
mean(e4^2, na.rm = TRUE)
## [1] 355.965

According to the results of cross-validation, the ETS model is the best. This result is different from the model we select in the previous question. The main reason is that the ETS model has better explanation of the entire time series data. However, when we focus on the testing data, the ETS Additive model works better. Moreover, the ETS Additive model also works well for the training data and in the cross-validation.

In conclusion, if the main objective is to make a good forecasting, the ETS Additive model is the best. However, if we would like to find a model that can best represent the original data, the ETS model is the best choice.

Question 3

Use ets() on the following series:

bicoal, chicken, dole,

Does it always give good forecasts?
Find an example where it does not work well. Can you figure out why?

autoplot(bicoal) + xlab("Year") + ylab("Bituminous Coal Production") + ggtitle("US Bituminous Coal Production (1920-1968)")

fit <- ets(bicoal)
fcast <- forecast(fit, h = 10)
autoplot(fcast) + xlab("Year") + ylab("Bituminous Coal Production") + ggtitle("Forecast for US Bituminous Coal Production")

fcast[['model']]

## ETS(M,N,N) 
## 
## Call:
##  ets(y = bicoal) 
## 
##   Smoothing parameters:
##     alpha = 0.8205 
## 
##   Initial states:
##     l = 542.665 
## 
##   sigma:  0.1262
## 
##      AIC     AICc      BIC 
## 595.2499 595.7832 600.9253

According to the time plot, this time series does not have a clear trend or any seasonality. Thus, the ETS model only has the level parameter, alpha. However, the forecast shown as a straight line fails to capture the variation of the original data. Moreover, the Prediction Interval is unreasonably large. Therefore, the ETS model is not a good choice for this time series.

autoplot(chicken) + xlab("Year") + ylab("Price(dollars)") + ggtitle("Price of Chicken in US (1924 - 1993)")

fcast <- forecast(ets(chicken), h = 10)

autoplot(fcast) + xlab("Year") + ylab("Price(dollars)") + ggtitle("Forecast for Price of Chicken in US")

fcast[['model']]

## ETS(M,N,N) 
## 
## Call:
##  ets(y = chicken) 
## 
##   Smoothing parameters:
##     alpha = 0.98 
## 
##   Initial states:
##     l = 159.8322 
## 
##   sigma:  0.1691
## 
##      AIC     AICc      BIC 
## 635.2382 635.6018 641.9836

According to the time plot, the time series has a clear decreasing trend after 1945, without seasonality. The ETS model assigns a large value to alpha, which means the forecast heavily depends on the most current observations. In general, this makes sense, because the original data has high variance over time, and the forecast fits the overal trend. However, the prediction interval is too wide for this forecast. The lower bound is even smaller than zero. Therefore, the ETS model is not a good choice for this data.

autoplot(dole) + xlab("Month") + ylab("Number of People") + ggtitle("Total of People on Unemployment Benefits in Australia (Jan 1965 - Jul 1992)")

fit <- ets(dole)
autoplot(fit)

fcast <- forecast(fit) 
autoplot(fcast) + xlab("Month") + ylab("Number of People") + ggtitle("Forecast for Total of People on Unemployment Benefits in Australia")

fcast[['model']]

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = dole) 
## 
##   Smoothing parameters:
##     alpha = 0.697 
##     beta  = 0.124 
##     gamma = 0.303 
##     phi   = 0.902 
## 
##   Initial states:
##     l = 2708.6621 
##     b = 836.017 
##     s = 1.0404 0.8893 0.9103 1.0301 1.0576 1.0584
##            0.9801 0.9632 1.021 0.9838 1.0145 1.0514
## 
##   sigma:  0.0935
## 
##      AIC     AICc      BIC 
## 10602.67 10604.30 10676.19

According to the time plot, the original data has an upward trend and strong seasonality. The forecast of the ETS model looks fine. However, there might be a concern of underestimation, because there is a dramatical increase since 1990. But the forecast fails to capture this trend. Moreover, the prediction interval is extremely large. This is attributed to the large variation of the seasonal component.

fit2 <- ets(dole, lambda = BoxCox.lambda(dole))
autoplot(fit2)

fcast2 <- forecast(fit2)
autoplot(fcast2) + xlab("Month") + ylab("Number of People") + ggtitle("Forecast for Total of People on Unemployment Benefits in Australia (Box-Cox Transformation")

fcast2[['model']]

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = dole, lambda = BoxCox.lambda(dole)) 
## 
##   Box-Cox transformation: lambda= 0.3291 
## 
##   Smoothing parameters:
##     alpha = 0.9988 
##     beta  = 0.1275 
##     gamma = 1e-04 
##     phi   = 0.8052 
## 
##   Initial states:
##     l = 48.412 
##     b = 1.0772 
##     s = 3.2062 -4.3768 -4.9978 -3.7341 -1.6098 0.147
##            0.7791 -0.2256 0.3576 0.4401 4.0227 5.9915
## 
##   sigma:  2.8435
## 
##      AIC     AICc      BIC 
## 3607.299 3608.927 3680.820

To better address the varying seasonality, Box-Cox transformation is performed before fitting the ETS model. This model gives a better forecast than the previous. The seasonal component is consistant over time, and the prediction interval is more reasonable.

In conclusion, the ETS model will work better when the time series has a clear trend and strong seasonality. In addition, Box-Cox transformation can help handle large variation of seasonality.

Question 4

Compare ets(), snaive() and stlf() on the following three time series. For stlf(), you might need to use a Box-Cox transformation. Use a test set of three years to decide what gives the best forecasts. bricksq, dole, a10.

p1 <- autoplot(bricksq) + xlab("Month") + ylab("Revenue Enplanements (millions)") + ggtitle("US Domestic Monthly Revenue Enplanements (1996 - 2000)")

##ets()
fit_ets <- ets(bricksq, lambda = BoxCox.lambda(bricksq))
fc1 <- forecast(fit_ets)

p2 <- autoplot(fc1, xlab = NULL, ylab = NULL) + ggtitle("ets() Model")

##snaive()
fc2 <- fit_naive_season <- snaive(bricksq)

p3 <- autoplot(fc2, xlab = NULL, ylab = NULL) + ggtitle("snaive() Model")

##stlf()
fc3 <- fit_ets_stl <- stlf(bricksq, t.window = 13, s.window = "periodic", robust = TRUE, lambda = BoxCox.lambda(bricksq))

p4 <- autoplot(fc3, xlab = NULL, ylab = NULL) + ggtitle("stlf() Model")

grid.arrange(p1, p2, p3, p4, ncol = 1)

writeLines("Evaluation of ets() model")

## Evaluation of ets() model

accuracy(fc1)

##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1.752939 21.37562 15.21701 0.4282784 3.762617 0.4274122
##                 ACF1
## Training set 0.12363

writeLines("Evaluation of snaive() model")

## Evaluation of snaive() model

accuracy(fc2)

##                    ME    RMSE      MAE      MPE     MAPE MASE    ACF1
## Training set 6.834437 48.7093 35.60265 1.513132 8.678982    1 0.80655

writeLines("Evaluation of stlf() model")

## Evaluation of stlf() model

accuracy(fc3)

##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1.783691 21.68564 15.20726 0.3737139 3.716862 0.4271384
##                    ACF1
## Training set 0.09877582

According to the results, the ets() model with Box-Cox transformation is the best. The stlf() model has very similar forecast.


p1 <- autoplot(dole) + xlab("Month") + ylab("Revenue Enplanements (millions)") + ggtitle("US Domestic Monthly Revenue Enplanements (1996 - 2000)")

##ets()
fit_ets <- ets(dole, lambda = BoxCox.lambda(dole))
fc1 <- forecast(fit_ets)

p2 <- autoplot(fc1, xlab = NULL, ylab = NULL) + ggtitle("ets() Model")

##snaive()
fc2 <- fit_naive_season <- snaive(dole)

p3 <- autoplot(fc2, xlab = NULL, ylab = NULL) + ggtitle("snaive() Model")

##stlf()
fc3 <- fit_ets_stl <- stlf(dole, t.window = 13, s.window = "periodic", robust = TRUE, lambda = BoxCox.lambda(dole))

p4 <- autoplot(fc3, xlab = NULL, ylab = NULL) + ggtitle("stlf() Model")

grid.arrange(p1, p2, p3, p4, ncol = 1)


writeLines("Evaluation of ets() model")
## Evaluation of ets() model
accuracy(fc1)
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1214.927 11120.48 6607.818 0.4749541 5.321139 0.1602151
##                   ACF1
## Training set 0.1952945

writeLines("Evaluation of snaive() model")
## Evaluation of snaive() model
accuracy(fc2)
##                   ME     RMSE      MAE      MPE     MAPE MASE      ACF1
## Training set 21826.9 73188.35 41243.41 4.481956 27.38484    1 0.9833331

writeLines("Evaluation of stlf() model")
## Evaluation of stlf() model
accuracy(fc3)
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 774.3351 9024.676 5651.671 0.1919982 5.690328 0.1370321
##                   ACF1
## Training set 0.1217972

According to the results, the stlf() model with Box-Cox transformation is the best.However, the prediction interval is fairly wide.


p1 <- autoplot(a10) + xlab("Month") + ylab("Revenue Enplanements (millions)") + ggtitle("US Domestic Monthly Revenue Enplanements (1996 - 2000)")

##ets()
fit_ets <- ets(a10, lambda = BoxCox.lambda(dole))
fc1 <- forecast(fit_ets)

p2 <- autoplot(fc1, xlab = NULL, ylab = NULL) + ggtitle("ets() Model")

##snaive()
fc2 <- fit_naive_season <- snaive(a10)

p3 <- autoplot(fc2, xlab = NULL, ylab = NULL) + ggtitle("snaive() Model")

##stlf()
fc3 <- fit_ets_stl <- stlf(a10, t.window = 13, s.window = "periodic", robust = TRUE, lambda = BoxCox.lambda(a10))

p4 <- autoplot(fc3, xlab = NULL, ylab = NULL) + ggtitle("stlf() Model")

grid.arrange(p1, p2, p3, p4, ncol = 1)


writeLines("Evaluation of ets() model")
## Evaluation of ets() model
accuracy(fc1)
##                      ME      RMSE       MAE        MPE    MAPE      MASE
## Training set 0.03287204 0.8193706 0.5282667 -0.3610049 4.69851 0.4078429
##                    ACF1
## Training set -0.1078884

writeLines("Evaluation of snaive() model")
## Evaluation of snaive() model
accuracy(fc2)
##                    ME     RMSE     MAE      MPE     MAPE MASE      ACF1
## Training set 1.229776 1.729655 1.29527 10.94158 11.44434    1 0.3727544

writeLines("Evaluation of stlf() model")
## Evaluation of stlf() model
accuracy(fc3)
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 0.02521363 0.7884153 0.5052229 -0.182305 4.374387 0.3900522
##                    ACF1
## Training set -0.1605709

According to the results, the stlf() model with Box-Cox transformation is the best.

Question 5

For this exercise use data set ukcars, the quarterly UK passenger vehicle production data from 1977Q1-2005Q1.

Plot the data and describe the main features of the series.

autoplot(ukcars) + xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("Quarterly UK Passenger Car Production (1st QT 1977 - 1st QT 2005)")

p1 <- ggsubseriesplot(ukcars, main = NULL) + xlab("Quarter") + ylab("Passenger Car Production (thousands)")

p2 <- ggAcf(ukcars, main = NULL)

grid.arrange(p1, p2, ncol = 2)

According to the time plot, the passenger car production has an increasing trend between 1980 and 2000. In 2000, there is a dramatical decrease. This time series has strong seasonality that the third quater usually has the lowest production, as shown in the subseries plot and correlogram.

Decompose the series using STL and obtain the seasonally adjusted data.

fit_stl <- stl(ukcars, s.window = 13, robust = TRUE)
autoplot(fit_stl) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("STL Decomposition of UK Passenger Car Production")

ssadj <- seasadj(fit_stl)

Forecast the next two years of the series using an additive damped trend method applied to the seasonally adjusted data. (This can be done in one step using stlf()with arguments etsmodel=“AAN”, damped=TRUE.)

fit_aan <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=TRUE)

fcast <- forecast(fit_aan, h = 8)

autoplot(fcast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("AAN Forecast for UK Passenger Car Production")

accuracy(fcast)

##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 2.404023 24.67648 19.6208 0.2220182 6.391966 0.6394325
##                    ACF1
## Training set 0.03554102

Forecast the next two years of the series using Holt’s linear method applied to the seasonally adjusted data (as before but with damped=FALSE).

fit_hlm <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=FALSE)

fcast <- forecast(fit_hlm, h = 8)

autoplot(fcast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("Holt's Linear Method for UK Passenger Car Production")

accuracy(fcast)

##                      ME     RMSE     MAE        MPE     MAPE      MASE
## Training set -0.1069739 24.79848 19.2994 -0.6297295 6.362326 0.6289582
##                    ACF1
## Training set 0.02121772

Now use ets() to choose a seasonal model for the data.

fit_ets <- ets(ukcars)

fcast <- forecast(fit_ets, h = 8)

fcast[["model"]]

## ETS(A,N,A) 
## 
## Call:
##  ets(y = ukcars) 
## 
##   Smoothing parameters:
##     alpha = 0.6199 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 314.2568 
##     s = -1.7579 -44.9601 21.1956 25.5223
## 
##   sigma:  25.9302
## 
##      AIC     AICc      BIC 
## 1277.752 1278.819 1296.844

autoplot(fcast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("ETS Method for UK Passenger Car Production")

accuracy(fcast)

##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1.313884 25.23244 20.17907 -0.1570979 6.629003 0.6576259
##                    ACF1
## Training set 0.02573334

Compare the RMSE of the ETS model with the RMSE of the models you obtained using STL decompositions. Which gives the better in-sample fits?

According to the evaluation results, the additive damped trend method gives the best fits.

Compare the forecasts from the three approaches? Which seems most reasonable?

In general, these three forecasts are very similar. Based on the RMSE, the additive damped trend method is the best.

Check the residuals of your preferred model.

fit_aan <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=TRUE)

fcast <- forecast(fit_aan, h = 8)

checkresiduals(fcast)

## Warning in checkresiduals(fcast): The fitted degrees of freedom is based on
## the model used for the seasonally adjusted data.

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(A,Ad,N)
## Q* = 19.758, df = 3, p-value = 0.0001905
## 
## Model df: 5.   Total lags used: 8

According to the residual plots, the residuals vary around zero. There is no trend or pattern observed over time and in the correlogram. The distribution of residuals is slightly left skewed. But overall, the residuals meet the requirements.

Question 6

For this exercise use data set eggs, the price of a dozen eggs in the United States from 1900-1993. Experiment with the various options in the holt() 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 argument is doing to the forecasts. Which model gives the best RMSE?


fcast1 <- holt(eggs)
fcast1[["model"]]
## Holt's method 
## 
## Call:
##  holt(y = eggs) 
## 
##   Smoothing parameters:
##     alpha = 0.8124 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 314.7232 
##     b = -2.7222 
## 
##   sigma:  27.1665
## 
##      AIC     AICc      BIC 
## 1053.755 1054.437 1066.472
p1 <- autoplot(fcast1) + autolayer(fitted(fcast1), series = "Fitted") + xlab("Year") + ylab("Dollars") + ggtitle("Forecast using Holt's Method")

fcast2 <- holt(eggs, damped = TRUE)
fcast2[["model"]]
## Damped Holt's method 
## 
## Call:
##  holt(y = eggs, damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.8462 
##     beta  = 0.004 
##     phi   = 0.8 
## 
##   Initial states:
##     l = 276.9842 
##     b = 4.9966 
## 
##   sigma:  27.2755
## 
##      AIC     AICc      BIC 
## 1055.458 1056.423 1070.718
p2 <- autoplot(fcast2) + autolayer(fitted(fcast2), series = "Fitted") + xlab("Year") + ylab("Dollars") + ggtitle("Forecast using Holt's Method with Damped")

fcast3 <- holt(eggs, lambda = BoxCox.lambda(eggs))
p3 <- autoplot(fcast3) + autolayer(fitted(fcast3), series = "Fitted") + xlab("Year") + ylab("Dollars") + ggtitle("Forecast using Holt's Method with Box-Cox Transformation")

fcast4 <- holt(eggs, lambda = BoxCox.lambda(eggs), damped = TRUE)
p4 <- autoplot(fcast4) + autolayer(fitted(fcast4), series = "Fitted") + xlab("Year") + ylab("Dollars") + ggtitle("Forecast using Holt's Method with Damped and Box-Cox Transformation")

grid.arrange(p1, p2, p3, p4, ncol = 2)

According to the forecasts, the Holt’s method with Box-Cox transformation gives the best fit with smaller prediction interval.

When comparing figure 1 vs. figure 2 and figure 3 vs figure 4, the “damped” parameter impacts the slope of the trend. That is, the forecasting trend is not linear over time after damped, which is more realistic.

Box-Cox transformation can help reduce the prediction interval significantly, when the data has large variance.


writeLines("RMSE for Holt's Method")
## RMSE for Holt's Method
accuracy(fcast1)[2]
## [1] 26.58219

writeLines("RMSE for Holt's Method with Damped")
## RMSE for Holt's Method with Damped
accuracy(fcast2)[2]
## [1] 26.54019

writeLines("RMSE for Holt's Method with Box-Cox Transformation")
## RMSE for Holt's Method with Box-Cox Transformation
accuracy(fcast3)[2]
## [1] 26.39376

writeLines("RMSE for Holt's Method with Damped and Box-Cox Transformation")
## RMSE for Holt's Method with Damped and Box-Cox Transformation
accuracy(fcast4)[2]
## [1] 26.53321

The Holt’s method with Box-Cox transformation gives the best forecast.