#############################################################################################################
#
# Example 3.8 (pages 106-7) – Figure 3.5 (page 107)
#
# Bootstrap filter applied to the first order dynamic linear model of Example 2.9 (pages 63-4).#
#
# For t=1,...,T
#
#   y(t) =              x(t)   + v(t)     v(t) ~ N(0,V)    
#   x(t) = alpha + beta*x(t-1) + w(t)     w(t) ~ N(0,W)
#
#   x(0) ~ N(m0,C0)
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
quant025 = function(x){quantile(x,0.025)}
quant975 = function(x){quantile(x,0.975)}
set.seed(9784152)

# first order DLM
# ---------------
T     = 50
alpha = -0.01328655
beta  = 0.98
sV    = 0.25
sW    = 0.1
W     = sW^2
V     = sV^2
y     = rep(0,2*T)
x     = rep(0,2*T)
for (t in 2:(2*T)){
  x[t] = rnorm(1,alpha+beta*x[t-1],sW)
  y[t] = rnorm(1,x[t],sV)
}
x = x[(T+1):(2*T)]
y = y[(T+1):(2*T)]

y[27] = 1.0 # Forcing the presence of an outlier!

# Exact calculations
# ------------------
m = rep(0,T)
C = rep(0,T)
m0   = 0
C0   = 100
R    = C0+W
Q    = R+V
A    = R/Q
m[1] = m0+A*(y[1]-m0)
C[1] = R-A^2*Q
for (t in 2:T){
  R    = C[t-1]+W
  Q    = R+V
  A    = R/Q
  m[t] = m[t-1]+A*(y[t]-m[t-1])
  C[t] = R-A^2*Q
}

# Sequential Monte Carlo Algorithm - Bootstrap Filter 
# ---------------------------------------------------
j = 0
letter = c("(a)","(b)","(c)","(d)")
par(mfrow=c(2,2))
for (M in c(100,500,1000,10000)){
  j = j + 1 
  h = rnorm(M,m0,sqrt(C0))
  hs = NULL
  for (t in 1:T){
    h1 = rnorm(M,alpha+beta*h,sW)
    w  = dnorm(y[t],h1,sV)
    w  = w/sum(w)
    h  = sample(h1,size=M,replace=T,prob=w)
    hs = cbind(hs,h)
  }
  h500 = apply(hs,2,median)
  h975 = apply(hs,2,quant975)
  h025 = apply(hs,2,quant025)
  h500 = apply(hs,2,mean)
  sd   = sqrt(apply(hs,2,var))
  h975 = h500+2*sd
  h025 = h500-2*sd
  
  # Figure 3.5 
  # ----------
  plot(22:32,rep(-1.1,11),ylim=c(-1.2,-0.27),col=0,xlab="time",ylab="",axes=F,main=letter[j])
  axis(2)
  axis(1,at=22:32)
  for (i in 22:32){
    segments(i,m[i]+qnorm(0.025)*sqrt(C[i]),i,m[i]+qnorm(0.975)*sqrt(C[i]))
    points(i,m[i],pch=15)
    segments(i+0.3,quantile(hs[,i],0.025),i+0.3,quantile(hs[,i],0.975),lty=3)
    points(i+0.3,mean(hs[,i],0.5),pch=16)
  }
}