Week 5 HW

HA 7.5 , 7.6 and 7.10

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.

    1. Plot the series and discuss the main features of the data.
library(fpp2)

autoplot(books)

We are working with a bivariate time series that contains the number of book sales for paperback and hardcover books. Both book format sales are partioned with a daily time interval over the course of 30 days.

Our basic plot reveals that both hardcover and paperback books have an overall positive trend, however we see large fluctuations happening daily. We are unable to tell if there exists some sort of seasonal influence because we are limited to 30 days.

decompose <- ts(books, frequency=7)

autoplot(decompose(decompose[,1]),main = "Paperback")

autoplot(decompose(decompose[,2]), main = "Hardcover")

The decomposition confirms our assumption of a mostly increasing trend for hardcover sales but fluctuation with paperback sales. We could also seeevidence of seasonality every two or three days. As mentioned before, we can’t see any longer term seasonality since we are restricted to 30 days.

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

We will consider the default alpha value and plot a forecase for 4 periods. We will leave out using optimal or toying with different values of alpha since it is not a problem requirement.

paperback <- ses(books[, 1], h = 4)

hardcover <- ses(books[, 2], h = 4)

autoplot(books[, 1], series = "Paperback") +   autolayer(paperback, series = "Paperback")+ ylab("Sales") +
ggtitle("Paperback Sales");

autoplot(books[, 2], series = "Hardcover") +   autolayer(hardcover, series = "Hardcover")+ ylab("Sales") +
ggtitle("Hardcover Sales");

summary(paperback);
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = books[, 1], h = 4) 
## 
##   Smoothing parameters:
##     alpha = 0.1685 
## 
##   Initial states:
##     l = 170.8271 
## 
##   sigma:  34.8183
## 
##      AIC     AICc      BIC 
## 318.9747 319.8978 323.1783 
## 
## Error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303
##                    ACF1
## Training set -0.2117522
## 
## Forecasts:
##    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
summary(hardcover)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = books[, 2], h = 4) 
## 
##   Smoothing parameters:
##     alpha = 0.3283 
## 
##   Initial states:
##     l = 149.2861 
## 
##   sigma:  33.0517
## 
##      AIC     AICc      BIC 
## 315.8506 316.7737 320.0542 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887
##                    ACF1
## Training set -0.1417763
## 
## Forecasts:
##    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

Both forecasts are underperforming. We can see that the trend is both non-increasing and non-decreasing.

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

Recall: https://www.statisticshowto.datasciencecentral.com/rmse/

To summarize, the RMSE can be computed with the following formula:

\[ RMSE=\sqrt(\frac{\sum_{i=1}^{N}{(Z_{fi}-Z_{sdi})^2}}{N}) \]

print(paste0("Paperback RMSE: ",sqrt(mean(paperback$residuals^2))));
## [1] "Paperback RMSE: 33.637686782912"
print(paste0("Hardcover RMSE: ",sqrt(mean(hardcover$residuals^2))))
## [1] "Hardcover RMSE: 31.9310149844547"

The RMSE is the standard deviation of prediction errors. We generally consider a smaller RmSE to be associated with a more accurate prediction or in our case forecast. IF we look at the numbers we obtained, the RMSE is less with the hardcover sales forecasts. If we examine our plots from part b, paperback book sales forecast tends to be flat.

7.6

    1. Now apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
hmodel_paperback <- holt(books[, 1], h = 4)

hmodel_hardcover <- holt(books[, 2], h = 4)

autoplot(books[, 1], series = "Paperback") +   autolayer(hmodel_paperback, series = "Paperback")+ ylab("Sales") +
ggtitle("Paperback Sales via Holt's Method");

autoplot(books[, 2], series = "Hardcover") +   autolayer(hmodel_hardcover, series = "Hardcover")+ ylab("Sales") +
ggtitle("Hardcover Sales via Holt's Method");

summary(hmodel_paperback);
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = books[, 1], h = 4) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 170.699 
##     b = 1.2621 
## 
##   sigma:  33.4464
## 
##      AIC     AICc      BIC 
## 318.3396 320.8396 325.3456 
## 
## Error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -3.717178 31.13692 26.18083 -5.508526 15.58354 0.6602122
##                    ACF1
## Training set -0.1750792
## 
## Forecasts:
##    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
summary(hmodel_hardcover)
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = books[, 2], h = 4) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 147.7935 
##     b = 3.303 
## 
##   sigma:  29.2106
## 
##      AIC     AICc      BIC 
## 310.2148 312.7148 317.2208 
## 
## Error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set -0.1357882 27.19358 23.15557 -2.114792 12.1626 0.6908555
##                     ACF1
## Training set -0.03245186
## 
## Forecasts:
##    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

