Data624 Assignment5

Chapter 7:

suppressMessages(suppressWarnings(library(fpp2)))
suppressMessages(suppressWarnings(library(readxl)))
suppressMessages(suppressWarnings(library(seasonal)))

7.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

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  95741
ses_pigs <- ses(pigs, h = 4)
#model
ses_pigs$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.665

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.

ses_pigs$upper[1, "95%"]
##      95% 
## 119020.8
ses_pigs$lower[1, "95%"]
##      95% 
## 78611.97
#using formula
s <- sd(ses_pigs$residuals)
ses_pigs$mean[1] + 1.96*s
## [1] 118952.8
ses_pigs$mean[1] - 1.96*s
## [1] 78679.97

We can see that there is difference in the value in the above 2 cases.

#plot
autoplot(ses_pigs) +  autolayer(ses_pigs$fitted)

7.5 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       168
#plot
autoplot(books)

The sales of books increased with time and thete is lot of variations with out any particular patern.

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

ses_paperback <- ses(books[, "Paperback"], h = 4)
ses_hardcover <- ses(books[, "Hardcover"], h = 4)
autoplot(books[, "Paperback"], series = "Paperback") +   autolayer(ses_paperback, series = "Paperback") + 
autolayer(books[, "Hardcover"], series = "Hardcover") + autolayer(ses_hardcover, series = "Hardcover", PI = FALSE) +   ylab("Sales") +
ggtitle("paperback and hardcover sales")

We can see the forecast is really flat.

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

sqrt(mean(ses_paperback$residuals^2))
## [1] 33.63769
sqrt(mean(ses_hardcover$residuals^2))
## [1] 31.93101

The residuals of hardcover sales are smaller than the one of paperback sales accprding to the RMSE values

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

holt_paperback <- holt(books[, "Paperback"], h = 4)
holt_hardcover <- holt(books[, "Hardcover"], h = 4)
#plot
autoplot(books[, "Paperback"]) +  autolayer(holt_paperback)

autoplot(books[, "Hardcover"]) +  autolayer(holt_hardcover)

There is linear trend in the forecasts.

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(holt_paperback$residuals^2))
s_hardcover <- sqrt(mean(holt_hardcover$residuals^2))
s_paperback
## [1] 31.13692
s_hardcover
## [1] 27.19358

RMSE values are now lower than the Holt’s method. If there isn’t any particular trend in data, it would be better to use SES method

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

Hardcover sales forcast are better than the paperback sales forcast.

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.

holt_paperback$upper[1, "95%"]
##      95% 
## 275.0205
holt_paperback$lower[1, "95%"]
##     95% 
## 143.913
holt_paperback$mean[1] + 1.96*s_paperback
## [1] 270.4951
holt_paperback$mean[1] - 1.96*s_paperback
## [1] 148.4384
holt_hardcover$upper[1, "95%"]
##      95% 
## 307.4256
holt_hardcover$lower[1, "95%"]
##      95% 
## 192.9222
holt_hardcover$mean[1] + 1.96*s_hardcover
## [1] 303.4733
holt_hardcover$mean[1] - 1.96*s_hardcover
## [1] 196.8745

The prediction interval for the first forecast for each series are almost same for all methods. It is different from the ses case, in which the PI was different when it was calculated by ses function and formula respectively.

7.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.92
autoplot(eggs)

We can see downward trend for the price of eggs here. Using holt function will be good here.

First, just use holt function without using any options.

holt_eggs <- holt(eggs, h = 100)
autoplot(holt_eggs) +   autolayer(holt_eggs$fitted)

he predicted price is going to be below 0.

Use holt function with damped option.

holt_damped_eggs <- holt(eggs, damped = TRUE, h = 100)
autoplot(holt_damped_eggs) +   autolayer(holt_damped_eggs$fitted)

The point forecasts didn’t reflect the existing trend.

Use holt function with Box-Cox transformation.

holt_BoxCox_eggs <- holt(eggs, lambda = BoxCox.lambda(eggs), h = 100)
autoplot(holt_BoxCox_eggs) +   autolayer(holt_BoxCox_eggs$fitted)

