############################################################################################
# 
# 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/
#
############################################################################################
set.seed(123456)

# Simulating 1st order DLM
# ------------------------
W = 0.5
V = 1.0
n = 100
x0 = 0.00
w  = rnorm(n,0,sqrt(W))
v  = rnorm(n,0,sqrt(V))
x  = y = rep(0,n)
x[1] = x0   + w[1]
y[1] = x[1] + v[1]
for (t in 2:n){
  x[t] = x[t-1] + w[t]
  y[t] = x[t] + v[t]
}
y1 = y
x1 = x

# Online fitting 1st order DLM
# -----------------------------
  y  = y1
  a1 = 0
  R1 = 100
  m  = rep(0,n)
  C  = rep(0,n)
  B  = rep(0,n-1)
  H  = rep(0,n-1)
  # time t=1
  a = a1
  R = R1
  f = a
  Q = R+V
  A = R/Q
  m[1] = a+A*(y[1]-f)
  C[1] = R-Q*A^2
  # forward filtering
  for (t in 2:n){
    a = m[t-1]
    R = C[t-1] + W
    f = a
    Q = R+V
    A = R/Q
    m[t] = a+A*(y[t]-f)
    C[t] = R-Q*A^2
}
m1 = m
C1 = C

# Smoothing the 1st order DLM
# ---------------------------
  y  = y1
  a1 = 0
  R1 = 100
  a  = rep(0,n)
  R  = rep(0,n)
  m  = rep(0,n)
  C  = rep(0,n)
  mm = rep(0,n-1)
  CC = rep(0,n-1)
  # time t=1
  a[1] = a1
  R[1] = R1
  f = a[1]
  Q = R[1]+V
  A = R[1]/Q
  m[1] = a[1]+A*(y[1]-f)
  C[1] = R[1]-Q*A^2
  # forward filtering
  for (t in 2:n){
    a[t] = m[t-1]
    R[t] = C[t-1] + W
    f = a[t]
    Q = R[t]+V
    A = R[t]/Q
    m[t] = a[t]+A*(y[t]-f)
    C[t] = R[t]-Q*A^2
  }
  mm[n] = m[n]
  CC[n] = C[n]
  for (t in (n-1):1){
    mm[t] = m[t] + C[t]/R[t+1]*(mm[t+1]-a[t+1])
    CC[t] = C[t] - (C[t]^2)/(R[t+1]^2)*(R[t+1]-CC[t+1])
  }

L1  = m-2*C
U1  = m+2*C
L2  = mm-2*CC
U2  = mm+2*CC

pdf(file="local-level-sim.pdf",width=10,height=7)
par(mfrow=c(1,1))
plot(y1,type="l",xlab="time",ylab="",lwd=3,ylim=c(min(L1,L2,y,x),max(U1,U2,y,x)))
lines(x1,col=2,lwd=3)
legend(70,8,legend=c("y(t)","x(t)"),lty=c(1,1),col=c(1,2),lwd=c(3,3),bty="n")
dev.off()

pdf(file="local-level-ffbs.pdf",width=10,height=7)
par(mfrow=c(1,1))
ts.plot(m,xlab="time",ylab="",ylim=c(min(L1,L2,y,x),max(U1,U2,y,x)),lwd=2)
lines(L1,lwd=2)
lines(U1,lwd=2)
lines(mm,col=2,lwd=2)
lines(L2,col=2,lwd=2)
lines(U2,col=2,lwd=2)
legend(70,8,legend=c("Forward filtering","Backward smoothing"),
       lty=c(1,1),col=c(1,2),lwd=c(3,3),bty="n")
dev.off()

N = 50
W = 0.5
V = 1.0
Vs = seq(0.3,1.5,length=N)
Ws = seq(0.1,1.5,length=N)
likes = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N){
  	V = Vs[i]
  	W = Ws[j]
  a1 = 0
  R1 = 100
  m  = rep(0,n)
  C  = rep(0,n)
  B  = rep(0,n-1)
  H  = rep(0,n-1)
  # time t=1
  a = a1
  R = R1
  f = a
  Q = R+V
  A = R/Q
  m[1] = a+A*(y[1]-f)
  C[1] = R-Q*A^2
  like = dnorm(y[1],f,sqrt(Q),log=TRUE)
  # forward filtering
  for (t in 2:n){
    a = m[t-1]
    R = C[t-1] + W
    f = a
    Q = R+V
    A = R/Q
    m[t] = a+A*(y[t]-f)
    C[t] = R-Q*A^2
    like = like + dnorm(y[t],f,sqrt(Q),log=TRUE)
 }
 likes[i,j] = like
}
pdf(file="local-level-contour.pdf",width=10,height=7)
contour(Vs,Ws,exp(likes),xlab=expression(sigma^2),ylab=expression(tau^2),drawlabels=FALSE)
dev.off()