Holt-Winters forecast using ggplot2

R has great support for Holt-Winter filtering and forecasting. I sometimes use this functionality, HoltWinter & predict.HoltWinter, to forecast demand figures based on historical data. Using the HoltWinter functions in R is pretty straightforward.

Let's say our dataset looks as follows;

demand <- ts(BJsales, start = c(2000, 1), frequency = 12)
plot(demand)

Now I pass the timeseries object to HoltWinter and plot the fitted data.

hw <- HoltWinters(demand)
plot(hw)

Next, we calculate the forecast for 12 months with a confidence interval of .95 and plot the forecast together with the actual and fitted values.

forecast <- predict(hw, n.ahead = 12, prediction.interval = T, level = 0.95)
plot(hw, forecast)

As you can see, this is pretty easy to accomplish. However, as I use ggplot2 to visualize a lot of my analyses, I would like to be able to do this in ggplot2 in order to maintain a certain uniformity in terms of visualization.

Therefore, I wrote a little function which extracts some data from the HoltWinter and predict.HoltWinter objects and feeds this to ggplot2;

#HWplot.R

library(ggplot2)
library(reshape)


HWplot<-function(ts_object,  n.ahead=4,  CI=.95,  error.ribbon='green', line.size=1){
	
	hw_object<-HoltWinters(ts_object)
	
	forecast<-predict(hw_object,  n.ahead=n.ahead,  prediction.interval=T,  level=CI)
	
	
	for_values<-data.frame(time=round(time(forecast),  3),  value_forecast=as.data.frame(forecast)$fit,  dev=as.data.frame(forecast)$upr-as.data.frame(forecast)$fit)
	
	fitted_values<-data.frame(time=round(time(hw_object$fitted),  3),  value_fitted=as.data.frame(hw_object$fitted)$xhat)
	
	actual_values<-data.frame(time=round(time(hw_object$x),  3),  Actual=c(hw_object$x))
	
	
	graphset<-merge(actual_values,  fitted_values,  by='time',  all=TRUE)
	graphset<-merge(graphset,  for_values,  all=TRUE,  by='time')
	graphset[is.na(graphset$dev),  ]$dev<-0
	
	graphset$Fitted<-c(rep(NA,  NROW(graphset)-(NROW(for_values) + NROW(fitted_values))),  fitted_values$value_fitted,  for_values$value_forecast)
	
	
	graphset.melt<-melt(graphset[, c('time', 'Actual', 'Fitted')], id='time')
	
	p<-ggplot(graphset.melt,  aes(x=time,  y=value)) + geom_ribbon(data=graphset, aes(x=time, y=Fitted, ymin=Fitted-dev,  ymax=Fitted + dev),  alpha=.2,  fill=error.ribbon) + geom_line(aes(colour=variable), size=line.size) + geom_vline(x=max(actual_values$time),  lty=2) + xlab('Time') + ylab('Value') + opts(legend.position='bottom') + scale_colour_hue('')
	return(p)

}

The above script is saved in a file called HWplot.R. If I load this file from R – via source() – I can directly call the function HWplot. The HWplot can be called as follows:

HWplot(ts_object, n.ahead=4, CI=.95, error.ribbon='green',line.size=1)

HWplot takes the following arguments;

• ts_object: the timeseries data
• n.ahead: number of periods to forecast
• CI: confidence interval
• error.ribbon: colour of the error ribbon
• line.size: size of the lines

source("HWplot.R")
demand <- ts(BJsales, start = c(2000, 1), frequency = 12)
HWplot(demand, n.ahead = 12)

It's also very easy to adjust the graph after it is returned by the function;

graph <- HWplot(demand, n.ahead = 12, error.ribbon = "red")

# add a title
graph <- graph + opts(title = "An example Holt-Winters (gg)plot")

# change the x scale a little
graph <- graph + scale_x_continuous(breaks = seq(1998, 2015))

# change the y-axis title
graph <- graph + ylab("Demand ($)")

# change the colour of the lines
graph <- graph + scale_colour_brewer("Legend", palette = "Set1")

# the result:
graph

The HWplot R code: HWplot.R

The post Holt-Winters forecast using ggplot2 appeared first on FishyOperations.