This a multipart series aiming to compare and contrast the various Holt Winters implementations in R. We intend to focus more on the practical and applied aspects of the implementations to get a better grip over the behaviour of models and predictions. So to begin with lets look at the *‘HoltWinters()’* function in stats package that comes built-in with R and the *‘hw()’* function in forecast package by Dr. Rob J Hyndman.

To install the forecast package, use the following command:

install.packages("forecast")

For this analysis, I have chosen to use a time series data available in rdatamarket. The data used is Monthly total number of pigs slaughtered in Victoria. Jan 1980 – August 1995. Since this is available in rdatamarket, we can directly load it in R using the *rdatamarket* package.

So lets load the package and download the data:

This is how the time series appears. On visual inspection, there clearly is a cyclicity and some trend elements involved. As you can see from the pattern, one can conclude that it follows an *additive pattern on both seasonal and trend components.*

Our approach here is to compare primarily on the accuracy of the forecasting performance of these functions rather than the closeness of the fit of training data. And the accuracy measure I have chosen for this analysis is Mean Absolute Percentage Deviation (MAPE).

Our aim here is not to conclude that one function predicts better than the other, but to come to a experiential understanding over how the functions behave depending on the nature of the timeseries. I also would like to point out that the selection of this data for this analysis is purely arbitrary and I encourage you to plug in this code to other time series data you may come across.

We have loaded the data and downloaded the packages, So lets begin..

About 30% of the available data is help for testing and we use the initial 70% of the data to build our Holt Winters models. Both the function HoltWinters() and hw() try to optimize the alpha, beta and gamma values by minimizing the residuals but both does not give the same values. Hence the differing fits and predictions.

The output of the following code stores the forecasting charts in your working directory. Each chart shows the actual, fit and forecast with 95% confidence drawn in orange lines.

The accuracy measures that gets printed on the console are computed for the training data and the MAPE is for the predictions made by the respective models.

Just by looking at the fitted curve of the training data, one can argue that HoltWinters() performed equally well if not better than hw(). But when you look at forecasted values, clearly in this case, the predictions of hw() performed better than HoltWinters().

Lets look at the accuracy measures..

*Interestingly, the MAPE of fitted values was slightly better for HoltWinters() over hw(), but the MAPE when calculated for forecast vs actual, the hw() function of ‘forecast’ package clearly outperforms. *

library(rdatamarket) library(forecast) pigs <- dmseries("http://data.is/H63F9L") data <- pigs data.ts <- as.ts(data) #plot(data.ts, type = 'b', col = "blue", main = colnames(data.ts)) #Training and Test Split split <- ceiling(0.7 * length(data)) dataTrain <- ts(data[1:split], frequency = 12, start = c(1980, 1)) dataTest <- ts(data[c((split+1) : nrow(data))], frequency = 12, start = c(1991, 1)) actual <- unclass(dataTrain) actualFull <- unclass(data) #metd values => HoltWinters, hw forecastGen <- function(metd){ if(metd == "HoltWinters"){ method <- "HoltWinters {stats}" mod <- HoltWinters(dataTrain, optim.start = c(alpha = 0.99, beta = 0.001, gamma =0.001)) model <- forecast.HoltWinters(mod, h=(length(data)-length(dataTrain))) actualFull <- actualFull[13:length(actualFull)] #drop the first 12 obs actual <- actual[13:length(actual)] #drop the first 12 obs i = 4 } if(metd == "hw"){ method <- "hw {forecast}" model <- hw(dataTrain, initial = "optimal",h=(length(data)-length(dataTrain))) i = 2 } filename <- paste(method, ".png", sep = "") fit <- as.data.frame(model$fitted)[,1] # Fitted Values fitV <- unclass(fit) fc <- unlist(model[i]) names(fc) <- NULL fitForecast <- (c(fit, fc)) comp <- cbind(dataTest, fc) compFull <- cbind(actualFull, fitForecast) cat("Accuracy of Model:",metd,"\n") print(accuracy(model)) mape <- round(mean(abs((comp[,2]-comp[,1])/comp[,1]))*100, 2) cat("\n",paste("MAPE of Predictions =" ,mape), "\n\n") upper <- c(actual,model$upper[,1]) lower <- c(actual,model$lower[,1]) #PLOT and SAVE png(file=filename,width=900,height=550) plot(actualFull, type = 'l', col = 'black', lwd = 2.5, main = method, ylab = colnames(data.ts), xlab = "Time Intervals", ylim = c(min(lower),max(upper))) grid() lines(fitForecast, type = 'l', col = 'red', lwd = 2) lines(fitV, type = 'l', col = 'blue', lwd = 3) lines(upper, col = "orange", lwd = 0.5) lines(lower, col = "orange", lwd = 0.5) legend("bottomright", inset=.01, c("Actual","Fitted", "Predicted"), fill=c(1,4,2), horiz=TRUE, cex=0.75 ) invisible(dev.off()) } #call the functions on HoltWinters() and hw() forecastGen(metd = "HoltWinters") forecastGen(metd = "hw")

If you are more interested in theoretical views you might find some related discussions here. If you are a beginner to R Programming language, this should help you to get upto speed.

Author: Selva Prabhakaran

Selva Prabhakaran Sanjeevi Julian