#####################################################################################
#
# Example 1: Local level model
#
# Lopes and Tsay (2011)
# Particle Filters and Bayesian Inference in Financial Econometrics.
# Journal of Forecasting.
# 
#####################################################################################
#
# Comparison of four particle filters for the local level model.
#
# For t = 1,...,n
#
#    y[t] ~ N(x[t],sig2)
#    x[t] ~ N(alpha*x[t-1],tau2)
# 
# with x[0] ~ N(m0,C0).
#
# Known quantities: (alpha,sig2,tau2,m0,C0).
#
# Particle filters: 
#    Bootstrap filter (Gordon et al, 1993)
#    Auxiliary particle filter (Pitt and Shephard, 1999)
#    Optimal bootstrap filter
#    Optimal auxiliary particle filter
#
######################################################################################
#
# Author : 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
#
#####################################################################################
rm(list=ls())

simDLM = function(n,alpha,sig,tau,x0,m0,C0){
	y  = rep(0,n+1)
	x  = rep(x0,n+1)
  for (t in 2:(n+1)){
    x[t] = rnorm(1,alpha*x[t-1],tau)
    y[t] = rnorm(1,x[t],sig)
  }
  x = x[2:(n+1)]
  y = y[2:(n+1)]
  return(list(x=x,y=y))
}

ffDLM = function(y,alpha,sig2,tau2,m0,C0){
  n  = length(y)
  m  = m0
  C  = C0
  qs = NULL
  for (t in 1:n){
    a = alpha*m
    R = (alpha^2)*C+tau2
    Q = R+sig2
    A = R/Q
    m = (1-A)*a+A*y[t]
    C = A*sig2
    quants = qnorm(c(0.25,0.5,0.75),m,sqrt(C))
    qs = rbind(qs,quants)
  }
  return(qs)
}

PF = function(y,alpha,sig2,tau2,m0,C0,N,index){
  sig = sqrt(sig2)
  tau = sqrt(tau2)
  if (index==1){output=BF(y,alpha,sig,tau,m0,C0,N)}
  if (index==2){output=OBF(y,alpha,sig2,tau2,m0,C0,N)}
  	if (index==3){output=APF(y,alpha,sig,tau,m0,C0,N)}
  	if (index==4){output=OAPF(y,alpha,sig2,tau2,m0,C0,N)}
  	return(output)
}
  	
BF = function(y,alpha,sig,tau,m0,C0,N){
  n  = length(y)
  xs = rnorm(N,m0,sqrt(C0))
  qs = NULL
  for (t in 1:n){
    x1   = rnorm(N,alpha*xs,tau)
    w    = dnorm(y[t],x1,sig)
    xs   = sample(x1,size=N,replace=TRUE,prob=w)
    quants = quantile(xs,c(0.25,0.5,0.75))
    qs   = rbind(qs,quants)
  }
  return(qs)
}

OBF = function(y,alpha,sig2,tau2,m0,C0,N){
  n   = length(y)
  xs  = rnorm(N,m0,sqrt(C0))
  var1 = 1/(1/sig2+1/tau2)
  sd  = sqrt(var1)
  sd1 = sqrt(sig2+tau2)
  qs = NULL
  for (t in 1:n){
  	  mean = var1*(alpha*xs/tau2+y[t]/sig2)
    x1 = rnorm(N,mean,sd)
    w  = dnorm(y[t],alpha*xs,sd1)
    xs = sample(x1,size=N,replace=TRUE,prob=w)
    quants = quantile(xs,c(0.25,0.5,0.75))
    qs     = rbind(qs,quants)
  }
  return(qs)
}

APF = function(y,alpha,sig,tau,m0,C0,N){
  n  = length(y)
  xs = rnorm(N,m0,sqrt(C0))
  qs = NULL
  for (t in 1:n){
    w    = dnorm(y[t],alpha*xs,sig)
    k    = sample(1:N,size=N,replace=TRUE,prob=w)
    x1   = rnorm(N,alpha*xs[k],tau)
    w1   = dnorm(y[t],x1,sig)/w[k]
    xs   = sample(x1,size=N,replace=TRUE,prob=w1)
    quants = quantile(xs,c(0.25,0.5,0.75))
    qs     = rbind(qs,quants)
  }
  return(qs)
}

