Ch.5: Q.2,5,8

Q.2 \(\textit{Relationship between Moving Average and Exponential Smoothing}\)

The relationship between the moving average window \(w\) and the smoothing constant \(\alpha\) is given by \[ w = \frac{2}{\alpha} - 1 \Rightarrow \alpha = \frac{2}{w+1}. \] Using a very short window for a moving average means we are using very recent information for forecasting purposes. For example a moving average window \(w=1\) yields a naive forecast. The smoothing constant would then be \(\alpha = \frac{2}{1+1} = 1\). So to achieve an equivalent result using simple exponential smoothing, we should use \(\alpha\) close to 1.

Q.5 \(\textit{Forecasting Department Store Sales}\)

  1. A plot of the series is reproduced below:
library(forecast)
# read in the data
store <-  read.csv("DeptStoreSales.csv")
# create time series with plot
salesTS <- ts(store$Sales,start=c(1,1),frequency=4)
yrange = range(salesTS)
plot(c(1,7),yrange,type="n",xlab="Years",ylab="Sales (thousands)",bty="l",xaxt="n",yaxt="n")
axis(1, at=seq(1,7,1), labels=format(seq(1,7,1)))
axis(2, at=seq(50000,110000,10000), labels=format(seq(50,110,10)), las=2)
lines(salesTS,bty="l")
title("Department Store Sales")

\(\bullet\) Moving average of raw series

Seasonality is easily observed as evidenced by the periodic peaks. Trend should also be assumed to be present due to the somewhat upward linear trend. Using a moving average of the raw series does not account for either seasonality or trend and therefore would not be suitable for forecasting this series.

\(\bullet\) Moving average of deseasonalized series

Trend is still present so for the same reason as above this method would not be suitable.

\(\bullet\) Simple exponential smoothing of the raw series

Similar to a moving average of a raw series that contains trend and seasonality, simple exponential smoothing of the raw series should not be used because this method should only be used when noise and level are present.

\(\bullet\) Double exponential smoothing of the raw series

This method accounts for trend but not seasonality and therefore is not suitable for forecasting this raw series. If this series has M seasons then we could apply lag-M differences which would remove the seasonality and would then make double exponential smoothing a viable method.

\(\bullet\) Holt-Winter’s exponential smoothing of the raw series

This method accounts for both trend and seasonality and is thus a suitable method.

  1. The plot of the series with Holt-Winter’s exponential smoothing method forecasts for the validation period is plotted below:
# validation period is the last 4 quarters (last year)
validLength <- 4
# training period length
trainLength <- length(salesTS) - validLength
# partitioned data
salesTrain <- window(salesTS,end=c(1,trainLength))
salesValid <- window(salesTS,start=c(1,trainLength+1))
sesSales <- ets(salesTrain,restrict=FALSE,model="ZZZ",alpha=0.2,beta=0.15,gamma=0.05)
# forecasts for the validation period
sesSalesPredictions <- forecast(sesSales,h=validLength,level=0)
# set up the plot
yrange = range(salesTS)
plot(c(1,7),yrange,type="n",xlab="Years",ylab="Sales (thousands)",bty="l",xaxt="n",yaxt="n")
axis(1, at=seq(1,7,1), labels=format(seq(1,7,1)))
axis(2, at=seq(50000,110000,10000), labels=format(seq(50,110,10)), las=2)
abline(v=6)
arrows(1, 100000, 6, 100000, code=3, length=0.1)
text(3.5, 102000, "Training")
abline(v=7)
arrows(6, 100000, 7, 100000, code=3, length=0.1)
text(6.5, 102000, "Validation")
# time series for training set
lines(salesTS,bty="l")
# add forecasts for validation period
lines(sesSalesPredictions$mean,lwd=2,col="red",lty=2)
legend(1,100000, c("Actual", "Forecast"),lwd=2, lty=c(1,2), col=c("black", "red"), bty="n")

  1. The forecast predictions for the validation period are:
sesSalesPredictions
##      Point Forecast      Lo 0      Hi 0
## 6 Q1       62115.16  62115.16  62115.16
## 6 Q2       65371.60  65371.60  65371.60
## 6 Q3       77076.87  77076.87  77076.87
## 6 Q4      102937.73 102937.73 102937.73
  1. The separate MAPE values for the forecast of quarters 21 and 22 were determined to be about 2.163% and 0.727%, respectively. Combining quarters 21 and 22, the MAPE value was about 1.44%.
