############################################################################################
# 
# FIRST ORDER DYNAMIC LINEAR MODEL
#
############################################################################################
#
# 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/
#
############################################################################################
rm(list=ls())
DLM.SIM = function(n,theta,x0){
  alpha = theta[1]
  beta  = theta[2]
  sig2  = theta[3]
  tau2  = theta[4]
  sig   = sqrt(sig2)
  tau   = sqrt(tau2)
  x     = rep(0,n)
  y     = rep(0,n)
  x[1]  = rnorm(1,alpha+beta*x0,tau)
  for (t in 2:n)
    x[t] = rnorm(1,alpha+beta*x[t-1],tau)
  y = rnorm(n,x,sig)  
  return(list(x=x,y=y))
}

DLM = function(y,theta,m0,C0){
  alpha = theta[1]
  beta  = theta[2]
  sig2  = theta[3]
  tau2  = theta[4]
  n     = length(y)
  m     = m0
  C     = C0
  mm    = rep(0,n)
  CC    = rep(0,n)
  for (t in 1:n){
    a = alpha+beta*m
    R = beta^2*C+tau2
    Q = R+sig2
    A = R/Q
    m = (1-A)*a+A*y[t]
    C = A*sig2
    mm[t] = m
    CC[t] = C
 }
 return(list(m=mm,C=CC))
} 

BF = function(y,theta,m0,C0,M){
  alpha = theta[1]
  beta  = theta[2]
  sig2  = theta[3]
  tau2  = theta[4]
  n     = length(y)
  x     = rnorm(M,m0,sqrt(C0))
  xs    = matrix(0,n,M)
  for (t in 1:n){
  	  x      = rnorm(M,alpha+beta*x,tau)
  	  w      = dnorm(y[t],x,sig)
    k      = sample(1:M,size=M,prob=w,replace=T)
    x      = x[k]
    xs[t,] = x 
  }
  return(xs)
}

APF = function(y,theta,m0,C0,M){
  alpha = theta[1]
  beta  = theta[2]
  sig2  = theta[3]
  tau2  = theta[4]
  n     = length(y)
  x     = rnorm(M,m0,sqrt(C0))
  xs    = matrix(0,n,M)
  for (t in 1:n){
  	  a      = alpha+beta*x
  	  w      = dnorm(y[t],a,sig)
    k      = sample(1:M,size=M,prob=w,replace=T)
  	  x      = rnorm(M,a[k],tau)
  	  w1     = dnorm(y[t],alpha+beta*x,sig)/w[k]
    x      = x[k]
    xs[t,] = x 
  }
  return(xs)
}

OBF = function(y,theta,m0,C0,M){
  alpha = theta[1]
  beta  = theta[2]
  sig2  = theta[3]
  tau2  = theta[4]
  n  = length(y)
  w2 = sig2+tau2
  A  = tau2/w2
  s  = sqrt(w2)
  d  = sqrt(A*sig2)
  x  = rnorm(M,m0,sqrt(C0))
  xs = matrix(0,n,M)
  for (t in 1:n){
  	  a      = alpha+beta*x
    m      = (1-A)*a+A*y[t]
    x1     = rnorm(M,m,d)
    w      = dnorm(y[t],a,s)
    k      = sample(1:M,size=M,prob=w,replace=T)
    x      = x1[k]
    xs[t,] = x 
  }
  return(xs)
}

OAPF = function(y,theta,m0,C0,M){
  n  = length(y)
  w2 = sig2+tau2
  A  = tau2/w2
  s  = sqrt(w2)
  d  = sqrt(A*sig2)
  x  = rnorm(M,m0,sqrt(C0))
  xs = matrix(0,n,M)
  for (t in 1:n){
  	  a = alpha+beta*x
    w = dnorm(y[t],a,s)
    k = sample(1:M,size=M,prob=w,replace=T)
    m = (1-A)*a[k]+A*y[t]
    x = rnorm(M,m,d)
    xs[t,] = x 
  }
  return(xs)
}

# SIMULATION EXERCISE
# -------------------
set.seed(1234)
ndat   = 20
nrep   = 20
M      = 1000
#n      = 100
n      = 1000
sig    = 1.0
taus   = c(0.05,0.25,0.75)
alpha  = 0.05
beta   = 0.95
x0     = alpha/(1-beta)
m0     = x0
C0     = 10
sig2   = sig^2
ntau   = length(taus)
qdlms  = array(0,c(ndat,ntau,n,3))
qbfs   = array(0,c(ndat,ntau,nrep,n,3))
qobfs  = array(0,c(ndat,ntau,nrep,n,3))
qapfs  = array(0,c(ndat,ntau,nrep,n,3))
qoapfs = array(0,c(ndat,ntau,nrep,n,3))
for (ll in 1:ndat){
  for (j in 1:ntau){
	  print(c(ll,j))
    tau   = taus[j]
    tau2  = tau^2
    theta = c(alpha,beta,sig2,tau2)
    sim   = DLM.SIM(n,theta,x0)
    y     = sim$y
    x     = sim$x
    dlm   = DLM(y,theta,m0,C0)
    qdlms[ll,j,,] = cbind(dlm$m+qnorm(0.025)*sqrt(dlm$C),dlm$m,dlm$m+qnorm(0.975)*sqrt(dlm$C))
    for (i in 1:nrep){
      qbfs[ll,j,i,,]   = t(apply(BF(y,theta,m0,C0,M),1,quantile,c(0.025,0.5,0.975)))
      qapfs[ll,j,i,,]  = t(apply(APF(y,theta,m0,C0,M),1,quantile,c(0.025,0.5,0.975)))
      qobfs[ll,j,i,,]  = t(apply(OBF(y,theta,m0,C0,M),1,quantile,c(0.025,0.5,0.975)))
      qoapfs[ll,j,i,,] = t(apply(OAPF(y,theta,m0,C0,M),1,quantile,c(0.025,0.5,0.975)))
    }
  }
}