As opposed to the SES method, we can see a positive linear trend in the forecast for both paperback book sales and hardcover book sales. The positive trend forecast is even greater for the Hardcover paperback sales.

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

Lets Recall the RMSE we generated using SES then compare to RMSE via Holt’s.

print(paste0("Paperback RMSE from SES: ",sqrt(mean(paperback$residuals^2))));
## [1] "Paperback RMSE from SES: 33.637686782912"
print(paste0("Hardcover RMSE from SES: ",sqrt(mean(hardcover$residuals^2))))
## [1] "Hardcover RMSE from SES: 31.9310149844547"
print(paste0("Paperback RMSE from Holt's Method: ",sqrt(mean(hmodel_paperback$residuals^2))));
## [1] "Paperback RMSE from Holt's Method: 31.1369230162347"
print(paste0("Hardcover RMSE from Holt's Method: ",sqrt(mean(hmodel_hardcover$residuals^2))))
## [1] "Hardcover RMSE from Holt's Method: 27.1935779818511"

The RSME generated from the forecast via Holt’s method overall is lower than when generated via SES. Hardcover sales forecast has the biggest decrease in RMSE.

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

Using both methods, hardcover sales forecast are better than paperback sales forecast. When we factor in the RMSE metric, Holt’s method produces a better forecast. By definition of RMSE, lower valuesindicate a better fit.

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

Lets look at the intervals found using r.

ses_paperback_one <- ses(books[, 1], h = 1)
s1<-sd((ses(books[,1], h=1))$residuals)

ses_hardcover_one <- ses(books[, 2], h = 1)
s2<-sd((ses(books[,2], h=1))$residuals)

holt_paperback_one<-holt(books[, 1], h = 1)
h1<-sd((ses(books[,1], h=1))$residuals)

holt_hardcover_one<-holt(books[, 2], h = 1)
h2<-sd((ses(books[,2], h=1))$residuals)
print("SES Paperback: 1 Forecast ")
## [1] "SES Paperback: 1 Forecast "
summary(ses_paperback_one)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = books[, 1], h = 1) 
## 
##   Smoothing parameters:
##     alpha = 0.1685 
## 
##   Initial states:
##     l = 170.8271 
## 
##   sigma:  34.8183
## 
##      AIC     AICc      BIC 
## 318.9747 319.8978 323.1783 
## 
## Error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303
##                    ACF1
## Training set -0.2117522
## 
## Forecasts:
##    Point Forecast    Lo 80    Hi 80   Lo 95    Hi 95
## 31       207.1097 162.4882 251.7311 138.867 275.3523

lower is 138.867 and higher is 275.3523. How does this compare to our own calculation?

print(paste0("lower Confidence Interval: ", ses(books[,1],h=1)$mean[1]-1.96*s1))
## [1] "lower Confidence Interval: 141.596371851814"
print(paste0("upper Confidence Interval: ", ses(books[,1],h=1)$mean[1]+1.96*s1))
## [1] "upper Confidence Interval: 272.622958138623"
print("SES Hardcover: 1 Forecast ")
## [1] "SES Hardcover: 1 Forecast "
summary(ses_hardcover_one)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = books[, 2], h = 1) 
## 
##   Smoothing parameters:
##     alpha = 0.3283 
## 
##   Initial states:
##     l = 149.2861 
## 
##   sigma:  33.0517
## 
##      AIC     AICc      BIC 
## 315.8506 316.7737 320.0542 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887
##                    ACF1
## Training set -0.1417763
## 
## Forecasts:
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 31       239.5601 197.2026 281.9176 174.7799 304.3403

The lower is 174.7799 and the upper is 304.3403. How does it compare to r our calculation?