# MAPE value for the forecast of quarter 21
accuracy.measures.Q21 <- accuracy(sesSalesPredictions$mean[1],salesValid[1])
accuracy.measures.Q21
##                 ME     RMSE      MAE       MPE     MAPE
## Test set -1315.156 1315.156 1315.156 -2.163086 2.163086
# MAPE value for the forecast of quarter 22
accuracy.measures.Q22 <- accuracy(sesSalesPredictions$mean[2],salesValid[2])
accuracy.measures.Q22
##                 ME     RMSE      MAE        MPE      MAPE
## Test set -471.6041 471.6041 471.6041 -0.7266627 0.7266627
# combined MAPE value for the forecast of quarters 21 and 22
accuracy.measures.Q21.Q22 <- accuracy(sesSalesPredictions$mean[1:2],salesValid[1:2])
accuracy.measures.Q21.Q22
##                 ME    RMSE      MAE       MPE     MAPE
## Test set -893.3803 987.939 893.3803 -1.444875 1.444875
  1. Plots of the fit and residuals are reproduced below.
plot(c(1,6),yrange,type="n",xlab="Years",ylab="Sales (thousands)",bty="l",xaxt="n",yaxt="n")
lines(salesTrain,bty="l")
axis(1, at=seq(1,6,1), labels=format(seq(1,6,1)))
axis(2, at=seq(50000,110000,10000), labels=format(seq(50,110,10)), las=2)
lines(sesSales$fitted,col="red",lty=2,lwd=2)
legend(1,100000, c("Actual", "Forecast"),lwd=2, lty=c(1,2), col=c("black", "red"), bty="n")
title("Exp Smoothing: Actual Vs Forecast (Training Data)")

# residual plot
# residuals
plot(sesSalesPredictions$residuals,xlab="Years",ylab="Forecast Error",bty="n",xaxt="n",yaxt="n")
axis(1,pos=c(0,0), at=seq(1,6,1), labels=format(seq(1,6,1)))
axis(2, at=seq(-0.04,0.06,0.02), labels=format(seq(-0.04,0.06,0.02)), las=2)
abline(h = 0)
title("Exp Smoothing Forecasting Errors (Training Data)")

As we learned previously, Holt-Winter’s method is suitable for series that contain both trend and seasonality which the Department Store Sales data does indeed have. We’ve seen a MAPE value of about 1.44% for quarters 21 and 22, which is a very low mean percentage deviation. We also observe relatively small forecast errors from the above plot. All in all, I would say this model is very suitable for forecasting quarters 21 and 22.

  1. We’ll first remove trend then seasonality and then the opposite order. To remove trend we apply a lag-1 difference and for seasonality a lag-4 difference (quarterly seasons, M=4).
par(mfrow=c(2,2))
# plot original series
plot(salesTS,ylab="Sales", xlab="Years",bty="l",main="Department Store Sales")
# lag-1 difference (de-trended)
lag1 <- diff(salesTS,lag=1)
# lag-12 difference (deseasonalized)
lag4 <- diff(salesTS,lag=4)

# double difference (trend then seasonality)
# plot lag-1 difference (de-trended)
plot(lag1,ylab="Lag-1",xlab="Years",bty="l",main="Lag-1 Difference")
# plot lag-4 difference (deseasonalized)
plot(lag4,ylab="Lag-4",xlab="Years",bty="l",main="Lag-4 Difference")
# double difference (de-trended then deseasonalized)
lag1ThenLag4 <- diff(diff(salesTS, lag=1),lag=4)
# plot the double difference
plot(lag1ThenLag4, ylab="Lag-1, then Lag-4", xlab="Years", bty="l", main="Twice-Differenced (Lag-1, Lag-4)")

Now we reverse the order by first removing seasonality then trend.