OAPF = function(y,alpha,sig2,tau2,m0,C0,N){
  n   = length(y)
  xs  = rnorm(N,m0,sqrt(C0))
  var1 = 1/(1/sig2+1/tau2)
  sd  = sqrt(var1)
  sd1 = sqrt(sig2+tau2)
  qs  = NULL
  for (t in 1:n){
    w      = dnorm(y[t],alpha*xs,sd1)
    xso    = xs
    k      = sample(1:N,size=N,replace=TRUE,prob=w)
    xs     = xs[k]    
  	  mean  = var1*(alpha*xs/tau2+y[t]/sig2)
    xs     = rnorm(N,mean,sd)
    quants = quantile(xs,c(0.25,0.5,0.75))
    qs     = rbind(qs,quants)
  }
  return(qs)
}

# -------------------------------------------------------------
# Monte Carlo study
# -------------------------------------------------------------
n     = 100
M     = 20
rep   = 20
Ns    = c(50,100)
alpha = 0.5
sig2  = 1.0
x0    = 25
m0    = x0
C0    = 10.0
tau2s = c(0.05,0.5,1.0,2.0)
taus  = sqrt(tau2s)
ntau  = length(tau2s)
sig   = sqrt(sig2)
sC0   = sqrt(C0)
nN    = length(Ns)
yy    = array(0,c(ntau,M,n))
xx    = array(0,c(ntau,M,n))
qtrue = array(0,c(ntau,M,n,3))
set.seed(123456)
for (i in 1:ntau){
  tau2 = tau2s[i]
  tau  = sqrt(tau2)
  for (j in 1:M){
    sim      = simDLM(n,alpha,sig,tau,x0,m0,C0)
    yy[i,j,] = sim$y
	  xx[i,j,] = sim$x
    qtrue[i,j,,] = ffDLM(yy[i,j,],alpha,sig2,tau2,m0,C0)
  }
}

filters = c("BF","OBF","APF","OAPF")
nfilter = length(filters)
pf  = array(0,c(4,ntau,nN,M,rep,n,3))
for (f in 1:nfilter){
  set.seed(654321)
  for (k in 1:ntau){
    for (i in 1:nN){
      print(paste(filters[f]," - N=",Ns[i]," - tau=",taus[k],sep=""))
      for (j in 1:M)
        for (l in 1:rep)
          pf[f,k,i,j,l,,] = PF(yy[k,j,],alpha,sig2,tau2s[k],m0,C0,Ns[i],f)
      }
    }
  }

mse = array(0,c(4,ntau,nN,3))
for (i in 1:4)
  for (j in 1:ntau)
    for (k in 1:nN)
      for (l in 1:3)
        for (m in 1:rep)
          mse[i,j,k,l] = mse[i,j,k,l]+sum((pf[i,j,k,,m,,l]-qtrue[j,,,l])^2)/(n*rep*Ns[k])
      
rmse = array(0,c(3,ntau,nN,3))
for (i in 1:3)
  rmse[i,,,] = mse[i+1,,,]/mse[1,,,]
  
pdf(file="example1-1.pdf",width=15,height=10)
par(mfrow=c(1,2))  
for (j in 1:nN){
  plot(taus,rmse[1,,j,1],xlab="tau",ylab="RMSE",pch=16,type="b",ylim=range(rmse),axes=F)
  axis(2);axis(1,at=taus,lab=round(taus,2));box()
  for (i in 1:3)
    for (l in 1:3)
      lines(taus,rmse[i,,j,l],col=i,pch=16,type="b")
  abline(h=1,lty=3)
  title(paste("N=",Ns[j],sep=""))
  if (j==1){
    legend(taus[2],max(rmse),legend=filters[2:4],col=1:3,lwd=rep(2,3),lty=rep(1,3),bty="n",cex=1.5)
  }
}
dev.off()