Now it reflects the existing trend.

Use holt function with Box-Cox transformation and damped option.

holt_BoxCox_damped_eggs <- holt(eggs,damped = TRUE, lambda = BoxCox.lambda(eggs), h = 100)
autoplot(holt_BoxCox_damped_eggs) +   autolayer(holt_BoxCox_damped_eggs$fitted)

These are decreasing trend and does not look good.

show RMSE values for each model

sqrt(mean(holt_eggs$residuals^2))
## [1] 26.58219
sqrt(mean(holt_damped_eggs$residuals^2))
## [1] 26.54019
sqrt(mean(holt_BoxCox_eggs$residuals^2))
## [1] 1.032217
sqrt(mean(holt_BoxCox_damped_eggs$residuals^2))
## [1] 1.039187

BoxCox transformation improves accuracy of the model. Holt’s method with damped option just prohibits the forecasts to be below 0.

The best model was the Box-Cox transformation with Holt’s linear method.

7.8 Recall your retail time series data (from Exercise 3 in Section 2.10).

retaildata <- readxl::read_excel("C:/Users/rites/Documents/GitHub/Data624_Assignment1/retail.xlsx", skip=1)
## readxl works best with a newer version of the tibble package.
## You currently have tibble v1.4.2.
## Falling back to column name repair from tibble <= v1.4.2.
## Message displays once per session.
myts <- ts(retaildata[,"A3349873A"],frequency=12, start=c(1982,4))
autoplot(myts) + xlab("Time") + ylab("Sales")

a. Why is multiplicative seasonality necessary for this series?

autoplot(myts)

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(myts, seasonal = "multiplicative")
ets_AAdM_retail <- hw(myts,seasonal = "multiplicative",damped = TRUE)
autoplot(ets_AAM_retail)

autoplot(ets_AAdM_retail)

When damped option was used the forcast was slow.

c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?

error_ets_AAM_retail <- tsCV(  myts, hw, h = 1, seasonal = "multiplicative")
error_ets_AAdM_retail <- tsCV( myts, hw, h = 1, seasonal = "multiplicative", damped = TRUE)
sqrt(mean(error_ets_AAM_retail^2, na.rm = TRUE))
## [1] 14.72762
sqrt(mean(error_ets_AAdM_retail^2, na.rm = TRUE))
## [1] 14.94306

RMSE are lmost same. So damped model will be good 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: 24

The residuals don’t look like white noise.

e. 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 3.7?

myts_train <- window(myts, end = c(2010, 12))
myts_test <- window(myts, start = 2011)

Holt-Winters’ method with damped option.

ets_AAdM_retail_train <- hw(myts_train, h = 36, seasonal = "multiplicative", damped = TRUE)
autoplot(ets_AAdM_retail_train)

accuracy(ets_AAdM_retail_train, myts_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(myts_train, h = 36, seasonal = "multiplicative")
autoplot(ets_AAM_retail_train)

accuracy(ets_AAM_retail_train, myts_test)
##                       ME      RMSE       MAE          MPE      MAPE
## Training set  0.03021223  9.107356  6.553533  0.001995484  3.293399
## Test set     78.34068365 94.806617 78.340684 19.945024968 19.945025
##                   MASE       ACF1 Theil's U
## Training set 0.4107058 0.02752875        NA
## Test set     4.9095618 0.52802701  1.613903

This has better accuracy than using the option.

7.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 <- myts_train %>% stlm( s.window = 13, robust = TRUE, method = "ets", lambda = BoxCox.lambda(myts_train) ) %>%
forecast( h = 36, lambda = BoxCox.lambda(myts_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, myts_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.679783
# try forecasting without doing transformation.

fc_stl_ets_retail_train_without_tr <- myts_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,  myts_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.549271

ETS forecasting after STL decomposition ‘without’ Box-Cox transformation yielded better result than when ETS(A, Ad, M) or ETS(A, A, M).