Data 624 HW 5 MK
Murali Kunissery
3/5/2020
Question 1
- Consider the pigs series - the number of pigs slaughtered in Victoria each month.
a. Use the ses function in R to find the optimal values of alpha and l0, and generate forecasts for the next four months.
str(pigs)## Time-Series [1:188] from 1980 to 1996: 76378 71947 33873 96428 105084 ...head(pigs)## Jan Feb Mar Apr May Jun
## 1980 76378 71947 33873 96428 105084 95741sespigs <- ses(pigs, h = 4)
# check Model
sespigs$model## Simple exponential smoothing
##
## Call:
## ses(y = pigs, h = 4)
##
## Smoothing parameters:
## alpha = 0.2971
##
## Initial states:
## l = 77260.0561
##
## sigma: 10308.58
##
## AIC AICc BIC
## 4462.955 4463.086 4472.665b.Compute a 95% prediction interval for the first forecast using
# 95% prediction interval
sespigs$upper[1, "95%"]## 95%
## 119020.8sespigs$lower[1, "95%"]## 95%
## 78611.97# calculate 95% prediction interval using formula
s <- sd(sespigs$residuals)
sespigs$mean[1] + 1.96*s## [1] 118952.8sespigs$mean[1] - 1.96*s## [1] 78679.97# plot the data, fitted values and forecasts.
autoplot(sespigs) + autolayer(sespigs$fitted)
- 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.
a. Plot the series and discuss the main features of the data.
str(books)## Time-Series [1:30, 1:2] from 1 to 30: 199 172 111 209 161 119 195 195 131 183 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "Paperback" "Hardcover"head(books)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## Paperback Hardcover
## 1 199 139
## 2 172 128
## 3 111 172
## 4 209 139
## 5 161 191
## 6 119 168autoplot(books)
# The sales of both paperback and hardcover have an increasing trend but there is no cyclic or seasonality factor seen.Question 5.
# b. Use the ses function to forecast each series, and plot the forecasts.
sespaperback <- ses(books[, "Paperback"], h = 4)
seshardcover <- ses(books[, "Hardcover"], h = 4)
autoplot(books[, "Paperback"], series = "Paperback") +
autolayer(sespaperback, series = "Paperback") +
autolayer(books[, "Hardcover"], series = "Hardcover") +
autolayer(seshardcover, series = "Hardcover", PI = FALSE) +
ylab("Sales amount") +
ggtitle("Sales of paperback and hardcover books")
Inference: The forecast is flat
c. Compute the RMSE values for the training data in each case.
sqrt(mean(sespaperback$residuals^2))## [1] 33.63769sqrt(mean(seshardcover$residuals^2))## [1] 31.93101Inference: Variance of the residuals of hardcover sales are lesser than paperback.
Question 6.
a. Now apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
holtpaperback <- holt(books[, "Paperback"], h = 4)
holthardcover <- holt(books[, "Hardcover"], h = 4)
autoplot(books[, "Paperback"]) +
autolayer(holtpaperback)
autoplot(books[, "Hardcover"]) +
autolayer(holthardcover)
Inference: The trend is linear.
b. 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.
s_paperback <- sqrt(mean(holtpaperback$residuals^2))
s_hardcover <- sqrt(mean(holthardcover$residuals^2))
s_paperback## [1] 31.13692s_hardcover## [1] 27.19358Inference: For both series, RMSE values became lower when Holt’s method was used. If there is linearly approximable trend in data, it would be better to use Holt’s linear method even if one more parameter is needed than SES. But if there isn’t any particular trend in data, it would be better to use SES method to make the model simpler.
c. Compare the forecasts for the two series using both methods. Which do you think is best? Inference: I think that the forecasts of hardcover sales were better than the ones of paperback sales. Because RMSE value is lower for hardcover sales. And because the forecasts of paperback sales couldn’t reflect the pattern in the data using Holt’s method.
d. 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("95% PI of paperback sales calculated by holt function")## 95% PI of paperback sales calculated by holt functionholtpaperback$upper[1, "95%"]## 95%
## 275.0205holtpaperback$lower[1, "95%"]## 95%
## 143.913writeLines("95% PI of paperback sales calculated by formula")## 95% PI of paperback sales calculated by formulaholtpaperback$mean[1] + 1.96*s_paperback## [1] 270.4951holtpaperback$mean[1] - 1.96*s_paperback## [1] 148.4384writeLines("95% PI of hardcover sales calculated by holt function")## 95% PI of hardcover sales calculated by holt functionholthardcover$upper[1, "95%"]## 95%
## 307.4256holthardcover$lower[1, "95%"]## 95%
## 192.9222writeLines("95% PI of hardcover sales calculated by formula")## 95% PI of hardcover sales calculated by formulaholthardcover$mean[1] + 1.96*s_hardcover## [1] 303.4733holthardcover$mean[1] - 1.96*s_hardcover## [1] 196.8745Inference: In this case, the prediction interval for the first forecast for each series was almost same regardless of calculating method. It is different from the ses case, in which the PI was different when it was calculated by ses function and formula respectively.
Question 7.
- 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.
[Hint: use h=100 when calling holt() so you can clearly see the differences between the various options when plotting the forecasts.]
Which model gives the best RMSE?
str(eggs)## Time-Series [1:94] from 1900 to 1993: 277 315 315 321 315 ...head(eggs)## Time Series:
## Start = 1900
## End = 1905
## Frequency = 1
## [1] 276.79 315.42 314.87 321.25 314.54 317.92autoplot(eggs)
holteggs <- holt(eggs, h = 100)
autoplot(holteggs) +
autolayer(holteggs$fitted)
holtdamped_eggs <- holt(eggs, damped = TRUE, h = 100)
autoplot(holtdamped_eggs) +
autolayer(holtdamped_eggs$fitted)
holtBoxCox_eggs <- holt(eggs,
lambda = BoxCox.lambda(eggs),
h = 100)
autoplot(holtBoxCox_eggs) +
autolayer(holtBoxCox_eggs$fitted)
holtBoxCox_damped_eggs <- holt(
eggs,
damped = TRUE,
lambda = BoxCox.lambda(eggs),
h = 100)
autoplot(holtBoxCox_damped_eggs) +
autolayer(holtBoxCox_damped_eggs$fitted)
# show RMSE values for each model
writeLines("RMSE when using holt function")## RMSE when using holt functionsqrt(mean(holteggs$residuals^2))## [1] 26.58219writeLines("RMSE when using holt function with damped option")## RMSE when using holt function with damped optionsqrt(mean(holtdamped_eggs$residuals^2))## [1] 26.54019writeLines("RMSE when using holt function with Box-Cox transformation")## RMSE when using holt function with Box-Cox transformationsqrt(mean(holtBoxCox_eggs$residuals^2))## [1] 1.032217writeLines("RMSE when using holt function with damped option and Box-Cox transformation")## RMSE when using holt function with damped option and Box-Cox transformationsqrt(mean(holtBoxCox_damped_eggs$residuals^2))## [1] 1.039187Inference: BoxCox transformation captures trend and reflects it to the forecasts. Therefore it improves accuracy of the model. Holt’s method with damped option just prohibits the forecasts to be below 0, not much improving accuracy . The best model was the Box-Cox transformation with Holt’s linear method. It gave plausible point forecasts and prediction intervals. For 100 years’ prediction, Box-Cox transformation did enough damping effect. With damping option together, the point forecast couldn’t follow the existing trend.
Question 8.
- Recall your retail time series data (from Exercise 3 in Section 2.10).
a. Why is multiplicative seasonality necessary for this series?
retail <- readxl::read_excel("retail.xlsx", skip=1)
ts_retail <- ts(retail[, "A3349873A"],frequency = 12,start = c(1982, 4))
autoplot(ts_retail)
Inference: The data show that the seasonality indices increased when the retail sales increased. Multiplicative seasonality can reflect the situation in the model, while additive seasonality can’t.
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
ets_AAM_retail <- hw(ts_retail,
seasonal = "multiplicative")
ets_AAdM_retail <- hw(ts_retail,
seasonal = "multiplicative",
damped = TRUE)
autoplot(ets_AAM_retail)
autoplot(ets_AAdM_retail)
# The forecasts increased more slowly when damped option was used than it wasn't used.c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
error_ets_AAM_retail <- tsCV(
ts_retail,
hw, h = 1, seasonal = "multiplicative"
)
error_ets_AAdM_retail <- tsCV(
ts_retail,
hw, h = 1, seasonal = "multiplicative", damped = TRUE
)
sqrt(mean(error_ets_AAM_retail^2, na.rm = TRUE))## [1] 14.72762sqrt(mean(error_ets_AAdM_retail^2, na.rm = TRUE))## [1] 14.94306Inference: When the RMSE values were compared, they were almost same. Therefore I prefer damped model because it will prohibit the limitless increase of sales forecast.
d. Check that the residuals from the best method look like white noise.
checkresiduals(ets_AAdM_retail)
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 42.932, df = 7, p-value = 3.437e-07
##
## Model df: 17. Total lags used: 24e. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naive approach from Exercise 7 in Section 3.7?
ts_retail_train <- window(ts_retail,
end = c(2010, 12))
ts_retail_test <- window(ts_retail,
start = 2011)
ets_AAdM_retail_train <- hw(ts_retail_train,
h = 36,
seasonal = "multiplicative",
damped = TRUE)
autoplot(ets_AAdM_retail_train)
accuracy(ets_AAdM_retail_train, ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.4556121 8.681456 6.24903 0.2040939 3.151257 0.3916228
## Test set 94.7346169 111.911266 94.73462 24.2839784 24.283978 5.9369594
## ACF1 Theil's U
## Training set -0.01331859 NA
## Test set 0.60960299 1.90013# Holt-Winters' method.
ets_AAM_retail_train <- hw(ts_retail_train,
h = 36,
seasonal = "multiplicative")
autoplot(ets_AAM_retail_train)
accuracy(ets_AAM_retail_train, ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.03021223 9.107356 6.553533 0.001995484 3.293399 0.4107058
## Test set 78.34068365 94.806617 78.340684 19.945024968 19.945025 4.9095618
## ACF1 Theil's U
## Training set 0.02752875 NA
## Test set 0.52802701 1.613903Inference: The accuracy is much better with Holt_winters without the damped setting.
Question 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?
fc_stl_ets_retail_train <- ts_retail_train %>%
stlm(
#made stl model first
s.window = 13,
robust = TRUE,
#designate that the seasonally adjusted data should be forecasted by ETS method.
method = "ets",
lambda = BoxCox.lambda(ts_retail_train)
) %>%
#forecast using stl model
forecast(
h = 36,
lambda = BoxCox.lambda(ts_retail_train)
)## Warning in InvBoxCox(fcast$mean, lambda, biasadj, fcast): biasadj information
## not found, defaulting to FALSE.autoplot(fc_stl_ets_retail_train)
accuracy(fc_stl_ets_retail_train, ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.6782982 8.583559 5.918078 -0.3254076 2.913104 0.3708823
## Test set 82.1015276 98.384220 82.101528 21.0189982 21.018998 5.1452516
## ACF1 Theil's U
## Training set 0.02704667 NA
## Test set 0.52161725 1.679783fc_stl_ets_retail_train_without_tr <-
ts_retail_train %>%
stlm(
s.window = 13,
robust = TRUE,
method = "ets"
) %>%
forecast(h = 36)
autoplot(fc_stl_ets_retail_train_without_tr)
accuracy(fc_stl_ets_retail_train_without_tr,
ts_retail_test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.5020795 10.00826 6.851597 -0.391432 3.489759 0.4293853
## Test set 74.2529959 91.04491 74.252996 18.837766 18.837766 4.6533890
## ACF1 Theil's U
## Training set 0.09741266 NA
## Test set 0.48917501 1.549271Inference: The ETS forcasting after STL decomposition without the Box-cox gives better result
```