# COMPUTING ROOT MEAN SQUARE ERRORS AND MEAN ABSOLUTE ERRORS
# ----------------------------------------------------------
rmse   = array(0,c(ndat,ntau,3,4,nrep))
rrmse  = array(0,c(ndat,ntau,3,3,nrep))
mae    = array(0,c(ndat,ntau,3,4,nrep))
rmae   = array(0,c(ndat,ntau,3,3,nrep))
for (ll in 1:ndat)
for (i in 1:ntau)
  for (j in 1:3){
    rmse1=qbfs[ll,i,,,j]-matrix(qdlms[ll,i,,j],nrep,n,byrow=TRUE)
    rmse2=qapfs[ll,i,,,j]-matrix(qdlms[ll,i,,j],nrep,n,byrow=TRUE)
    rmse3=qobfs[ll,i,,,j]-matrix(qdlms[ll,i,,j],nrep,n,byrow=TRUE)
    rmse4=qoapfs[ll,i,,,j]-matrix(qdlms[ll,i,,j],nrep,n,byrow=TRUE)
    rmse[ll,i,j,1,]=sqrt(apply(rmse1^2,1,mean))
    rmse[ll,i,j,2,]=sqrt(apply(rmse2^2,1,mean))
    rmse[ll,i,j,3,]=sqrt(apply(rmse3^2,1,mean))
    rmse[ll,i,j,4,]=sqrt(apply(rmse4^2,1,mean))
    mae[ll,i,j,1,]=apply(abs(rmse1),1,mean)
    mae[ll,i,j,2,]=apply(abs(rmse2),1,mean)
    mae[ll,i,j,3,]=apply(abs(rmse3),1,mean)
    mae[ll,i,j,4,]=apply(abs(rmse4),1,mean)
    rrmse[ll,i,j,1,] = rmse[ll,i,j,1,]/rmse[ll,i,j,4,]
    rrmse[ll,i,j,2,] = rmse[ll,i,j,2,]/rmse[ll,i,j,4,]
    rrmse[ll,i,j,3,] = rmse[ll,i,j,3,]/rmse[ll,i,j,4,]
    rmae[ll,i,j,1,]  = mae[ll,i,j,1,]/mae[ll,i,j,4,]
    rmae[ll,i,j,2,]  = mae[ll,i,j,2,]/mae[ll,i,j,4,]
    rmae[ll,i,j,3,]  = mae[ll,i,j,3,]/mae[ll,i,j,4,]
}

quants = c("2.5th","50th","97.5th")
filters = c("BF","APF","OBF","OAPF")

pdf(file="longrun2-a.pdf",width=20,height=12)
par(mfrow=c(3,ntau))
for (j in 1:3)
  for (i in 1:ntau){
    boxplot(as.numeric(rmse[,i,j,1,]),as.numeric(rmse[,i,j,2,]),as.numeric(rmse[,i,j,3,]),
            as.numeric(rmse[,i,j,4,]),names=filters,outline=FALSE,ylab="Root MSE")
    title(paste("tau=",taus[i]," - ",quants[j]," percentile",sep=""))
  }
dev.off()


pdf(file="longrun2-b.pdf",width=20,height=12)
par(mfrow=c(3,ntau))
for (j in 1:3)
  for (i in 1:ntau){
    boxplot(as.numeric(mae[,i,j,1,]),as.numeric(mae[,i,j,2,]),as.numeric(mae[,i,j,3,]),
            as.numeric(mae[,i,j,4,]),names=filters,outline=FALSE,ylab="MAE")
    title(paste("tau=",taus[i]," - ",quants[j]," percentile",sep=""))
  }
dev.off()

pdf(file="longrun2-c.pdf",width=20,height=12)
par(mfrow=c(3,ntau))
for (j in 1:3)
  for (i in 1:ntau){
    boxplot(as.numeric(rrmse[,i,j,1,]),as.numeric(rrmse[,i,j,2,]),as.numeric(rrmse[,i,j,3,]),
            names=filters[1:3],outline=FALSE,ylab="Relative Root MSE")
    title(paste("tau=",taus[i]," - ",quants[j]," percentile",sep=""))
    abline(h=1,col=2)
  }
dev.off()

pdf(file="longrun2-d.pdf",width=20,height=12)
par(mfrow=c(3,ntau))
for (j in 1:3)
  for (i in 1:ntau){
    boxplot(as.numeric(rmae[,i,j,1,]),as.numeric(rmae[,i,j,2,]),as.numeric(rmae[,i,j,3,]),
            names=filters[1:3],outline=FALSE,ylab="Relative MAE")
    title(paste("tau=",taus[i]," - ",quants[j]," percentile",sep=""))
    abline(h=1,col=2)
  }
dev.off()