Data 624 Homework 5
Anthony Pagan
3/7/2020
1
Consider the pigs series — the number of pigs slaughtered in Victoria each month.
(a)
Use the ses() function in R to and the optimal values of @ and lo , and generate forecasts for the next four months.
## 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.8
## Jan 1996 98816.41 83448.52 114184.3 75313.25 122319.6
## Feb 1996 98816.41 82955.06 114677.8 74558.56 123074.2
## Mar 1996 98816.41 82476.49 115156.3 73826.66 123806.2
## Apr 1996 98816.41 82011.54 115621.3 73115.58 124517.2
## May 1996 98816.41 81559.12 116073.7 72423.66 125209.2
## Jun 1996 98816.41 81118.26 116514.6 71749.42 125883.4
(b)
Compute a 95% prediction interval for the first forecast using y ± 1.96s where is the standard deviation of the residuals. Compare your interval with the interval produced by R.
## Lower95 Upper95
## 78679.97 118952.84## Lower95 Upper95
## 78611.6 119021.25
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.
(b)
Use the ses() function to forecast each series, and plot the forecasts.
(c)
Compute the RMSE values for the training data in each case.
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303 -0.2117522## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887 -0.14177636
(a)
Now apply Holt’s linear method to the paperback and hardback series and compute four-day forecasts in each case.
(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.
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7.175981 33.63769 27.8431 0.4736071 15.57784 0.7021303 -0.2117522## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 9.166735 31.93101 26.77319 2.636189 13.39487 0.7987887 -0.1417763## 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## 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(c)
Compare the forecasts for the two series using both methods. Which do you think is best?
The simple exponential smoothing has a higher RMSE value compare to the holts method. The books data has a clear trend and the holts method was specifically built to expend simple exponential smoothing to allow forecasting for data with trends.
(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 .
## Lower95 Upper95
## 148.4384 270.4951## Lower95 Upper95
## 156.1674 262.76627
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?
The holt method with damped option has a lower RSME and therefore is the best model.
## ME RMSE MAE MPE MAPE MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
## ACF1
## Training set -0.003195358## ME RMSE MAE MPE MAPE MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
## ACF1
## Training set 0.01348202## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.306667 26.97157 20.27226 -2.497352 10.42787 1 -0.1466138
Recall your retail time series data (from Exercise 3 in Section 2.10).
(a)
Why is multiplicative seasonality necessary for this series?
The increasing trend tells us that the multiplicative series is necessary.
(b)
Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
(c)
Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
The lower RMSE of the multiplicative method tells us that it would be the prefered method. In fact , all multiplicative accuracy values are smaller and are preferred.
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2585256 15.37847 10.58919 -0.1132325 5.464699 0.5592564
## ACF1
## Training set 0.2405126## ME RMSE MAE MPE MAPE MASE
## Training set 0.1170648 13.29378 8.991856 -0.1217735 3.918351 0.4748948
## ACF1
## Training set 0.08635577## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.892255 15.70405 10.81895 0.7781595 5.563244 0.5713906 0.2174256## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.414869 13.30494 9.042151 0.6105987 3.959617 0.4775511 0.04077895(d)
Check that the residuals from the best method look like white noise.
The residuals for the multiplicative method border are probably not white noise as 2 bars are past the boundaries and 2 more touching or are barely past the boundaries.
(e)
Now end 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?
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4504682 8.655091 6.242194 0.2405149 3.136605 0.3936879
## ACF1
## Training set 0.06874494## ME RMSE MAE MPE MAPE MASE
## Training set 0.4785645 8.873172 6.442556 0.3066717 3.256778 0.4063245
## ACF1
## Training set -0.009321209## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 8.458154 20.1662 15.85569 4.997393 8.130524 1 0.718433## ETS(M,Ad,M)
##
## Call:
## ets(y = myts2, model = "ZZZ", damped = NULL, alpha = NULL, beta = NULL,
##
## Call:
## gamma = NULL, phi = NULL, additive.only = FALSE, lambda = NULL,
##
## Call:
## biasadj = FALSE, restrict = TRUE, allow.multiplicative.trend = FALSE)
##
## Smoothing parameters:
## alpha = 0.6052
## beta = 0.0071
## gamma = 1e-04
## phi = 0.98
##
## Initial states:
## l = 62.7069
## b = 0.803
## s = 0.9419 0.9066 0.9363 1.5171 1.0882 1.0126
## 0.9654 0.9422 0.9193 0.8925 0.9601 0.9179
##
## sigma: 0.0419
##
## AIC AICc BIC
## 3346.794 3348.945 3415.556
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4504682 8.655091 6.242194 0.2405149 3.136605 0.3936879
## ACF1
## Training set 0.068744949
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?
The ETS on the BoxCox transformed data has a lower RMSE and is better forecast
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4504682 8.655091 6.242194 0.2405149 3.136605 0.3936879
## ACF1
## Training set 0.06874494## ME RMSE MAE MPE MAPE MASE
## Training set 0.4785645 8.873172 6.442556 0.3066717 3.256778 0.4063245
## ACF1
## Training set -0.009321209## ME RMSE MAE MPE MAPE MASE
## Training set -0.5887334 8.629737 6.260496 -0.3091327 3.127894 0.3948422
## ACF1
## Training set 0.09273547## ETS(A,A,A)
##
## Call:
## ets(y = myts2, model = "ZZZ", damped = NULL, alpha = NULL, beta = NULL,
##
## Call:
## gamma = NULL, phi = NULL, additive.only = FALSE, lambda = lambda,
##
## Call:
## biasadj = FALSE, restrict = TRUE, allow.multiplicative.trend = FALSE)
##
## Box-Cox transformation: lambda= 0.084
##
## Smoothing parameters:
## alpha = 0.5782
## beta = 0.0068
## gamma = 2e-04
##
## Initial states:
## l = 4.9688
## b = 0.0103
## s = -0.0807 -0.1428 -0.0932 0.6708 0.1525 0.0421
## -0.0367 -0.0743 -0.1136 -0.1595 -0.0469 -0.1178
##
## sigma: 0.0637
##
## AIC AICc BIC
## 122.7920 124.7105 187.7334
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.5887334 8.629737 6.260496 -0.3091327 3.127894 0.3948422
## ACF1
## Training set 0.09273547APPENDIX
Code used in analysis
knitr::opts_chunk$set(
echo = FALSE,
message = FALSE,
warning = FALSE
)
#knitr::opts_chunk$set(echo = TRUE)
require(knitr)
library(ggplot2)
library(tidyr)
library(MASS)
library(psych)
library(kableExtra)
library(dplyr)
library(faraway)
library(gridExtra)
library(reshape2)
library(leaps)
library(pROC)
library(caret)
library(naniar)
library(pander)
library(pROC)
library(mlbench)
library(e1071)
library(fpp2)
data(pigs)
ses(pigs)
pc<-ses(pigs, h=4)
autoplot(pc)+
autolayer(fitted(pc), series="Fitted")
pcsd<-sqrt(var(residuals(pc)))
p<-head(pc[["mean"]],1)
c(Lower95 = p - 1.96 * pcsd,
Upper95 = p + 1.96 * pcsd)
pcsd <- sqrt(sum(residuals(pc)^2)/(length(pigs)-2))
c(Lower95 = p - 1.96 * pcsd,
Upper95 = p+ 1.96 * pcsd)
data(books,)
autoplot(books)
pb<-c(books[,'Paperback'])
hc<-c(books[,'Hardcover'])
bkspb<-ses(pb,h=4)
bkshc<-ses(hc,h=4)
autoplot(books)+
autolayer(bkspb, series="Paperback", PI=FALSE)+
autolayer(bkshc, series="Hardcover", PI=FALSE)
accuracy(bkspb)
accuracy(bkshc)
pba<-holt(pb, h=4)
hca<-holt(hc,h=4)
autoplot(books)+
autolayer(pba, series="Paperback", PI=FALSE)+
autolayer(hca, series="Hardcover", PI=FALSE)
accuracy(bkspb)
accuracy(bkshc)
accuracy(pba)
accuracy(hca)
pbasd<-sqrt(pba$model$mse)
c(Lower95 = pba$mean[1] - 1.96 * pbasd,
Upper95 = pba$mean[1] + 1.96 * pbasd)
hcasd<-sqrt(hca$model$mse)
c(Lower95 = pba$mean[1] - 1.96 * hcasd,
Upper95 = pba$mean[1] + 1.96 * hcasd)
autoplot(eggs)
ht100<-holt(eggs, h=100)
ht100d<-holt(eggs, h=100, damped = TRUE)
autoplot(ht100)+
autolayer(ht100, series = "Holts Method", PI=FALSE)+
autolayer(ht100d, series = "Damped Holts method" , PI=FALSE)+
ggtitle("Eggs")+ xlab("Year")+
ylab("Eggs")+
guides(colour=guide_legend(title="Forecast"))
lambda<-BoxCox.lambda(eggs)
egg100<-rwf(eggs,h=100, lambda=lambda)
autoplot(BoxCox(eggs, lambda))+
autolayer(egg100, series="egg100")
accuracy(ht100d)
accuracy(ht100)
accuracy(egg100)
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
myts <- ts(retaildata[,"A3349873A"],
frequency=12, start=c(1982,4))
checkresiduals(myts)
myts %>% decompose(type="multiplicative")%>%
autoplot(frequency = frequency(retaildata)) + xlab("Year") +
ggtitle("Classical multiplicative decomposition of retail datas")
fit1 <- hw(myts,seasonal="additive")
fit2 <- hw(myts,seasonal="multiplicative")
autoplot(myts) +
autolayer(fit1, series="HW additive forecasts", PI=FALSE) +
autolayer(fit2, series="HW multiplicative forecasts",
PI=FALSE) +
xlab("Year") +
ylab("Retail Sales") +
ggtitle("Retail Data") +
guides(colour=guide_legend(title="Forecast"))
fit3 <- hw(myts,seasonal="additive", damped = TRUE)
fit4 <- hw(myts,seasonal="multiplicative", damped = TRUE)
autoplot(myts) +
autolayer(fit3, series="HW additive forecasts", PI=FALSE) +
autolayer(fit4, series="HW multiplicative forecasts",
PI=FALSE) +
xlab("Year") +
ylab("Retail Sales") +
ggtitle("Retail Data Damped") +
guides(colour=guide_legend(title="Forecast Damped"))
accuracy(fit1,h=1)
accuracy(fit2, h=1)
accuracy(fit3, h=1)
accuracy(fit4, h=1)
ggAcf(residuals(fit2))
myts2 <- ts(retaildata[,"A3349873A"],
frequency=12, start=c(1982,4), end=c(2010,4))
emyts2<-ets(myts2, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL,
gamma=NULL, phi=NULL, lambda=NULL, biasadj=FALSE,
additive.only=FALSE, restrict=TRUE,
allow.multiplicative.trend=FALSE)
emyts2 %>% forecast(h=24) %>%
autoplot() +
ylab("Retail Sales")
fit2 <- hw(myts2,seasonal="multiplicative")
autoplot(fit2)
cbind('Residuals' = residuals(emyts2), 'Forecast errors' =residuals(emyts2, type='response')) %>%
autoplot(facet=TRUE)+ xlab("Year")+ylab("")
accuracy(emyts2)
accuracy(fit2)
accuracy(snaive(myts2))
summary(emyts2)
lambda<-BoxCox.lambda(myts2)
bmyts<-BoxCox(myts2, lambda)
ebmyts<-ets(myts2, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL,
gamma=NULL, phi=NULL, lambda=lambda, biasadj=FALSE,
additive.only=FALSE, restrict=TRUE,
allow.multiplicative.trend=FALSE)
ebmyts %>% forecast(h=24) %>%
autoplot() +
ylab("Retail Sales")
accuracy(emyts2)
accuracy(fit2)
accuracy(ebmyts)
summary(ebmyts)