# This version fits the 1st order polynomial DLM
#  Computes (theta_t|D_t) and (theta_t|D_N)
conj.update.dlm=function ( Y, delta, m.0, C.0, n.0, S.0 ) 
{
  N <- length(y)+1
  m = rep ( NA, length=N )
  n = rep ( NA, length=N )
  C = rep ( NA, length=N )
  R = rep ( NA, length=N )
  Q = rep ( NA, length=N )
  S = rep ( NA, length=N )
  f = rep ( NA, length=N )
  A = rep ( NA, length=N )
  e = rep ( NA, length=N )
#
  aRet= rep(NA,length=N)
  RRet= rep(NA,length=N)

  Y = c(NA,Y)
  
  C[1] <- C.0
  m[1] <- m.0
  S[1] <- S.0
  n[1] <- n.0

  for ( t in 2:N ) {
    n[t] <- n[t-1] + 1
    R[t] <- C[t-1]/delta
    f[t] <- m[t-1]
    Q[t] <- R[t] + S[t-1]
    A[t] <- R[t] / Q[t]
    e[t] <- Y[t] - f[t]
    S[t] <- S[t-1] + ( S[t-1]/n[t] ) * ( e[t]^2/Q[t] - 1 )
    m[t] <- m[t-1] + A[t]*e[t]
    C[t] <- A[t]*S[t]
  }
# Filtering recurrances.  This version
# only valid for G=I and discount factor delta.
# In this case, a_t= m_{t-1}.  See exercise 4.11 (d)
   aRet[N]=m[N]
   RRet[N]=C[N]
  for (t in 1:(N-2)){
    aRet[N-t]= m[N-t] + delta*(aRet[N-t+1]-m[N-t])
    RRet[N-t]= C[N-t] + delta^2*(RRet[N-t+1]-R[N-t+1])
     }   
return(list(m=m,C=C,R=R,f=f,Q=Q,n=n,S=S,aRet=aRet,RRet=RRet))
 }
# Exchange rate data
y = c ( 1.35,  1.00, -1.96, -2.17, -1.78, -4.21,
         -3.30, -1.43, -1.35, -0.34, -1.38,  0.30,
         -0.10, -4.13, -5.12, -2.13, -1.17, -1.24,
         -0.18,  0.29,  0.23,  0.00,  0.06,  0.17,
          0.98,  0.17,  1.59,  2.62,  1.96,  4.28,
          0.26, -1.66, -3.03, -1.80,  1.04,  3.06,
         -0.05,  1.68,  1.70, -0.73,  2.59,  6.77,
         -0.98, -1.71, -2.53, -0.61,  3.14,  2.96,
          1.01, -3.69,  0.45,  3.89,  1.38,  1.57,
         -0.08,  1.30,  0.62, -0.87, -2.11,  2.48,
         -4.73, -2.70, -2.45, -4.17, -5.76, -5.09,
         -2.92, -0.22,  1.42,  3.26,  0.05, -0.95,
         -2.14, -2.19, -1.96,  2.18, -2.97, -1.89,
          0.12, -0.76, -0.94, -3.90, -0.86, -2.88,
         -2.58, -2.78,  3.30,  2.06, -1.54, -1.30,
         -1.78, -0.13, -0.20, -1.35, -2.82, -1.97,
          2.25,  1.17, -2.29, -2.49, -0.87, -4.15, -0.53 )
#
# Plot of data with posterior means m_t
#
plot(y, type="p", pch="*",ylab=" ", main="Estimated level of exchange rate")
for ( d in 1:4) {
  lines(conj.update.dlm (y,delta=.6+.1*d,m.0=0,C.0=1,n.0=1,S.0=1.0)$m,lty=d )
}
legend(80,6, legend=paste(delta=.6+.1*(1:4)), lty=1:4, bty="n" )
text (80, 6.5, "Discount factor", adj=0 )
#
# Plot of data with posterior means a_N(-k)
#
plot(y, type="p", pch="*",ylab=" ", main="Estimated level of exchange rate (Restrospective)" )
for ( d in 1:4) {
  lines(conj.update.dlm (y,delta=.6+.1*d,m.0=0,C.0=1,n.0=1,S.0=1.0)$aRet,lty=d )
}
legend (80,6,legend=paste(delta=.6+.1*(1:4)), lty=1:4, bty="n" )
text(80, 6.5, "Discount factor", adj=0 )
#
# Comparison of m_t and a_N(-t)
plot(y, type="p", pch="*",ylab=" ", main="Estimated level of exchange rate" )
rs=conj.update.dlm(y,delta=0.8,	m.0=0,C.0=1,n.0=1,S.0=1)
lines(rs$m,col=2)
lines(rs$aRet,col=4)
legend(80,6,legend=c("forward","retrospective"),lty=1,col=c(2,4),bty="n")
# Prob. interval around retrospective mean.  Use Corollary #4.3
# (theta_t|D_N) ~ t_nt(a_N(-t), S_N/S_t R_N(-t) )
m=rs$aRet
N=length(y)+1
C=(rs$S[N]/rs$S)*rs$RRet
plot(y,type='p',pch="*",ylab=" ")
lines(rs$aRet,col=4)
lines(rs$aRet + qt(.975,rs$n)*sqrt(C),lty=2,col=6)
lines(rs$aRet - qt(.975,rs$n)*sqrt(C),lty=2,col=6)