par(mfrow=c(2,2))
# plot original series
plot(salesTS,ylab="Sales", xlab="Years",bty="l",main="Department Store Sales")
# plot lag-4 difference (deseasonalized)
plot(lag4,ylab="Lag-4",xlab="Years",bty="l",main="Lag-4 Difference")
# plot lag-1 difference (de-trended)
plot(lag1,ylab="Lag-1",xlab="Years",bty="l",main="Lag-1 Difference")
# double difference (deseasonalized then de-trended)
lag4ThenLag1 <- diff(diff(salesTS, lag=4),lag=1)
# plot the double difference
plot(lag4ThenLag1, ylab="Lag-4, then Lag-1", xlab="Years", bty="l", main="Twice-Differenced (Lag-4, Lag-1)")

The order of differencing does not matter as evidenced from the two double-differenced plots (they are identical).

  1. Now we use the average of the double-differenced series obtained from (d) to forecast quarters 21-22. The point forecasts are
pointForecasts <- meanf(diff(diff(salesTrain, lag=4), lag=1), h=4)
pointForecasts
##      Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 6 Q1          569.2 -2116.935 3255.335 -3714.114 4852.514
## 6 Q2          569.2 -2116.935 3255.335 -3714.114 4852.514
## 6 Q3          569.2 -2116.935 3255.335 -3714.114 4852.514
## 6 Q4          569.2 -2116.935 3255.335 -3714.114 4852.514

Now converting back for the original time series…

realForecasts <- vector()

for (i in 1:validLength) {
  if(i == 1) {
    realForecasts[i] <- pointForecasts$mean[i] + salesTrain[(trainLength+i)-validLength] + (salesTrain[trainLength] - salesTrain[trainLength - validLength])
  } else {
    realForecasts[i] <- pointForecasts$mean[i] + salesTrain[(trainLength+i)-validLength] + (realForecasts[i-1] - salesTrain[trainLength+i-1-validLength])
  }
}
#plot
par(mfrow = c(1,1))
plot(realForecasts, type="l", bty="l")

realForecasts
## [1]  63982.2  68177.4  80201.6 101467.8

We are interested in the forecasts for quarters 21-22, which are 63982.2 and 68177.4, respectively.

  1. Let’s first plot the double-differenced forecasts against the original series.
plot(c(1,7),yrange,type="n",xlab="Years",ylab="Sales (thousands)",bty="l",xaxt="n",yaxt="n")
axis(1, at=seq(1,7,1), labels=format(seq(1,7,1)))
axis(2, at=seq(50000,110000,10000), labels=format(seq(50,110,10)), las=2)
abline(v=6)
arrows(1, 100000, 6, 100000, code=3, length=0.1)
text(3.5, 102000, "Training")
abline(v=7)
arrows(6, 100000, 7, 100000, code=3, length=0.1)
text(6.5, 102000, "Validation")
# time series for training set
lines(salesTS,bty="l")
# add forecasts for validation period
lines(seq(6, 7-1/validLength, 1/validLength), realForecasts, col="red", lwd=2, lty=2)
legend(1,100000, c("Actual", "Double-Differenced Forecasts"),lwd=2, lty=c(1,2), col=c("black", "red"), bty="n")

We can also look at some accuracy measures.

accuracy(realForecasts, salesValid)
##                ME    RMSE     MAE       MPE     MAPE        ACF1 Theil's U
## Test set -1948.75 2942.41 2883.35 -3.159236 4.063656 -0.07650809 0.1948234

Notice the MAPE of about 4.06% for the double differenced series is greater than the MAPE we obtained from Holt-Winter’s method. I would use Holt-Winter’s method due to its lower MAPE value for the validation period in addition to its simple execution in R. We also have more flexibility using Holt-Winter since we can tune the 3 parameters (alpha, beta, gamma) to our desire.

  1. Comparing Holt-Winter’s method to a seasonal naive method is an even simpler approach. Let’s plot the seasonal naive forecast against the original series.