print(paste0("lower Confidence Interval: ", ses(books[,2],h=1)$mean[1]-1.96*s2))
## [1] "lower Confidence Interval: 178.584828931457"
print(paste0("upper Confidence Interval: ", ses(books[,2],h=1)$mean[1]+1.96*s2))
## [1] "upper Confidence Interval: 300.535355072066"
print("Holt Paperback: 1 Forecast ")
## [1] "Holt Paperback: 1 Forecast "
summary(holt_paperback_one)
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = books[, 1], h = 1) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 170.699 
##     b = 1.2621 
## 
##   sigma:  33.4464
## 
##      AIC     AICc      BIC 
## 318.3396 320.8396 325.3456 
## 
## Error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -3.717178 31.13692 26.18083 -5.508526 15.58354 0.6602122
##                    ACF1
## Training set -0.1750792
## 
## Forecasts:
##    Point Forecast    Lo 80    Hi 80   Lo 95    Hi 95
## 31       209.4668 166.6035 252.3301 143.913 275.0205

lower is 143.913 and upper is 275.0205. How does it compare to our own?

print(paste0("lower Confidence Interval: ", holt(books[,1],h=1)$mean[1]-1.96*h1))
## [1] "lower Confidence Interval: 143.953472377646"
print(paste0("upper Confidence Interval: ", holt(books[,1],h=1)$mean[1]+1.96*h1))
## [1] "upper Confidence Interval: 274.980058664455"
print("Holt Hardcover: 1 Forecast ")
## [1] "Holt Hardcover: 1 Forecast "
summary(holt_hardcover_one)
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = books[, 2], h = 1) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 147.7935 
##     b = 3.303 
## 
##   sigma:  29.2106
## 
##      AIC     AICc      BIC 
## 310.2148 312.7148 317.2208 
## 
## Error measures:
##                      ME     RMSE      MAE       MPE    MAPE      MASE
## Training set -0.1357882 27.19358 23.15557 -2.114792 12.1626 0.6908555
##                     ACF1
## Training set -0.03245186
## 
## Forecasts:
##    Point Forecast   Lo 80    Hi 80    Lo 95    Hi 95
## 31       250.1739 212.739 287.6087 192.9222 307.4256

low is 192.9222 and high is 307.4256. How does it compare?

print(paste0("lower Confidence Interval: ", holt(books[,2],h=1)$mean[1]-1.96*h2))
## [1] "lower Confidence Interval: 189.198609142052"
print(paste0("upper Confidence Interval: ", holt(books[,2],h=1)$mean[1]+1.96*h2))
## [1] "upper Confidence Interval: 311.14913528266"

7.10

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

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

We are working with Quarterly UK passenger car production from 1977 to 2005. Our data is time partioned by year.

autoplot(ukcars)

We can observe evidence of seasonal effects within the series. We see a decreasing trend up until the 1980. 1980 is an inflection point since the trend begins to increase up until the sudden drop in 2000. 2000 is another inflection point since we observe production to begin increasing again.

    1. Decompose the series using STL and obtain the seasonally adjusted data.
decompose <- stl(ukcars, s.window="periodic", robust=TRUE)

seasonal <- decompose$time.series[,1]

ukcars_adjusted <- ukcars - seasonal

autoplot(ukcars_adjusted)

    1. 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.)
mod1 <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=TRUE)

mod1_forecast <- forecast(mod1, h = 8)

autoplot(mod1_forecast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("STLF Forecast for UK Passenger Car Production")+   autolayer(mod1_forecast, series = "Car Production");

summary(mod1)
## 
## Forecast method: STL +  ETS(A,Ad,N)
## 
## Model Information:
## ETS(A,Ad,N) 
## 
## Call:
##  ets(y = x, model = etsmodel, damped = TRUE, allow.multiplicative.trend = allow.multiplicative.trend) 
## 
##   Smoothing parameters:
##     alpha = 0.608 
##     beta  = 1e-04 
##     phi   = 0.9127 
## 
##   Initial states:
##     l = 335.0123 
##     b = -8.9945 
## 
##   sigma:  25.2412
## 
##      AIC     AICc      BIC 
## 1270.717 1271.509 1287.081 
## 
## Error measures:
##                    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
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005 Q2       423.1304 390.7825 455.4783 373.6585 472.6023
## 2005 Q3       364.9947 327.1359 402.8535 307.0947 422.8947
## 2005 Q4       407.1140 364.4492 449.7789 341.8638 472.3643
## 2006 Q1       430.6043 383.6215 477.5872 358.7502 502.4584
## 2006 Q2       423.1287 372.1917 474.0657 345.2272 501.0301
## 2006 Q3       364.9931 310.3868 419.5994 281.4800 448.5062
## 2006 Q4       407.1126 349.0678 465.1573 318.3408 495.8843
## 2007 Q1       430.6030 369.3119 491.8941 336.8663 524.3397
    1. 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).
