install.packages("bsts")
library("bsts")

data = matrix(scan("pollutiondata-NO3.txt"),342,22,byrow=T)

par(mfrow=c(2,4))
ind = sort(sample(1:22,replace=F,size=8))
for (i in ind)
  ts.plot(data[,i],main=paste("Station ",i,sep=""),ylab="")

y = data[,1]
par(mfrow=c(1,1))
ts.plot(y)

ss   = AddLocalLinearTrend(list(),y)

ss   = AddSeasonal(ss, nseasons = 52,y)

ssm  = bsts(y,state.specification = ss,niter=2000,family="student")

trend = apply(ssm$state.contributions[,1,],2,median)
seasonality = apply(ssm$state.contributions[,2,],2,median)
error = y-trend-seasonality

par(mfrow=c(2,4))
ts.plot(y)
acf(y,lag.max=100)
ts.plot(y)
lines(trend,col=2,lwd=2)
ts.plot(y)
lines(trend+seasonality,col=2,lwd=2)
ts.plot(error)
acf(error)

hist(ssm$sigma.obs,prob=TRUE,xlab="Standard deviation",main="Observational");box()
boxplot(ssm$sigma.trend.level,ssm$sigma.trend.slope,ssm$sigma.seasonal.52,
names=c("Level","Slope","Seasonality"),ylab="Standard deviation",outline=FALSE,main="Evolutional")