############################################################################################
# 
# AR(1) PLUS NOISE MODEL - PURE FILTER
#
#      Observation equation: y(t)|x(t),theta   ~ N(x(t),sig2)
#      Evolution equation:   x(t)|x(t-1),theta ~ N(alpha+beta*x(t-1),tau2)
#  
# where x(0) ~ N(m0,C0) and theta=(alpha,beta,sig2,tau2) is known.
#
############################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoBooth.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
############################################################################################

rm(list=ls())
source("ar1plusnoisemodel-purefilter-routines.R")

# Simulating the data
set.seed(135792468)
M      = 500
n      = 1000
sig2   = 1.0
tau2   = 0.1
alpha  = 0.05
beta   = 0.95
x0     = alpha/(1-beta)
m0     = x0
C0     = 10
theta  = c(alpha,beta,sig2,tau2)
sim    = DLM.SIM(n,theta,x0)
y      = sim$y
xtrue  = sim$x

# Plotting the simulated data
plot(y,type="b")
lines(xtrue,col=2)

# Running kalman filter
dlm   = DLM(y,theta,m0,C0)
L = dlm$m+qnorm(0.05)*sqrt(dlm$C)
U = dlm$m+qnorm(0.95)*sqrt(dlm$C)

plot(dlm$m,ylim=range(L,U),type="l",xlab="Time",ylab="State")
lines(L,lty=2)
lines(U,lty=2)
points(xtrue,col=2,pch=16)


# Running the 4 particles filters
qbf   = t(apply(BF(y,theta,m0,C0,M),1,quantile,c(0.05,0.5,0.95)))
qapf  = t(apply(APF(y,theta,m0,C0,M),1,quantile,c(0.05,0.5,0.95)))
qobf  = t(apply(OBF(y,theta,m0,C0,M),1,quantile,c(0.05,0.5,0.95)))
qoapf = t(apply(OAPF(y,theta,m0,C0,M),1,quantile,c(0.05,0.5,0.95)))
LL    = min(qbf,qapf,qobf,qoapf,L)
UU    = max(qbf,qapf,qobf,qoapf,U)

par(mfrow=c(2,2))
ts.plot(qbf,ylim=c(LL,UU),xlab="Time",ylab="",main="BF")
lines(dlm$m,col=2);lines(L,col=2);lines(U,col=2)
ts.plot(qapf,ylim=c(LL,UU),xlab="Time",ylab="",main="APF")
lines(dlm$m,col=2);lines(L,col=2);lines(U,col=2)
ts.plot(qobf,ylim=c(LL,UU),xlab="Time",ylab="",main="OBF")
lines(dlm$m,col=2);lines(L,col=2);lines(U,col=2)
ts.plot(qoapf,ylim=c(LL,UU),xlab="Time",ylab="",main="OAPF")
lines(dlm$m,col=2);lines(L,col=2);lines(U,col=2)