ema-predictions.r

# Exponential-moving average for visitors.csv dataset.

# load dataset (assumes in same directory)
data1<-read.csv("visitors.csv");

# reverse order of dataset to ascending
data<-data1[order(data1$Date, decreasing = FALSE), ];

ema<-function(x,p,intervalratio)
	# where 
	# x = actual values
	# p = number of periods
	# intervalratio = user defined variable applied to residuals to calculate forecast range
{
	# create a dataset for collecting all usefull infromation
	result2<-data.frame(date=data[,1],act.val=x,cal.val=x,residual=rep(NA,nrow(data)),min.val=x,max.val=x,if.outlier=as.logical(rep(0,nrow(data))));

	# define alpha
	a<-2/(p+1)	# proof by Steven Nahmias

	# loop to calculate ema values and detect outliers
	for (i in 2:nrow(data)) # start with second row
	{
		result2$cal.val[i]<-a*result2$act.val[i]+(1-a)*result2$cal.val[i-1];
		result2$residual[i]<-abs(result2$act.val[i]-result2$cal.val[i]);
		residual<-mean(result2$residual,na.rm=T);
		result2$min.val[i]<-result2$cal.val[i]-intervalratio*residual;
		result2$max.val[i]<-result2$cal.val[i]+intervalratio*residual;
		if(result2$act.val[i]<result2$min.val[i]||result2$act.val[i]>result2$max.val[i]){result2$if.outlier[i]<-TRUE}
	}
	result2$out.val<-as.numeric(ifelse(result2$if.outlier=="TRUE",data[,2],""));
	return(result2);
}

# run ema function
result2<-ema(data[,2],15,5); 
	# in this example, we run the ema function using the values from 
	# second row of dataset 'data' (Total_Number_of_Visitors) using 
	# a 15-day exponential moving average and an intervalratio of 5. 
	# a 15-day period was used to 'weight' the recent values more
	# heavily and a 5 intervalratio was found by trial and error
	# (the larger the interval ratio, the wider the forecast range).
	# 'p' and 'intervalratio' can be modified as desired.

# plot ema
plot(result2$act.val,type="l",col="grey",lwd=2);
points(result2$cal.val,type="l",col="blue",lwd=2);
points(result2$max.val,type="l",col="green",lwd=2);
points(result2$min.val,type="l",col="red",lwd=2);
points(result2$out.val,col="orange",lwd=3);

holt-winters-predictions.r

# Holt-Winters computation for visitors.csv dataset.

# load dataset (assumes in same directory)
data1<-read.csv("visitors.csv");

#reverse the order of dataset to ascending
data<-data1[order(data1$Date, decreasing = FALSE), ];

# define outlier function
outlier<-function(n,f,c)
	# where
	# n = number of observations
	# f = time series frequency
	# c = confidence interval
{
	# create a time series with user submitted frequency, defined 
	# as the frequency with which the data is recorded
	ts<-ts(data,frequency=f); # holt-winters function needs time series data

	# create a dataset for collecting all usefull infromation
	result<-data.frame(time=as.numeric(ts[,1]),date=data[,1],act.val=as.ts(as.numeric(ts[,2])),pred.val=rep(NA,nrow(data)),min.forecast=rep(NA,nrow(data)),max.forecast=rep(NA,nrow(data)));

	# create a whether outlier logical list
	result$if.outlier<-as.logical(rep(0,nrow(data)));

	# loop for every observation after initially selected observation for model
	for (i in 1:n)
	{
		# subset on which next prediction is based
		sample<-as.ts(subset(ts[,2],ts[,1]<((f*2)+i))); # holt-winters needs two full historical periods,
								# which is two years in this example

		# compute holt-winters filtering 
		# gamma is false because no seasonal variation was detected otherwise
		hw<-HoltWinters(sample,gamma=FALSE,beta=FALSE,alpha=.1); # note: R is case-sensitive, so holt-winters
									 # function needs capital 'H' & 'W'

		# forecast next value for holt-winters
		forecast<-predict(hw,n.ahead=1,prediction.interval=TRUE,level=c);

		# find next actual value
		z<-ts[(f*2)+i,2];

		# detect if outlier
		if(z<min(forecast)||z>max(forecast)){result$if.outlier[(f*2)+i]<-"TRUE"}

		# feed data for future use in result dataset
		result$pred.val[(f*2)+i]<-forecast[1];
		result$min.forecast[(f*2)+i]<-min(forecast);
		result$max.forecast[(f*2)+i]<-max(forecast);
	}
	result$pred.val<-as.ts(result$pred.val,frequency=1200);
	result$min.forecast<-as.ts(result$min.forecast,frequency=1200);
	result$max.forecast<-as.ts(result$max.forecast,frequency=1200);
	result$if.outlier.n<-as.numeric(as.logical(result$if.outlier));
	result$out.val<-as.numeric(ifelse(result$if.outlier=="TRUE",ts[,2],""));
	return(result);
}

# run function
result<-outlier((nrow(data)-(365*2)), 365, 0.99);
	# in this example, the number of observations is equal to the 
	# total rows in dataset less the two historical periods required
	# by holt-winters. the frequency is daily (or 365) and we assumed
	# a 99% confidency interval. 'c' can be modified as desired.

# plot hw
plot(result$act.val,col="grey",xlab="TIME",ylab="VISITORS",ylim=c(0,20000),lwd=2);
par(new=T);
plot(result$min.forecast,col="red",axes=F,xlab="",ylab="",ylim=c(0,20000),lwd=2);
par(new=T);
plot(result$pred.val,col="blue",axes=F,xlab="",ylab="",ylim=c(0,20000),lwd=2);
par(new=T);
plot(result$max.forecast,col="green",axes=F,xlab="",ylab="",ylim=c(0,20000),lwd=2);
par(new=T);
points(result$out.val,col="purple",ylim=c(0,20000),lwd=3);

results.csv

visitors.csv