Homework 5
Gurpreet Singh
March 7, 2019
7.1 Consider the pigs series - the number of pigs slaughtered in Victoria each month.
Pigs is the time series data for monthly number of pigs slaughtered in Victoria, Australia (Jan 1980 -Aug 1995)
class(pigs)## [1] "ts"summary(pigs)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 33873 79080 91662 90640 101493 120184a.Use the ses() function in R to find the optimal values of \(\alpha\) and \(\l_0\), and generate forecasts for the next four months.
fc <- ses(pigs, h=4)
## smoothing parameter alpha
fc$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## forecast for next four months
forecast(fc)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 1995 98816.41 85605.43 112027.4 78611.97 119020.8
## Oct 1995 98816.41 85034.52 112598.3 77738.83 119894.0
## Nov 1995 98816.41 84486.34 113146.5 76900.46 120732.4
## Dec 1995 98816.41 83958.37 113674.4 76092.99 121539.8b. Compute a 95% prediction interval for the first forecast using y ± 1.96 s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
s <- sd(fc$residuals)
mean_fc <- fc$mean[1]
int_fc <- c(mean_fc - (1.96*s),mean_fc + (1.96*s))
int_fc## [1] 78679.97 118952.84The 95% prediction interval for the first forecast using y ± 1.96 s are almost similar as generated by R.
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.
class(books)## [1] "mts" "ts"summary(books)## Paperback Hardcover
## Min. :111.0 Min. :128.0
## 1st Qu.:167.2 1st Qu.:170.5
## Median :189.0 Median :200.5
## Mean :186.4 Mean :198.8
## 3rd Qu.:207.2 3rd Qu.:222.0
## Max. :247.0 Max. :283.0a. Plot the series and discuss the main features of the data.
autoplot(books)
Both the series shows fluctuations with increasing trends. We do not see any seasonality in the data. No indications of cyclic pattern either. Hardcover has lower sales in the begining and tends to increase as time increases as compared to paperback.
b. Use the ses() function to forecast each series, and plot the forecasts.
fc_pb <- ses(books[,1], h=4)
fc_hc <- ses(books[,2], h=4)
forecast(fc_pb)## 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.2005forecast(fc_hc)## 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.0792pb1 <- autoplot(fc_pb) +
autolayer(fitted(fc_pb), series = "Pred. Paperback") +
xlab("Time") + ylab("Sales")
hc1 <- autoplot(fc_hc) +
autolayer(fitted(fc_hc), series = "Pred. Hardcover") +
xlab("Time") + ylab("Sales")
grid.arrange(pb1,hc1)
There are fluctuations in the hardcover predictions as compared to smoother forecasts in the paperback. In addition, the hardcover is following the tremnd more closer than the paperback, providing the inference that it might represent the actual sales.
c. Compute the RMSE values for the training data in each case.
rmse_pb <- accuracy(fc_pb)[2]
rmse_pb## [1] 33.63769rmse_hc <- accuracy(fc_hc)[2]
rmse_hc## [1] 31.93101The rmse for hardcover was better than the paperback which supports the inferences we obtained from the plots.
6. a.Now apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
holt_hc <- holt(books[,1], h=4)
holt_pb <- holt(books[,2], h=4)
holt_hc## 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.7735holt_pb## 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.3334b.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.
rmse_h_hc <- accuracy(holt_hc)[2]
rmse_h_hc## [1] 31.13692rmse_h_pb <- accuracy(holt_pb)[2]
rmse_h_pb## [1] 27.19358The RMSE values for holt method are better in both series. In case of hardcover, holt method outperforms the ses method.
c.Compare the forecasts for the two series using both methods. Which do you think is best?
pb2 <- autoplot(holt_pb) +
autolayer(fitted(fc_pb), series = "Pred. Paperback") +
xlab("Time") + ylab("Sales")
grid.arrange(pb1,pb2)
hc2 <- autoplot(holt_hc) +
autolayer(fitted(fc_hc), series = "Pred. Hardcover") +
xlab("Time") + ylab("Sales")
grid.arrange(hc1,hc2)
Comparing both models, in terms of RMSE, holt is the best model. In terms of fitted plot, holt outperforms the hardcover, however if we look at the paperback in case of holt method the forecast does not seems to represent the actual series, which is a little concerning in an attempt to choose the best model.
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.
### from formula
s_fc_pb <- sd(fc_pb$residuals)
mean_fc_pb <- fc_pb$mean[1]
int_fc_pb <- c(mean_fc_pb - (1.96*s_fc_pb),mean_fc_pb + (1.96*s_fc_pb))
int_fc_pb## [1] 141.5964 272.6230#from model
forecast(fc_pb)$lower[1, "95%"]## 95%
## 138.867forecast(fc_pb)$upper[1, "95%"]## 95%
## 275.3523s_fc_hc <- sd(fc_hc$residuals)
mean_fc_hc <- fc_hc$mean[1]
int_fc_hc <- c(mean_fc_hc - (1.96*s_fc_hc),mean_fc_hc + (1.96*s_fc_hc))
int_fc_hc## [1] 178.5848 300.5354#from model
forecast(fc_hc)$lower[1, "95%"]## 95%
## 174.7799forecast(fc_hc)$upper[1, "95%"]## 95%
## 304.3403s_holt_pb <- sd(holt_pb$residuals)
mean_holt_pb <- holt_pb$mean[1]
int_holt_pb <- c(mean_holt_pb - (1.96*s_holt_pb),mean_holt_pb + (1.96*s_holt_pb))
int_holt_pb## [1] 195.9640 304.3838#from model
forecast(holt_pb)$lower[1, "95%"]## 95%
## 192.9222forecast(holt_pb)$upper[1, "95%"]## 95%
## 307.4256s_holt_hc <- sd(holt_hc$residuals)
mean_holt_hc <- holt_hc$mean[1]
int_holt_hc <- c(mean_holt_hc - (1.96*s_holt_hc),mean_holt_hc + (1.96*s_holt_hc))
int_holt_hc## [1] 147.8390 271.0945#from model
forecast(holt_hc)$lower[1, "95%"]## 95%
## 143.913forecast(holt_hc)$upper[1, "95%"]## 95%
## 275.0205Lower and Upper intervals for both the time series using ses and holt methods are comparitively similar.
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?
holt_r<-holt(eggs,h=100)
holt_damped<-holt(eggs,damped = TRUE,h=100)
holt_bc <-holt(eggs, lambda=BoxCox.lambda(eggs),h=100)
holt_bc_damped <-holt(eggs, lambda=BoxCox.lambda(eggs),damped=TRUE,h=100)
a1 <- autoplot(holt_r)
a2 <- autoplot(holt_damped)
a3 <- autoplot(holt_bc)
a4<- autoplot(holt_bc_damped)
grid.arrange(a1, a2, a3,a4, nrow = 2)
accuracy(holt_r)## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202accuracy(holt_damped)## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358accuracy(holt_bc)## ME RMSE MAE MPE MAPE MASE
## Training set 0.7736844 26.39376 18.96387 -1.072416 9.620095 0.9354593
## ACF1
## Training set 0.03887152accuracy(holt_bc_damped)## ME RMSE MAE MPE MAPE MASE
## Training set -0.8200445 26.53321 19.45654 -2.019718 9.976131 0.9597618
## ACF1
## Training set 0.005852382Four models were used for the eggs time series data. Holt Method - Holt method was applied, the model was predicting negative forecasts.
Holt with Damped - The damped method combined with holt was resolving the issue of negative forecast, the trend however does not seem to follow the actual pattern as the trend line seems to be straighten a lot.
Holt with BOxCox - Holt method was combined with the BoxCox method and it seems to provide better model in terms of forecast and trend.
Holt with BoxCox and Damped - Holt was combined with boxcox and damped, the model was providing better inference than model one and model two. The results seems to be comparable with model 3 (boxcox alone)
Comparing the accuracy of all four models reveals that the RMSE was alsmost similar, however holt with box-cox have slight better value which supports our exploration form autoplots. Model 3 with boxcox will be our best model.
8.Recall your retail time series data (from Exercise 3 in Section 2.10).
retaildata <- readxl::read_excel("C:/Users/Gurpreet/Documents/Data624/retail.xlsx", skip=1)
myts <- ts(retaildata[,"A3349335T"],
frequency=12, start=c(1982,4))autoplot(myts)
a.Why is multiplicative seasonality necessary for this series?
It is clear from the graph that seasonality variations are changing with increase in time. In that case, multiplicative seasonality is the best approach because seasonal variations are not constant and additive method can handle constant seasonal variations only.
b.Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
hw_myts <- hw(myts, seasonal = "multiplicative", h=100)
hw_myts_damped <- hw(myts, damped =TRUE, seasonal = "multiplicative", h=100)
autoplot(myts) + autolayer(hw_myts, series='Retail Data Multiplicative', PI=FALSE) +
autolayer(hw_myts_damped, series='Retail Data Damped', PI=FALSE) 
We were not able to make a clear comparison from the autoplot. In order to compare two methods closely, the autoplots for next 100 forecasts were plotted. In damped method, forecasts were steadily increasing as compared to the multiplicative method.
c.Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
os_hw <- forecast(hw_myts, h=1)
accuracy(os_hw)## ME RMSE MAE MPE MAPE MASE
## Training set 0.9212824 25.20381 18.77683 0.06856226 1.979315 0.3016982
## ACF1
## Training set -0.1217931os_hw_damped <- forecast(hw_myts_damped, h=1)
accuracy(os_hw_damped)## ME RMSE MAE MPE MAPE MASE
## Training set 2.9059 25.10059 18.74334 0.2366077 1.993423 0.3011601
## ACF1
## Training set -0.1394101The RMSE for both the models is almost similar. In terms of RMSE, it might not be a better approach to decide the best fit model. If we look up the autoplot in earlier step, we can assume that the damped model is trying to prevent over-forecast done by the multiplicative model. In addition, for longer forecast horizons, damped method is the best. If in our case in part b and c (this part), we choose the forecast value to be 100 and we were clearly able to get an inference that the damped method is dampening the trend. In this case for h = 1, we can not make a decision for the best model. If we look up the autoplot, damped is the best model.
d.Check that the residuals from the best method look like white noise.
checkresiduals(hw_myts_damped)
##
## Ljung-Box test
##
## data: Residuals from Damped Holt-Winters' multiplicative method
## Q* = 285.62, df = 7, p-value < 2.2e-16
##
## Model df: 17. Total lags used: 24#whitenoise.test(hw_myts_damped)It is evident that the residuals from our damped method do not look like white noise. ACF plot reveals we the the number of spikes lying outside the bounds seems more than 5%.
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 8 in Section 3.7?
myts_train <- window(myts, end = c(2010, 12))
myts_test <- window(myts, start = 2011)myts_train_hw <- hw(myts_train, damped = TRUE, seasonal = "multiplicative")
accuracy(myts_train_hw, myts_test)## ME RMSE MAE MPE MAPE MASE
## Training set 4.089724 25.68057 18.92510 0.371369 2.057716 0.3068053
## Test set -27.787782 41.08034 32.12435 -1.279976 1.487404 0.5207858
## ACF1 Theil's U
## Training set -0.05009924 NA
## Test set 0.13519074 0.3334222myts_train_sn <- snaive(myts_train)
accuracy(myts_train_sn,myts_test)## ME RMSE MAE MPE MAPE MASE
## Training set 61.56787 72.20702 61.68438 6.388722 6.404105 1.000000
## Test set 97.44583 109.62545 100.02917 4.629852 4.751209 1.621629
## ACF1 Theil's U
## Training set 0.6018274 NA
## Test set 0.2686595 0.9036205We can not beat the snaive model with our current damped model. Comparison of test set RMSE reveals that snaive is far more better than the damped model.
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?
lambda1 <- BoxCox.lambda(myts_train)
fc_bc<-BoxCox(myts_train,lambda1)
stl_bc <- stlm(fc_bc)
forecast(stl_bc)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2011 17.30171 17.19937 17.40405 17.14519 17.45822
## Feb 2011 16.94176 16.83622 17.04730 16.78036 17.10316
## Mar 2011 17.30992 17.20128 17.41856 17.14377 17.47607
## Apr 2011 17.14459 17.03293 17.25625 16.97382 17.31536
## May 2011 17.21550 17.10090 17.33011 17.04023 17.39078
## Jun 2011 17.05559 16.93811 17.17307 16.87592 17.23526
## Jul 2011 17.23865 17.11837 17.35893 17.05469 17.42261
## Aug 2011 17.29598 17.17296 17.41901 17.10783 17.48414
## Sep 2011 17.22419 17.09847 17.34990 17.03193 17.41645
## Oct 2011 17.45532 17.32697 17.58366 17.25903 17.65160
## Nov 2011 17.50675 17.37582 17.63767 17.30651 17.70698
## Dec 2011 17.92089 17.78743 18.05435 17.71678 18.12500
## Jan 2012 17.53678 17.40083 17.67272 17.32886 17.74469
## Feb 2012 17.17683 17.03844 17.31522 16.96518 17.38848
## Mar 2012 17.54499 17.40419 17.68578 17.32966 17.76032
## Apr 2012 17.37966 17.23649 17.52282 17.16071 17.59861
## May 2012 17.45057 17.30508 17.59606 17.22806 17.67308
## Jun 2012 17.29066 17.14287 17.43845 17.06464 17.51668
## Jul 2012 17.47372 17.32367 17.62377 17.24424 17.70320
## Aug 2012 17.53105 17.37877 17.68333 17.29816 17.76394
## Sep 2012 17.45926 17.30478 17.61373 17.22301 17.69551
## Oct 2012 17.69039 17.53374 17.84703 17.45082 17.92995
## Nov 2012 17.74182 17.58303 17.90060 17.49897 17.98466
## Dec 2012 18.15596 17.99506 18.31686 17.90988 18.40204summary(stl_bc)## Length Class Mode
## stl 1380 mstl numeric
## model 19 ets list
## modelfunction 1 -none- function
## lambda 0 -none- NULL
## x 345 ts numeric
## series 1 -none- character
## m 1 -none- numeric
## fitted 345 ts numeric
## residuals 345 ts numericstl_bc%>%forecast()%>%autoplot()
## ets after extracting seasonal component
fc_dec<-decompose(fc_bc)
fc_seas<-fc_dec$x - fc_dec$seasonal
fc_ets<-ets(fc_seas,model="ZZZ")
forecast_ets <-forecast(fc_ets,h=50, lambda=lambda1)## Warning in forecast.ets(fc_ets, h = 50, lambda = lambda1): biasadj
## information not found, defaulting to FALSE.autoplot(fc_ets)
The plot with stl to boxcox transformation appears to be forecasting and represent the actual series. Due to some error in stl function, we extrcated the seasonal component manually and applied the ets method to the boxcox transformed series, after extracting seasonal adjusted data, the rmse of this method was better than all the methods. The ets method was the best method outperforming our previous best models.