dble.exp <- function(X, alpha, beta) {
	n <- length(X)
	L <- rep(NA, n)
	B <- rep(NA, n)
	Fc <- rep(NA, n+1)
	L[1] <- X[1]
	B[1] <- X[2] - X[1]
	Fc[1:2] <- X[1:2]
	for (t in 2:n) {
		L[t] <- alpha*X[t] + (1 - alpha)*(L[t-1] + B[t-1])
		B[t] <- beta*(L[t] - L[t-1]) + (1 - beta)*B[t-1]
		Fc[t+1] <- L[t] + B[t]
	}
	return(Fc)
}

# setwd("D:/")
x <- read.csv("demand.csv")
x <- x$Observed
plot(c(1,25), c(0.9*min(x), 1.1*max(x)), type='n', xlab='t', ylab='X(t)')
lines(1:24, x)
points(1:24, x)

F.exp <- dble.exp(x, .3, .7)
mse.exp <- mean((F.exp[13:24] - x[13:24])^2)
cat('mse for double exp smoother is', mse.exp, '\n')
cat('forecast for t = 25 is', F.exp[25], '\n')
lines(13:25, F.exp[13:25], col='blue')
points(13:25, F.exp[13:25], col='blue')