624_5
Albina Gallyavova
3/8/2020
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 α and ℓ0, and generate forecasts for the next four months.
autoplot(pigs)
fc <- ses(pigs, h=4)
#fc$model
summary(fc)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## 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
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 385.8721 10253.6 7961.383 -0.922652 9.274016 0.7966249 0.01282239
##
## Forecasts:
## 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.8THe optimal values for alpha and level is 0.2971 and 77260.0561 respectively.
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.
autoplot(fc) +
autolayer(fitted(fc), series="Fitted")
# using model results
fc$upper[1, "95%"]## 95%
## 119020.8fc$lower[1, "95%"]## 95%
## 78611.97# using ^y±1.96s
fc$mean[1] + 1.96*sd(fc$residuals)## [1] 118952.8fc$mean[1] - 1.96*sd(fc$residuals)## [1] 78679.97Interval produced by R
- Upper: 1.190208410^{5}
- Lower: 7.86119710^{4}
Interval produced using ^y±1.96s
- Upper: 1.189528410^{5}
- Lower: 7.86799710^{4}
Interval produced by R is wider than when using 1.96sd.
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.
autoplot(books)
There is a visible upward trend in the data for both categories, however, there is no apparent seasonality or cyclicality.
b. Use the ses() function to forecast each series, and plot the forecasts.
paper_fc_ses <- ses(books[,'Paperback'], h=4)
hard_fc_ses <- ses(books[,'Hardcover'], h=4)
autoplot(books[, "Paperback"], series = "Paperback") +
autolayer(paper_fc_ses, series = "Paperback") +
ggtitle('Paperback sales and forecast')
autoplot(books[, "Hardcover"], series = "Hardcover") +
autolayer(hard_fc_ses, series = "Hardcover") +
ggtitle('Hardcover sales and forecast')
c. Compute the RMSE values for the training data in each case.
RMSE for paperback training: 33.64
RMSE for hardcover training: 31.93
round(forecast::accuracy(paper_fc_ses),2)[2]## [1] 33.64round(forecast::accuracy(hard_fc_ses),2)[2]## [1] 31.936
We will continue with the daily sales of paperback and hardcover books in data set books.
#### a. Apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
paper_fc_holt <- holt(books[,'Paperback'], h=4)
hard_fc_holt <- holt(books[,'Hardcover'], h=4)
summary(paper_fc_holt)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = books[, "Paperback"], 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.7735summary(hard_fc_holt)##
## Forecast method: Holt's method
##
## Model Information:
## Holt's method
##
## Call:
## holt(y = books[, "Hardcover"], 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.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.
Paper forecast produced with ses of 33.64 is greater than that producted by holt of 31.14. The same can be said about hardcover sales forecast.
round(forecast::accuracy(paper_fc_ses),2)[2]## [1] 33.64round(forecast::accuracy(hard_fc_ses),2)[2]## [1] 31.93round(forecast::accuracy(paper_fc_holt),2)[2]## [1] 31.14round(forecast::accuracy(hard_fc_holt),2)[2]## [1] 27.19c. Compare the forecasts for the two series using both methods. Which do you think is best?
a <- autoplot(books[, "Paperback"], series = "Paperback") +
autolayer(paper_fc_ses, series = "Paperback_ses", PI=FALSE) +
autolayer(paper_fc_holt, series = "Paperback_holt", PI=FALSE) +
ggtitle('Paperback sales and forecast')
b <- autoplot(books[, "Hardcover"], series = "Hardcover") +
autolayer(hard_fc_ses, series = "Hardcover_ses", PI=FALSE) +
autolayer(hard_fc_holt, series = "Hardcover_holt", PI=FALSE) +
ggtitle('Hardcover sales and forecast')
grid.arrange(a, b, ncol = 1)
Just visually, holt method seems to work better than ses as it includes the upward trend in the forecast.
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.
paper_rmse_up <- paper_fc_holt$mean[1] + round(forecast::accuracy(paper_fc_holt)[2],2)*1.96
paper_rmse_low <- paper_fc_holt$mean[1] - round(forecast::accuracy(paper_fc_holt)[2],2)*1.96
hard_rmse_up <- hard_fc_holt$mean[1] + round(forecast::accuracy(hard_fc_holt)[2],2)*1.96
hard_rmse_low <- hard_fc_holt$mean[1] - round(forecast::accuracy(hard_fc_holt)[2],2)*1.96
# holt PI
paper_holt_up <- paper_fc_holt$upper[1, "95%"]
paper_holt_low <- paper_fc_holt$lower[1, "95%"]
hard_holt_up <- hard_fc_holt$upper[1, "95%"]
hard_holt_low <- hard_fc_holt$lower[1, "95%"]
# ses PI
paper_ses_up <- paper_fc_ses$upper[1, "95%"]
paper_ses_low <- paper_fc_ses$lower[1, "95%"]
hard_ses_up <- hard_fc_ses$upper[1, "95%"]
hard_ses_low <- hard_fc_ses$lower[1, "95%"]
df <- data.frame("Method" = c('Paper using RMSE','Paper Holt','Paper ses',
'Hardcover using RMSE','Hardcover Holt','Hardcover ses'),
"Upper" = c(paper_rmse_up,paper_holt_up,paper_ses_up,
hard_rmse_up,hard_holt_up,hard_ses_up),
"Lower" = c(paper_rmse_low,paper_holt_low,paper_ses_low,
hard_rmse_low,hard_holt_low,hard_ses_low))
datatable(df) %>%
formatRound('Upper',3) %>%
formatRound('Lower',3)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.]
autoplot(eggs)
eggs_holt <- holt(eggs,h=100)
eggs_holt_damp <- holt(eggs,damped=TRUE, h=100)
eggs_holt_bc <- holt(eggs,lambda=BoxCox.lambda(eggs),h=100)
eggs_holt_damp_bc <- holt(eggs,damped=TRUE, lambda=BoxCox.lambda(eggs),h=100)
autoplot(eggs) +
autolayer(eggs_holt, series='Holt') +
autolayer(eggs_holt_damp, series='Damped') +
autolayer(eggs_holt_bc, series='BoxCox transformed')+
autolayer(eggs_holt_damp_bc, series='BoxCox transformed and damped')+
theme(legend.position="bottom") 
Which model gives the best RMSE?
res <- data.frame(rbind(forecast::accuracy(eggs_holt)[2],
forecast::accuracy(eggs_holt_damp)[2],
forecast::accuracy(eggs_holt_bc)[2],
forecast::accuracy(eggs_holt_damp_bc)[2]))
colnames(res) <- "RMSE"
rownames(res) <- c('Holt','Holt damped','Holt with Box-Cox','Holt damped with Box-Cox')
datatable(res) %>% formatRound('RMSE',3)Holt with Box-Cox produces lowest RMSE. Although it is the lowest, looking at the prediction intervals, Box-Cox transformed combined with damping logically make more sense as damping flattens out a downward trend preventing prices from going negative which is more realistic.
8
Recall your retail time series data (from Exercise 3 in Section 2.10).
download.file('https://otexts.com/fpp2/extrafiles/retail.xlsx','retail.xlsx')
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
myts <- ts(retaildata[,"A3349398A"],
frequency=12, start=c(1982,4))
autoplot(myts)
a. Why is multiplicative seasonality necessary for this series?
Time series exhibits seasonal variation that increases proportionally to the time series making multiplicative method more appropriate.
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
hw <- hw(myts, seasonal='multiplicative')
hw_damped <- hw(myts, damped=TRUE, seasonal='multiplicative')
autoplot(myts) +
autolayer(hw,series='HW',PI=FALSE) +
autolayer(hw_damped,series='HW_damped',PI=FALSE)
Damped method seems to slow down the trend a little which could be beneficial for longer term forecasts.
c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
holt_rmse <- forecast::accuracy(hw(myts, seasonal='multiplicative',h=1))[2]
holt_damped_rmse <- forecast::accuracy(hw(myts, damped=TRUE, seasonal='multiplicative',h=1))[2]Holt-Winters’ works better given lower RMSE.
Holt-Winters’ RMSE: 29.4305131 Holt-Winters’ damped RMSE: 29.6308699
d. Check that the residuals from the best method look like white noise.
checkresiduals(hw)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 244.86, df = 8, p-value < 2.2e-16
##
## Model df: 16. Total lags used: 24Although generally residuals appear to be noise and normal, as evident from distribution and lack of runs in lags plot, the variation is slightly wider in the beginning of the series.
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)
hw_fc <- hw(myts.train,h=36)
hw_test_rmse <- rmse(myts.test,hw_fc$mean)
hw_damped_fc <- hw(myts.train,damped=TRUE,h=36)
hw_d_test_rmse <- rmse(myts.test,hw_damped_fc$mean)
naive_fc <- snaive(myts.train,h=36)
naive_rmse <- rmse(myts.test,naive_fc$mean)
autoplot(myts) +
autolayer(naive_fc,series='Naive') +
autolayer(hw_fc,series='HW') +
autolayer(hw_damped_fc,series='HW damped') +
theme(legend.title=element_blank())
checkresiduals(fc)
##
## Ljung-Box test
##
## data: Residuals from Simple exponential smoothing
## Q* = 55.356, df = 22, p-value = 0.0001057
##
## Model df: 2. Total lags used: 24hw_test_rmse## [1] 49.39973hw_d_test_rmse ## [1] 100.5575naive_rmse## [1] 180.1991Naive test RMSE: 180.1990574
Both Holt-Winters’ and Holt-Winters’ damped perform significantly better than Naive due to its ability to capture trend. As evident from residuals plots, Naive’s errors increase due to projections staying relatively flat whereas time series keep increasing.
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?
stlm takes a time series y, applies an STL decomposition, and models the seasonally adjusted data using the model passed as modelfunction or specified using method. It returns an object that includes the original STL decomposition and a time series model fitted to the seasonally adjusted data. This object can be passed to the forecast.stlm for forecasting.
fit <- stlm(myts.train,method='ets',lambda='auto',allow.multiplicative.trend = TRUE) #,lambda='auto',biasadj=TRUE
fc <- forecast(fit, h = 36)
summary(fc)##
## Forecast method: STL + ETS(A,A,N)
##
## Model Information:
## ETS(A,A,N)
##
## Call:
## ets(y = x, model = etsmodel, allow.multiplicative.trend = allow.multiplicative.trend)
##
## Smoothing parameters:
## alpha = 0.3216
## beta = 1e-04
##
## Initial states:
## l = 8.7058
## b = 0.0118
##
## sigma: 0.0477
##
## AIC AICc BIC
## -78.17328 -77.99629 -58.95556
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1.077737 26.34835 19.95571 -0.02085 1.675349 0.2655976
## ACF1
## Training set -0.04737883
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2011 2607.340 2543.931 2672.142 2510.920 2707.019
## Feb 2011 2388.081 2326.492 2451.108 2294.461 2485.064
## Mar 2011 2607.981 2538.377 2679.267 2502.200 2717.697
## Apr 2011 2514.255 2444.074 2586.212 2407.629 2625.036
## May 2011 2544.683 2471.034 2620.266 2432.815 2661.076
## Jun 2011 2450.722 2376.956 2526.506 2338.707 2567.456
## Jul 2011 2557.418 2478.274 2638.790 2437.260 2682.785
## Aug 2011 2597.863 2515.117 2683.008 2472.264 2729.072
## Sep 2011 2566.704 2482.444 2653.483 2438.838 2700.463
## Oct 2011 2712.660 2621.725 2806.372 2574.687 2857.130
## Nov 2011 2758.187 2663.489 2855.851 2614.532 2908.779
## Dec 2011 3148.159 3039.043 3260.726 2982.646 3321.743
## Jan 2012 2759.345 2659.963 2861.996 2608.646 2917.691
## Feb 2012 2528.761 2434.700 2626.022 2386.171 2678.836
## Mar 2012 2760.019 2656.179 2867.433 2602.620 2925.777
## Apr 2012 2661.467 2558.840 2767.716 2505.943 2825.465
## May 2012 2693.464 2587.686 2803.047 2533.192 2862.638
## Jun 2012 2594.650 2490.367 2702.776 2436.677 2761.612
## Jul 2012 2706.855 2596.595 2821.236 2539.850 2883.498
## Aug 2012 2749.381 2635.606 2867.478 2577.079 2931.792
## Sep 2012 2716.619 2602.121 2835.549 2543.255 2900.352
## Oct 2012 2870.063 2747.863 2997.045 2685.056 3066.255
## Nov 2012 2917.918 2791.929 3048.908 2727.202 3120.334
## Dec 2012 3327.637 3183.889 3477.095 3110.039 3558.591
## Jan 2013 2919.134 2789.222 3054.368 2722.542 3128.176
## Feb 2013 2676.721 2554.667 2803.905 2492.069 2873.371
## Mar 2013 2919.843 2786.139 3059.189 2717.577 3135.308
## Apr 2013 2816.249 2684.977 2953.165 2617.701 3028.000
## May 2013 2849.885 2715.446 2990.177 2646.575 3066.887
## Jun 2013 2746.003 2614.221 2883.625 2546.751 2958.918
## Jul 2013 2863.962 2725.434 3008.681 2654.528 3087.878
## Aug 2013 2908.661 2766.476 3057.271 2693.725 3138.626
## Sep 2013 2874.225 2731.833 3023.142 2659.010 3104.703
## Oct 2013 3035.491 2884.270 3193.681 2806.948 3280.337
## Nov 2013 3085.774 2930.564 3248.209 2851.229 3337.220
## Dec 2013 3516.104 3339.876 3700.505 3249.787 3801.541plot(forecast(fit, h = 36))
lambda <- fc$lambda
pred <- fc$mean
#pred_transf <- InvBoxCox(pred, lambda=lambda, biasadj = FALSE, fvar = NULL)
actual <- myts.test
forecast::accuracy(fc$mean ,myts.test)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -135.2801 153.7568 135.6001 -5.045165 5.057412 0.6869809 0.7791571ets_rmse <- rmse(actual, pred)ETS achieves better test RMSE of 153.7567905 than Naive of 180.1990574, however, regular Holt-Winters still outperforms all other methods for these data.
df <- data.frame(rbind(hw_test_rmse,hw_d_test_rmse ,naive_rmse,ets_rmse))
colnames(df) <- "RMSE"
rownames(df) <- c('Holt-Winters','Holt-Winters damped','Naive','ETS')
datatable(df) %>% formatRound('RMSE',3)