mod2 <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=FALSE)

mod2_forecast <- forecast(mod2, h = 8)

autoplot(mod2_forecast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("Holt's Linear Method for UK Passenger Car Production") + autolayer(mod2_forecast, series = "Car Production");

summary(mod2)
## 
## Forecast method: STL +  ETS(A,A,N)
## 
## Model Information:
## ETS(A,A,N) 
## 
## Call:
##  ets(y = x, model = etsmodel, damped = FALSE, allow.multiplicative.trend = allow.multiplicative.trend) 
## 
##   Smoothing parameters:
##     alpha = 0.6458 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 333.1107 
##     b = 0.7271 
## 
##   sigma:  25.2494
## 
##      AIC     AICc      BIC 
## 1269.832 1270.392 1283.469 
## 
## Error measures:
##                      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
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005 Q2       424.1567 391.7983 456.5151 374.6688 473.6446
## 2005 Q3       366.7474 328.2264 405.2684 307.8346 425.6602
## 2005 Q4       409.5931 365.7661 453.4201 342.5654 476.6207
## 2006 Q1       433.8097 385.2517 482.3677 359.5466 508.0728
## 2006 Q2       427.0603 374.1917 479.9289 346.2047 507.9158
## 2006 Q3       369.6509 312.7965 426.5054 282.6995 456.6024
## 2006 Q4       412.4966 351.9168 473.0764 319.8478 505.1455
## 2007 Q1       436.7133 372.6233 500.8033 338.6961 534.7305
    1. Now use ets() to choose a seasonal model for the data.
mod3 <- ets(ukcars)

mod3_forecast <- forecast(mod3, h = 8)

autoplot(mod3_forecast) +  xlab("Quarter") + ylab("Passenger Car Production (thousands)") + ggtitle("ETS Method for UK Passenger Car Production") + autolayer(mod3_forecast, series = "Car Production");

summary(mod3_forecast)
## 
## Forecast method: ETS(A,N,A)
## 
## Model Information:
## 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 
## 
## Error measures:
##                    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
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005 Q2       427.4885 394.2576 460.7195 376.6662 478.3109
## 2005 Q3       361.3329 322.2353 400.4305 301.5383 421.1275
## 2005 Q4       404.5358 360.3437 448.7280 336.9497 472.1219
## 2006 Q1       431.8154 383.0568 480.5741 357.2455 506.3854
## 2006 Q2       427.4885 374.5571 480.4200 346.5369 508.4401
## 2006 Q3       361.3329 304.5345 418.1313 274.4672 448.1986
## 2006 Q4       404.5358 344.1174 464.9542 312.1338 496.9378
## 2007 Q1       431.8154 367.9809 495.6500 334.1890 529.4419
    1. 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?
print("Additive")
## [1] "Additive"
accuracy(mod1_forecast);
##                    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
print("Holt")
## [1] "Holt"
accuracy(mod2_forecast);
##                      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
print("ETS")
## [1] "ETS"
accuracy(mod3_forecast)
##                    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

Based on the RMSE scores from the three models, the additive model has the lowest RMSE. By definition, the lower the RMSE, the better the forecast. It should be noted that overall, there is no major difference between RMSE scores. It is more accurate to say that the additive model is marginally better than the Holts linear model or ETs model.

    1. Compare the forecasts from the three approaches? Which seems most reasonable?
plot(mod1, main="Additive Model");

plot(mod2, main="Holt Linear Model");

plot(mod3_forecast, main="ETS Model")

Out of the three models we have fit, the additive model is still the better of the three, even if it is marginally better than the Holt linear method and ETS model.

    1. Check the residuals of your preferred model.
selected_fit <- stlf(ukcars, s.window = 13, robust = TRUE, etsmodel="AAN", damped=TRUE)

selected_forecast <- forecast(selected_fit, h = 8)

checkresiduals(selected_forecast)

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

Aside from the residuals having a minor left skew, we can tell that the underlying assumption of residuals is met overall. There is certainly the potnetial for improvement, however we leave that as a suggestion. The rest of the output reveals that there is no clear and cut trend.