Comparing Holt Winters Implementations in R – Part 1

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:

library(rdatamarket)
library(forecast)
 
pigs <- dmseries("http://data.is/H63F9L")

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.

Monthly Total Number of Pigs Slaughtered in Victoria

Monthly Total Number of Pigs Slaughtered in Victoria

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.

 

Holt Winters - hw()

Forecast Projections of Holt Winters Implementation – hw {forecast}

 

Holt Winters Projections from base 'stats' Package

Holt Winters Projections from base ‘stats’ Package

 

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..

Compare accuracy of HoltWinters() vs hw()

Compare accuracy of HoltWinters() vs hw()

 

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s