salesSnaive <- snaive(salesTrain,h=validLength)
salesSnaive$mean
##    Qtr1  Qtr2  Qtr3  Qtr4
## 6 56405 60031 71486 92183
plot(c(1,7),yrange,type="n",xlab="Years",ylab="Sales (thousands)",bty="l",xaxt="n",yaxt="n")
axis(1, at=seq(1,7,1), labels=format(seq(1,7,1)))
axis(2, at=seq(50000,110000,10000), labels=format(seq(50,110,10)), las=2)
abline(v=6)
arrows(1, 100000, 6, 100000, code=3, length=0.1)
text(3.5, 102000, "Training")
abline(v=7)
arrows(6, 100000, 7, 100000, code=3, length=0.1)
text(6.5, 102000, "Validation")
# time series for training set
lines(salesTS,bty="l")
# add forecasts for validation period
lines(salesSnaive$mean, col="blue", lwd=3, lty=3)
lines(sesSalesPredictions$mean,lwd=2,col="red",lty=2)
legend(1,100000, c("Actual","Holt-Winter", "Seasonal Naive Forecast"),lwd=2, lty=c(1,2,3), col=c("black", "red","blue"), bty="n")

Note the seasonal naive forecast is always underpredicting. Let’s examine some accuracy measures:

accuracy(salesSnaive,salesValid)
##                   ME     RMSE     MAE      MPE     MAPE     MASE
## Training set 2925.25 3878.784 3116.50 4.344358 4.737739 1.000000
## Test set     6482.25 7032.176 6482.25 8.170540 8.170540 2.079978
##                    ACF1 Theil's U
## Training set 0.54855812        NA
## Test set     0.01334402 0.4705349

We obtain a MAPE of about 4.74%, which is actually pretty good considering the simplicity of the method (only about 0.68% higher percentage deviation than the double-difference method). The decision to choose a method really depends on the forecasting goal, and if the accuracy measures obtained from the seasonal naive forecast are adequate, it may be preferred due to its simplicity.

Q.8 \(\textit{Forecasting Australian Sales}\)

  1. It is easily observed that all the series have trend and/or seasonality (most have both), therefore I would use Holt-Winter’s smoothing method to account for both trend and seasonality.

  2. Forecasts for the next two months (January 1994, February 1994) are

# read in the data
aus.wine <- read.csv("AustralianWines.csv")
fortwineTS <- ts(aus.wine$Fortified,start=c(1980,1),frequency=12)
validLength <- 12
trainLength <- length(fortwineTS) - validLength
fortwineTrain <- window(fortwineTS,end=c(1980,trainLength))
fortwineValid <- window(fortwineTS,start=c(1980,trainLength+1),end=c(1980,trainLength+validLength))
sesFortwine <- ets(fortwineTrain,model="ZZM",restrict=FALSE,alpha = 0.2,beta = 0.01,gamma = 0.01,allow.multiplicative.trend = TRUE)
fortwinePred <- forecast(sesFortwine,h=validLength,level=0)
fortwinePred$mean[1:2]
## [1] 2034.799 1899.206
plot(fortwinePred,bty="l")

Looking at some accuracy measures for the validation period…

accuracy(fortwinePred$mean[1:2],fortwineValid[1:2])
##               ME     RMSE     MAE      MPE     MAPE
## Test set 248.997 281.9171 248.997 10.78575 10.78575

The MAPE is about 10.786 percent for these next two forecasted months.

  1. The residual plot is produced below.
plot(fortwinePred$residuals,xlab="Year",ylab="Forecast Error",bty="n",xaxt="n")
axis(1, at=seq(1980,1995,1), labels=format(seq(1980,1995,1)))
abline(h = 0)
title("Exp Smoothing Forecasting Errors (Training Data)")

\(\bullet\) Decembers (month 12) are not captured well by the model.

As observed from the residual plot, there appears to be both positive and negative peak errors occuring several times in the December months. Thus I would say this statement is reasonable.

\(\bullet\) There is a strong correlation between sales on the same calendar month

Seasonality seems to be present in both the raw series and the series produced from Holt-Winter’s method in addition to some apparent seasonality in the residual plot. I would say this is a reasonable statement to make.

\(\bullet\) The model does not capture the seasonality well.

This is reasonable. Seasonality still seems to be present even after applying Holt-Winter’s method.

\(\bullet\) We should first deseasonalize the data and then apply Holt-Winter’s exponential smoothing.

This is unnecessary since Holt-Winter’s method captures the seasonality. We could however use Holt’s linear trend model.

  1. I used the parameter values of alpha=0.2, beta=0.01, and gamma=0.01, so various tuning of these (specifically gamma) may be able to capture seasonality better. We could also try additive seasonality instead of multiplicative.