##########################################################################################
#
#  http://people.ischool.berkeley.edu/~hal/Papers/2013/pred-present-with-bsts.pdf
#
##########################################################################################
install.packages("bsts")
library(bsts)
data(AirPassengers)
y = log(AirPassengers)

pdf(file="bsts-airline.pdf",width=6,height=4)
plot(y)
dev.off()

############################################################
# Model with linear trend and seasonal dummies
############################################################
ss = AddLocalLinearTrend(list(), y)
ss = AddSeasonal(ss, y, nseasons = 12)
model = bsts(y, state.specification = ss, niter = 10000)

sigmas = cbind(model$sigma.obs,model$sigma.trend.level,
model$sigma.trend.slope,model$sigma.seasonal.12)

round(cbind(apply(sigmas[5001:10000,],2,mean),
apply(sigmas[5001:10000,],2,median),
apply(sigmas[5001:10000,],2,quantile,0.025),
apply(sigmas[5001:10000,],2,quantile,0.975)),4)

pdf(file="bsts-graph1.pdf",width=10,height=7)
par(mfrow=c(3,4))
for (i in 1:4)
  ts.plot(sigmas[,i],ylab="")
for (i in 1:4)
  acf(sigmas[5001:10000,i],main="")
for (i in 1:4)
  hist(sigmas[5001:10000,i],xlab="",main="",prob=TRUE)
dev.off()


components = cbind.data.frame(
  colMeans(model$state.contributions[-(1:5000),"trend",]),                               
  colMeans(model$state.contributions[-(1:5000),"seasonal.12.1",]),
  as.Date(time(y)))  
names(components) = c("Trend", "Seasonality", "Date")
components = melt(components, id="Date")
names(components) = c("Date", "Component", "Value")

pdf(file="bsts-graph2.pdf",width=10,height=7)
ggplot(data=components, aes(x=Date, y=Value)) + geom_line() + 
  theme_bw() + theme(legend.title = element_blank()) + ylab("") + xlab("") + 
  facet_grid(Component ~ ., scales="free") + guides(colour=FALSE) + 
  theme(axis.text.x=element_text(angle = -90, hjust = 0))
dev.off()


############################################################
# Model with linear trend and trigonometric seasonality
############################################################
data(AirPassengers)
y <- log(AirPassengers)
ss <- AddLocalLinearTrend(list(), y)
ss <- AddTrig(ss, y, period = 12, frequencies = 1)
model1 = bsts(y, state.specification = ss, niter = 10000)

sigmas1 = cbind(model1$sigma.obs,model1$sigma.trend.level,model1$sigma.trend.slope,model1$trig.coefficient.sd[,1],model1$trig.coefficient.sd[,2])

round(cbind(apply(sigmas1[5001:10000,],2,mean),
apply(sigmas1[5001:10000,],2,median),
apply(sigmas1[5001:10000,],2,quantile,0.025),
apply(sigmas1[5001:10000,],2,quantile,0.975)),4)

par(mfrow=c(3,5))
for (i in 1:5)
  ts.plot(sigmas1[,i],ylab="")
for (i in 1:5)
  acf(sigmas1[5001:10000,i],main="")
for (i in 1:5)
  hist(sigmas1[5001:10000,i],xlab="",main="",prob=TRUE)


components = cbind.data.frame(
  colMeans(model1$state.contributions[-(1:5000),1,]),                               
  colMeans(model1$state.contributions[-(1:5000),2,]),
  as.Date(time(y)))  
names(components) = c("Trend", "Seasonality", "Date")
components = melt(components, id="Date")
names(components) = c("Date", "Component", "Value")


pdf(file="bsts-graph3.pdf",width=10,height=7)
ggplot(data=components, aes(x=Date, y=Value)) + geom_line() + 
  theme_bw() + theme(legend.title = element_blank()) + ylab("") + xlab("") + 
  facet_grid(Component ~ ., scales="free") + guides(colour=FALSE) + 
  theme(axis.text.x=element_text(angle = -90, hjust = 0))
dev.off()