####################################################################
# 
# LIU AND WEST FILTER
#
#   y(t) = x(t)/(1+x(t)^2) + v(t)
#   x(t) = x(t-1) + w(t)
#
#   v(t) ~ N(0,1)
#   w(t) ~ N(0,0.5)
#   x(0) ~ N(1,10)
#
# Reference: Liu and West (2001)
#
####################################################################
#
# 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@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
####################################################################
q025 = function(x){quantile(x,0.025)}
q975 = function(x){quantile(x,0.975)}

set.seed(1235)
n    = 500
sig2 = 0.5
tau2 = 0.5
x    = rep(0,n)
y    = rep(0,n)
sig  = sqrt(sig2)
x[1] = 2
y[1] = rnorm(1,x[1]/(1+x[1]^2),sig)
for (t in 2:n){
  x[t] = x[t-1] + rnorm(1,0,sqrt(tau2))
  y[t]     = rnorm(1,x[t]/(1+x[t]^2),sig)
}
par(mfrow=c(1,2))
ts.plot(x)
ts.plot(y)
xtrue = x
tau2true = tau2

# Prior hyperparameters
a0 = 10
b0 = (a0-1)*tau2true


sig2 = 0.25
sig = sqrt(sig2)

# Particle filter setup
set.seed(923579)
N    = 1000
xs   = matrix(0,n,N)
ts   = matrix(0,n,N)
x    = rnorm(N,1,sqrt(10))
tau2 = 1/rgamma(N,a0,b0)
a    = 0.99
pred = rep(0,n)
for (t in 1:n){
  # computing sequential predictive: p(y(t)|y(1)....y(t-1))
  xt    = x + rnorm(N,0,sqrt(tau2)) 
  pred[t] = mean(dnorm(y[t],xt,sig))
   	
  theta = log(tau2)
  thbar = mean(theta)
  V     = var(theta)
  m     = a*theta + (1-a)*thbar
  gx    = x
  my    = gx/(1+gx^2)
  w     = dnorm(y[t],my,sig)
  k     = sample(1:N,size=N,replace=TRUE,prob=w)
  theta = rnorm(N,m[k],sqrt((1-a^2)*V))
  tau2  = exp(theta)
  x1    = rnorm(N,x[k],sqrt(tau2))
  my1   = x1/(1+x1^2)
  w1    = dnorm(y[t],my1,sig,log=TRUE)-dnorm(y[t],my[k],sig,log=TRUE)
  w1    = exp(w1-max(w1))
  k     = sample(1:N,size=N,replace=TRUE,prob=w1)
  x     = x1[k]
  tau2  = tau2[k]
  xs[t,] = x
  ts[t,] = tau2
  hist(tau2,xlab="",ylab="",main=t)
}

mx = apply(xs,1,median)
lx = apply(xs,1,q025)
ux = apply(xs,1,q975)
Lx = min(lx)
Ux = max(ux)
mt = apply(ts,1,median)
lt = apply(ts,1,q025)
ut = apply(ts,1,q975)
Lt = min(lt)
Ut = max(ut)

par(mfrow=c(1,2))
plot(mx,type="l",ylim=c(Lx,Ux),xlab="Time")
title(paste("N=",N,sep=""))
lines(lx)
lines(ux)
lines(xtrue,col=2)
plot(mt,type="l",ylim=c(Lt,Ut),xlab="Time")
title(paste("N=",N,sep=""))
lines(lt)
lines(ut)
abline(h=tau2true,col=2)

par(mfrow=c(1,1))   
ts.plot(pred,main="PREDICTIVE")

sum(log(pred))

model 1: -4053.285  (sig2=0.5)
model 2: -2522.560  (sig2=1.0)
model 3: -1170.387  (sig2=0.25)