shampoo_ts
2024-04-25
#install required packages
#Q1. ARIMA – shampoo
#We have a dataset, which includes monthly sales for shampoo at a retailer store from January 2001 to December 2003 (shampoo.csv). Each observation includes two values: month-year pair, and sales in that particular month. First we read the csv file.
data <- read.csv(file = "C:/Documents and Settings/mt/OneDrive - Erasmus University Rotterdam/Desktop/CLASSES/bam_sca/2024/assignments/shampoo.csv",header=TRUE)
#We then convert sales into a time series object. We specify frequency as 12 because it is a monthly data.
#frequency indicates the number of observations per unit of time. It is typically 7 for daily data (the natural time period is a week), and 12 for monthly data (the natural time period is a year).
shampoo <- ts(data[, 2], frequency = 12, start = c(2001, 1))
print(shampoo,calendar=TRUE)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2001 266.0 145.9 183.1 119.3 180.3 168.5 231.8 224.5 192.8 122.9 336.5 185.9
## 2002 194.3 149.5 210.1 273.3 191.4 287.0 226.0 303.6 289.9 421.6 264.5 342.3
## 2003 339.7 440.4 315.9 439.3 401.3 437.4 575.5 407.6 682.0 475.3 581.3 646.9#some useful function to get insights about the data set
summary(shampoo)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 119.3 192.4 280.1 312.6 411.1 682.0head(shampoo)## Jan Feb Mar Apr May Jun
## 2001 266.0 145.9 183.1 119.3 180.3 168.5#a. With time series analysis, the first step is always to visually inspect time series. A couple of functions are readily available in R.
plot.ts(shampoo)
autoplot(shampoo)
#alternatively, ggtsdisplay() plots acf and pacf along with the time series plot. But note that the graph of acf does not start from lag 0.
ggtsdisplay(shampoo)
#Judging from the plots, it is clear that there is an upward trend in the sales, but there is no clear evidence of seasonality.#b. Is the time series stationary? If not, what needs to be done to make it stationary?
Due to the trend in time series, it is non-stationary. So we need to stationarize time series and we can get rid of trend by taking first-order difference.
#diff() function returns lagged differences and to take the first-order difference, we should set differences=1
shampoo.diff1 <- diff(shampoo, differences = 1)
autoplot(shampoo.diff1)
ggtsdisplay(shampoo.diff1)
#We plot the time series after the difference: it looks stationary.
##It has a random band scattered around point 0, variance is sort of constant, there is no consistent pattern, no parts of data is below zero or above zero. There does not seem to be trend or seasonality. but we are not sure! Then what can we do?
#Stationary tests...
#We run ADF, PP and KPSS tests to formally test the stationarity of time series after taking the first-order difference
###stationarity tests
adf.test(shampoo.diff1)## Warning in adf.test(shampoo.diff1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: shampoo.diff1
## Dickey-Fuller = -5.9245, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationarypp.test(shampoo.diff1)## Warning in pp.test(shampoo.diff1): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: shampoo.diff1
## Dickey-Fuller Z(alpha) = -54.269, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationarykpss.test(shampoo.diff1)## Warning in kpss.test(shampoo.diff1): p-value greater than printed p-value##
## KPSS Test for Level Stationarity
##
## data: shampoo.diff1
## KPSS Level = 0.32152, Truncation lag parameter = 3, p-value = 0.1####Conclusion: all the tests suggest that the time series is stationary.#There are also two automatic functions, ndiffs() and nsdiffs() that would tell us how many first-order differences (what is d?), and how many seasonal differences, (what is D?) respectively, that we need to take to make the time series stationary.
#first order difference
ndiffs(shampoo)## [1] 1#seasonal stationarity
nsdiffs(shampoo)## [1] 0#Test results suggest that we only need to take first-order difference once to obtain a stationary time series.#c. Once we have a stationary time series, the next step is to determine the optimal orders of MA and AR components. We first plot the ACF and PACF of the time series.
### determine optimal p and q
# acf plot
acf(shampoo.diff1, lag.max = 20) # q <= 2 
#shows the same ACF plot but starting from the lag 1 instead of 0
ggAcf(shampoo.diff1)
#Based on the plots, we know that the order of MA is not greater than 2 (as the first two lags of ACF are significant)
#pacf plot
pacf(shampoo.diff1, lag.max = 20) # p <= 1
ggPacf(shampoo.diff1)
#the order of AR is not greater than 1 (as only the first lag of PACF is significant).
#no seasonal factors are detected, perhaps very slight, rather close to blue band for ACF at lag 12.
#or, in a neater way
ggtsdisplay(shampoo.diff1)
#d. Estimation: Use auto.arima() to evaluate the potential models based on BIC criteria. \(auto.arima()\) evaluates the quality of various ARIMA models. Potentially, we can use our earlier analysis as input: according to maximum values for p and q, the order constraints can be set. But for this simple exercise, we just use the default setting from \(auto.arima()\), and run a stepwise search.
# choose optimal p and q based on information criteria
auto.arima(shampoo, trace = TRUE, ic = 'bic') ##
## ARIMA(2,1,2)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0) with drift : 433.9838
## ARIMA(1,1,0)(1,0,0)[12] with drift : 410.9315
## ARIMA(0,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(0,1,0) : 430.7841
## ARIMA(1,1,0) with drift : 410.6883
## ARIMA(1,1,0)(0,0,1)[12] with drift : 409.1295
## ARIMA(1,1,0)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(2,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] with drift : 407.3023
## ARIMA(1,1,1) with drift : 408.688
## ARIMA(1,1,1)(1,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(1,0,0)[12] with drift : 408.7365
## ARIMA(2,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,2)(0,0,1)[12] with drift : 410.1852
## ARIMA(0,1,2)(0,0,1)[12] with drift : Inf
## ARIMA(2,1,2)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] : 413.8951
##
## Best model: ARIMA(1,1,1)(0,0,1)[12] with drift## Series: shampoo
## ARIMA(1,1,1)(0,0,1)[12] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## -0.5479 -0.4823 -0.7441 12.8027
## s.e. 0.1736 0.1753 0.6477 2.4503
##
## sigma^2 = 3435: log likelihood = -194.76
## AIC=399.53 AICc=401.59 BIC=407.3###gives the same result even if you specify the no. of difference:
auto.arima(shampoo, trace = TRUE, d=1, ic = 'bic') ##
## ARIMA(2,1,2)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0) with drift : 433.9838
## ARIMA(1,1,0)(1,0,0)[12] with drift : 410.9315
## ARIMA(0,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(0,1,0) : 430.7841
## ARIMA(1,1,0) with drift : 410.6883
## ARIMA(1,1,0)(0,0,1)[12] with drift : 409.1295
## ARIMA(1,1,0)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(2,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] with drift : 407.3023
## ARIMA(1,1,1) with drift : 408.688
## ARIMA(1,1,1)(1,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(1,0,0)[12] with drift : 408.7365
## ARIMA(2,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,2)(0,0,1)[12] with drift : 410.1852
## ARIMA(0,1,2)(0,0,1)[12] with drift : Inf
## ARIMA(2,1,2)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] : 413.8951
##
## Best model: ARIMA(1,1,1)(0,0,1)[12] with drift## Series: shampoo
## ARIMA(1,1,1)(0,0,1)[12] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## -0.5479 -0.4823 -0.7441 12.8027
## s.e. 0.1736 0.1753 0.6477 2.4503
##
## sigma^2 = 3435: log likelihood = -194.76
## AIC=399.53 AICc=401.59 BIC=407.3auto.arima(shampoo, trace = TRUE, d=1, max.p=1 ,max.q=2,ic = 'bic') ##
## ARIMA(1,1,2)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0) with drift : 433.9838
## ARIMA(1,1,0)(1,0,0)[12] with drift : 410.9315
## ARIMA(0,1,1)(0,0,1)[12] with drift : Inf
## ARIMA(0,1,0) : 430.7841
## ARIMA(1,1,0) with drift : 410.6883
## ARIMA(1,1,0)(0,0,1)[12] with drift : 409.1295
## ARIMA(1,1,0)(1,0,1)[12] with drift : Inf
## ARIMA(0,1,0)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] with drift : 407.3023
## ARIMA(1,1,1) with drift : 408.688
## ARIMA(1,1,1)(1,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(1,0,0)[12] with drift : 408.7365
## ARIMA(1,1,2)(0,0,1)[12] with drift : 410.1852
## ARIMA(0,1,2)(0,0,1)[12] with drift : Inf
## ARIMA(1,1,1)(0,0,1)[12] : 413.8951
##
## Best model: ARIMA(1,1,1)(0,0,1)[12] with drift## Series: shampoo
## ARIMA(1,1,1)(0,0,1)[12] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## -0.5479 -0.4823 -0.7441 12.8027
## s.e. 0.1736 0.1753 0.6477 2.4503
##
## sigma^2 = 3435: log likelihood = -194.76
## AIC=399.53 AICc=401.59 BIC=407.3#Based on the output of $auto.arima()$, a couple of models have similar BICs. Suppose that we choose the two models with the lowest BICs, as the candidate models that we would like to evaluate further.
# Best model: ARIMA(1,1,1)(0,0,1)[12] with drift (BIC=407.3), the coefficients of the model are also provided in the output.
# Second best: ARIMA(1,1,1) with drift (BIC=408.688)
#Note: infinity means R function tried to evaluate this model and it is close to nonstationary model, so basically R suggests not to choose it.
#Both outcomes give drift term, this is the slope of a linear line.#In sample forecast performance – training set performance The only difference in the two models is the seasonal MA part. Recall that our preliminary analysis with ACF and PACF plots does not reveal seasonal factors, so we will have to determine whether or not to include the seasonal MA part in the model.
# two candidate models
shampoo.m1 <- Arima(shampoo, order = c(1, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 12), include.drift = TRUE)
shampoo.m2 <- Arima(shampoo, order = c(1, 1, 1), include.drift = TRUE)
coef(shampoo.m2)## ar1 ma1 drift
## -0.5673814 -0.5132668 12.0965205#A couple of functions are proved to be useful for us to evaluate the in-sample performance/fit of the model. One is $accuracy()$ function, which summarizes various measures of fitting errors.
#in-sample one-step forecast performance
accuracy(shampoo.m1)## ME RMSE MAE MPE MAPE MASE
## Training set -9.343786 54.38305 41.40145 -10.19758 17.43117 0.269709
## ACF1
## Training set -0.01143663accuracy(shampoo.m2)## ME RMSE MAE MPE MAPE MASE
## Training set -3.690698 65.69307 50.94426 -10.03182 21.12789 0.3318754
## ACF1
## Training set 0.01636144###m1 model, including seasonal MA part seems to perform slightly better
###ARIMAX####
#adding covariates
data$year<-c(rep(2001,12),rep(2002,12),rep(2003,12))
#1 being January
data$month<-rep(1:12,3)
#add year to shampoo data specifically
shampoo.m3year<-Arima(shampoo, order = c(1, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 12), include.drift = TRUE, xreg=data$year)
accuracy(shampoo.m3year)## ME RMSE MAE MPE MAPE MASE
## Training set -6.601165 50.66932 39.98636 -8.665808 16.72924 0.2604904
## ACF1
## Training set -0.1129301shampoo.m3month<-Arima(shampoo, order = c(1, 1, 1),
seasonal = list(order = c(0, 0, 1), period = 12), include.drift = TRUE, xreg=data$month)
accuracy(shampoo.m3month)## ME RMSE MAE MPE MAPE MASE
## Training set -8.967601 48.98702 37.88462 -9.592053 15.95265 0.2467987
## ACF1
## Training set -0.07381429#which one is the best model among the 4?
#we can plot what is the model prediction and how the fit compares to the provided data set
plot(shampoo,type="l")
lines(fitted(shampoo.m1),col="red")
lines(fitted(shampoo.m3month),col="green")
#fitted gives the predictions of one-step ahead forecast; we can manually calculate the error by taking the difference of fitted values from the observed values
accuracy(shampoo.m1)## ME RMSE MAE MPE MAPE MASE
## Training set -9.343786 54.38305 41.40145 -10.19758 17.43117 0.269709
## ACF1
## Training set -0.01143663#mean errors
mean(shampoo-fitted(shampoo.m1))## [1] -9.343786#e. Validation: Perform a residual analysis of the best model chosen. Are there any warning signals? In the post-estimation analysis, we should also check out the residual plots, including residuals’ time series plot and ACFs to make sure that there is no warning message. In particular, residuals should have a zero mean, constant variance, and should be distributed symmetrically around mean zero. ACF of any lag greater than zero should be statistically insignificant.
# residual analysis
autoplot(shampoo.m1$residuals)
ggAcf(shampoo.m1$residuals)
checkresiduals(shampoo.m1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(0,0,1)[12] with drift
## Q* = 5.1294, df = 4, p-value = 0.2743
##
## Model df: 3. Total lags used: 7#Residuals appear to be uncorrelated
#have zero mean, constant variance
#histogram of residuals mimics the std. normal distr. So there is no warning message.
#Ljung Box test: more rigorous way of testing residuals whether coming from white-noise process.
#H0 = residuals exhibits no autocorrelation for a fixed number of lags L
#Since p value is large (not close to 0.05), we fail to reject the null.
#also check the models with covariates
checkresiduals(shampoo.m3month)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(1,1,1)(0,0,1)[12] errors
## Q* = 6.213, df = 4, p-value = 0.1838
##
## Model df: 3. Total lags used: 7checkresiduals(shampoo.m3year)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(1,1,1)(0,0,1)[12] errors
## Q* = 4.4207, df = 4, p-value = 0.3521
##
## Model df: 3. Total lags used: 7#Forecasting future observations: Once we have selected a model, forecasting part is easy - it is implemented by \(forecast()\) function from the forecast package.
# h-step ahead (in this case, 6-month ahead) forecast
shampoo.f <- forecast(shampoo.m1, h = 6)
autoplot(shampoo.f)
#Out of sample performance Lastly, let’s discuss how we can evaluate the true performance of ARIMA models on a real-life setting (out-of-sample performance). Since we haven’t observed the future yet, we can not measure the performance with those unobserved data points! Instead, what we can do closest to that is, we divide the data into 2 sets: training set and test (validation) set. For time series analysis, we divide data based on time indices of observations. Earlier observations are used for training, and more recent observations are used for testing This is due to intertemporal correlations in time series data, so we don’t randomly sample them!
We have 36 months of observations = 3 years of data –> 80% corresponds to nearly 2.5 years of data to be allocated to training. Training: Jan 2001 up to June 2003 Test: July 2003 up to Dec 2003 Train the model on 2.5 year data: Based on \(auto.arima()\), we choose two candidate models with the lowest BICs (here, we are skipping all the steps on stationarity tests, ACF and PACF plots etc.; but they are definitely important parts of the analysis).
### model evaluation
# Fit model with the first 2.5-year of data
auto.arima(window(shampoo, end = c(2003, 6)), d = 1, trace = TRUE, ic = 'bic')##
## ARIMA(2,1,2) with drift : 337.6606
## ARIMA(0,1,0) with drift : 349.1833
## ARIMA(1,1,0) with drift : 333.9026
## ARIMA(0,1,1) with drift : 332.9748
## ARIMA(0,1,0) : 345.9445
## ARIMA(1,1,1) with drift : 332.667
## ARIMA(2,1,1) with drift : 335.9243
## ARIMA(1,1,2) with drift : 334.3087
## ARIMA(0,1,2) with drift : 332.4953
## ARIMA(0,1,3) with drift : 335.848
## ARIMA(1,1,3) with drift : 337.6186
## ARIMA(0,1,2) : 332.9803
##
## Best model: ARIMA(0,1,2) with drift## Series: window(shampoo, end = c(2003, 6))
## ARIMA(0,1,2) with drift
##
## Coefficients:
## ma1 ma2 drift
## -1.1532 0.4352 8.6419
## s.e. 0.1951 0.1991 3.1521
##
## sigma^2 = 3722: log likelihood = -159.51
## AIC=327.03 AICc=328.69 BIC=332.5# two candidate models:
#best model: ARIMA(0,1,2) with drift, BIC: 332.5,
#the second best: ARIMA(1,1,1) with drift, BIC:332.6
#To understand the performance of each model, we fit these candidate models separately with Arima on the training set
m1 <- Arima(window(shampoo, end = c(2003, 6)), order = c(1, 1, 1), include.drift = TRUE)
m2 <- Arima(window(shampoo, end = c(2003, 6)), order = c(0, 1, 2), include.drift = TRUE)##How to evaluate out-of-sample performance of the candidates’ models? Two different measures to evaluate out-of-sample performance of models: 1. based on one-step ahead forecast, 2. based on multi-step ahead forecast. With one-step ahead forecast, at time period \(t\) (\(t\) is from June 2003 to Nov 2003), we have all the information up to time \(t\), and generate a forecast of sales in time period \(t+1\). ##One step ahead forecast
# Apply fitted model m1 and m2 to later (test) data: July 2003 up to Dec 2003
m1.f <- Arima(window(shampoo, start = c(2003, 7)), model = m1)
m2.f <- Arima(window(shampoo, start = c(2003, 7)), model = m2)
#$accuracy()$ function summarizes various measures of one-step ahead forecasting errors.
###out-of-sample one-step ahead forecasts
accuracy(m1.f)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 11.71466 79.46943 61.1882 -0.3220726 11.25665 NaN -0.4886267accuracy(m2.f)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 10.23239 77.45575 60.33499 -0.6108842 11.22599 NaN -0.5231363#Conclusion: in almost all errors the second model seems to outperform the first model.##Alternatively, we can apply multi-step ahead forecast. At time June 2003, we generate forecasts of monthly sales for each month from July 2003 to Dec 2003. Then again, we use \(accuracy()\) to report various measures of forecasting errors across this 6-month period.
# out-of-sample multi-step ahead forecasts
accuracy(forecast(m1, h = 6), window(shampoo, start = c(2003, 7)))## ME RMSE MAE MPE MAPE MASE
## Training set -5.309904 57.14505 45.87963 -9.950254 21.76062 0.3817294
## Test set 114.664852 147.28821 124.98553 18.004335 20.53640 1.0399092
## ACF1 Theil's U
## Training set 0.08956709 NA
## Test set -0.70731992 0.8573852accuracy(forecast(m2, h = 6), window(shampoo, start = c(2003, 7)))## ME RMSE MAE MPE MAPE MASE
## Training set -4.677904 56.79416 45.75745 -9.574285 21.79676 0.3807128
## Test set 108.864315 142.05917 120.18644 16.984726 19.76248 0.9999796
## ACF1 Theil's U
## Training set 0.09609745 NA
## Test set -0.71161499 0.824693#Conclusion: again the second model seems to perform better both in test and training set
m1.f6<-as.data.frame(forecast(m1,h=6))
mean(shampoo[31:36]-m1.f6$`Point Forecast`)## [1] 114.6649In practice whether you want to evaluate the performance of models using one-step or multi-step ahead forecasts depends on the application. – If your firm has a very responsive SC, and is able to source products quickly (say less than a month), then one-step ahead forecast would be sufficient. – However, in many cases, firms have only one chance to place order, and need to determine order quantities for the entire selling season well in advance. (say, for fashion companies, you want to forecast demand for the whole season e.g. the whole winter, from Oct to Feb). In this case, firms prefer to minimize forecasting errors for the entire selling season, at the time when the order is placed. So, a model with lowest multi-period ahead forecasting errors is preferred.
#Q2. Holt-Winter’s model For time series analysis, the first step is always to visually inspect the time series. In this regard, the \(stl()\) function is quite useful.
#pipeline operator
shampoo %>% stl(s.window = "periodic") %>% autoplot
It decomposes the original time series into trend, seasonal factors, and random error terms. Linear increasing trend and yearly seasonality.
How to interpret the relative importance of different components? – Grey bars
For this data set, the grey bar of the trend panel is only slightly larger than that of the original time series panel –> trend component contributes to a great proportion of variations in the time series.
On the other hand, the grey bar of the seasonal panel is very large, even larger than the grey bar of random error term, which indicates that the contribution of seasonal factors to the variation in the original time series is marginal –> no seasonality in the data.
##a. Model estimation, using Holt-Winters
#We first fit the data with $HoltWinters()$ function. As we observe no seasonality in the data, we specify $gamma = F$ in the function;
#smoothing parameters
#alpha: level
#beta: trend
#gamma: seasonal component
shampoo.HW <- HoltWinters(shampoo, gamma = FALSE)
shampoo.HW## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = shampoo, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.3563906
## beta : 1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 612.67141
## b 23.92159#lower value of alpha/beta/gamma means it is less responsive to recent observations.
#values near 1 means the latest value has more weight (beta=1 means trend component reacts more quickly to changes in the data)
#we also can see the estimated level and trend components given by coefficients a and b.
#Fitting the model with $HoltWinters()$ is straightforward, but one thing we need to keep in mind is that, as in any optimization problem, results from $HoltWinters()$ are sensitive to starting values.
shampoo.HW2 <- HoltWinters(shampoo, gamma = FALSE, optim.start = c(alpha = 0, beta = 0), l.start = 200, b.start = -2.6)
shampoo.HW2## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = shampoo, gamma = FALSE, l.start = 200, b.start = -2.6, optim.start = c(alpha = 0, beta = 0))
##
## Smoothing parameters:
## alpha: 0.05794248
## beta : 1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 598.15635
## b 27.91208#alpha: smoothing parameter for level
#beta: smoothing parameter for trend
#l.start=Start value for level (a[0]).
#b.start=Start value for trend (b[0])
#plot of both fitted models
plot(shampoo.HW)
#shows data in black and fitted values in red
##on top of it plot the fitted values from the second model in blue
lines(fitted(shampoo.HW2)[,2], col = "blue", lty = 2)
##accuracy() does not work for Holt-Winters fnc.!!!
#in-sample one-step forecast: standard dev. of errors --> variance=SSE/(n-2)
#RMSE
sqrt(shampoo.HW$SSE/(length(shampoo)-2))## [1] 95.95459sqrt(shampoo.HW2$SSE/(length(shampoo)-2))## [1] 67.11598We estimate two models - one with the default starting values and one with manually supplied values - and we notice that the fit differs significantly. –> S.E. reduces by around one-third after changing the initial values.
In practice, how to guess the “correct” starting points?
Not easy, but luckily this problem is solved by \(ets()\).
With \(ets()\), initial states and smoothing parameters are jointly estimated by maximizing the likelihood function. We need to specify the model in \(ets()\) using three letters: the error type, the trend type, and the seasonal variation type.
Steps to follow: (1) check out time series plot, and see if there is any trend and seasonality, and whether they are additive (linear trend, constant variations in seasonal factors) or multiplicative; (2) run \(ets()\) with model = “ZZZ”, and see whether the best model is consistent with your expectation; (3) if they are consistent, it gives us confidence that our model specification is correct; otherwise try to figure out why there is a discrepancy.
In this problem: error terms are also most of the time additive, the trend seems to be linear, which suggests an additive pattern, and we have no seasonality.
##b. model est. using ets()
shampoo.ets <- ets(shampoo, model = "AAN")
shampoo.ets## ETS(A,A,N)
##
## Call:
## ets(y = shampoo, model = "AAN")
##
## Smoothing parameters:
## alpha = 0.0563
## beta = 0.0563
##
## Initial states:
## l = 176.6845
## b = -1.237
##
## sigma: 70.3984
##
## AIC AICc BIC
## 441.0668 443.0668 448.9844shampoo.ets2 <- ets(shampoo, model = "ZZZ")
#Observe that R automatically comes up with the same model to ours.
shampoo.ets2## ETS(A,A,N)
##
## Call:
## ets(y = shampoo, model = "ZZZ")
##
## Smoothing parameters:
## alpha = 0.0563
## beta = 0.0563
##
## Initial states:
## l = 176.6845
## b = -1.237
##
## sigma: 70.3984
##
## AIC AICc BIC
## 441.0668 443.0668 448.9844After estimation, we can use \(accuracy()\) function to determine in-sample fit and \(forecast()\) function to generate forecast.
Let’s compare HoltWinter’s and ets function! #accuracy and forecast
# in-sample one-step forecast errors
accuracy(shampoo.ets)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 14.84085 66.37224 50.5648 0.2848111 17.5656 0.3294034 -0.4839511#6-step ahead forecast
shampoo.HW.f <- forecast(shampoo.HW, h = 6)
shampoo.ets.f <- forecast(shampoo.ets, h = 6)
#We plot forecasts from $HoltWinters()$ in blue and $ets()$ in red on the same plot...
plot(shampoo.HW.f)
lines(fitted(shampoo.HW.f), col = "blue")
lines(fitted(shampoo.ets), col = "red", lty = 2)
lines(shampoo.ets.f$mean, col = "red", lty = 2)
and notice that, even though the in-sample fit from the two functions differ significantly, the forecasting results are sort of comparable.
Above is an illustration of \(HoltWinters()\) and \(ets()\) functions.
In practice, we need to first divide samples into training and test sets, and then evaluate models by comparing out-of-sample performance. Training: use the first 2.5 years of data for training Test: the remaining 0.5 year for evaluating the performance of models from \(HoltWinters()\) and \(ets()\).
#out-of-sample performance
### model evaluation
# Fit model with first 2.5 years of data
m1 <- HoltWinters(window(shampoo, end = c(2003, 6)), gamma = FALSE)
m2 <- ets(window(shampoo, end = c(2003, 6)), model = "AAN")
# out-of-sample multi-step forecasts
accuracy(forecast(m1, h = 6), window(shampoo, start = c(2003, 7)))## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 12.56592 87.02346 72.12966 2.647523 32.75038 0.6001358 -0.1434291
## Test set 62.15626 108.21178 94.39078 8.636494 16.15520 0.7853536 -0.7637161
## Theil's U
## Training set NA
## Test set 0.6349905accuracy(forecast(m2, h = 6), window(shampoo, start = c(2003, 7)))## ME RMSE MAE MPE MAPE MASE
## Training set 12.07206 58.35153 43.79584 0.5196929 17.89491 0.3643918
## Test set 54.47860 104.04027 92.37836 7.1909678 16.01970 0.7686098
## ACF1 Theil's U
## Training set -0.3082067 NA
## Test set -0.7591104 0.6137215The results show that the model from \(ets()\) is consistently better across all different measures.