rfirst.order <- function ( n, V, W, start ) {
   ### n: length of time series
   ### V: observational variances
   ### W: evolution variances
   ### start: starting level

  mu _ start + cumsum ( rnorm ( n, sd=sqrt(W) ) )
  return ( list ( level=mu, observation = mu + rnorm ( n, sd=sqrt(V) ) ) )
}

examp1 _ rfirst.ord ( 200, 1, .05, 25 )
examp2 _ rfirst.ord ( 200, 1, .5, 25 )
examp3 _ rfirst.ord ( 200, 1, 5, 25 )

par ( mfrow=c(3,1) )
tsplot ( examp1$obs, examp1$level, lty=1:2 )
tsplot ( examp2$obs, examp2$level, lty=1:2 )
tsplot ( examp3$obs, examp3$level, lty=1:2 )

par ( mfrow=c(3,2) )
tsplot ( bonds.yield[,1], main="daily bond yields", xlab="day",
  ylab="yield" )
tsplot( co2, main="Mauna Loa CO2 concentration", xlab="Year",
  ylab="concentration" )
tsplot ( tsmatrix ( corn.rain, corn.yield ),
  main="Rainfall and Corn Yield", xlab="Year",
  ylab="rain (inches), yield (bushels per acre)" )
tsplot ( lynx, main="Annual Number of Lynx Trappings", xlab="Year",
  ylab="Number of Lynx" )
tsplot ( ozone.median, main="Ozone levels",
  sub="Median daily maxima for June to August", xlab="Day" )
tsplot ( sunspots, main="Monthly Mean Relative Sunspot Numbers",
  xlab="Year" )

### which might be modeled by a first-order polynomial? ###
### bond yields, ozone levels and crop yields look like fluctuations
### around a slowly changing level.  The others have cycles which
### cannot be modeled yet.
par ( mfrow=c(2,3) )
tsplot ( examp1$obs, main="V/W = 20" )
tsplot ( examp2$obs, main="V/W = 2" )
tsplot ( examp3$obs, main="V/W = 0.2" )
tsplot ( bonds.yield[,1], main="daily bond yields", xlab="day",
  ylab="yield" )
tsplot ( tsmatrix ( corn.rain, corn.yield ),
  main="Rainfall and Corn Yield", xlab="Year",
  ylab="rain (inches), yield (bushels per acre)" )
tsplot ( ozone.median, main="Ozone levels",
  sub="Median daily maxima for June to August", xlab="Day" )
### V/W=0.2 looks good for bond yields
### replot ;with different lengths for ozone and crop
tsplot ( examp1$obs[1:40], main="V/W = 20" )
tsplot ( examp2$obs[1:40], main="V/W = 2" )
tsplot ( examp3$obs[1:40], main="V/W = 0.2" )
tsplot ( bonds.yield[,1], main="daily bond yields", xlab="day",
  ylab="yield" )
tsplot ( tsmatrix ( corn.yield ),
  main="Corn Yield", xlab="Year",
  ylab="rain (inches), yield (bushels per acre)" )
tsplot ( ozone.median, main="Ozone levels",
  sub="Median daily maxima for June to August", xlab="Day" )

update.dlm _ function ( Y, V, W, start ) {
  n <- length(y)
  m _ rep ( NA, length=n )
  C _ rep ( NA, length=n )
  R _ rep ( NA, length=n )
  Q _ rep ( NA, length=n )
  f _ rep ( NA, length=n )
  A _ rep ( NA, length=n )
  e _ rep ( NA, length=n )

  Y _ c(NA,Y)
  
  C[1] <- start$C
  m[1] <- start$m

  for ( t in 2:n ) {
    R[t] <- C[t-1] + W[t]
    f[t] <- m[t-1]
    Q[t] <- R[t] + V[t]
    A[t] <- R[t] / Q[t]
    e[t] <- Y[t] - f[t]
    m[t] <- m[t-1] + A[t]*e[t]
    C[t] <- A[t]*V[t]
  }

  return ( list ( m=m, C=C, R=R, f=f, Q=Q ) )
}

    ### KURIT example
    y _ c ( 150, 136, 143, 154, 135, 148, 128, 149, 146 )
    V _ rep ( 100, length(y) )
    W _ rep ( 5, length(y) )
    start _ list ( m=130, C=400 )
    result _ update.dlm ( y, V, W, start )

conj.update.dlm _ function ( Y, W=NA, delta, m.0, C.0, n.0, S.0 ) {
  N <- length(y)
  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 )

  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
       ### to specify the model via a discount factor instead of W
       ### we use the following line
    if ( missing(W) & !missing(delta) ) W[t] _ C[t-1] * (1-delta) / delta
    R[t] <- C[t-1] + W[t]
    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]
  }

  return ( list ( m=m, C=C, R=R, f=f, Q=Q, n=n, S=S ) )
}

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 )

y _ ts ( y, start=1975, frequency=12 )

par(mfrow=c(1,1))
tsplot(y)
fit.low _ lowess ( time(y), y )
lines ( fit.low )

tsplot ( y, type="p", pch="*", main="Estimated level of exchange rate" )
for ( d in 1:4) {
  lines ( time(y),
          conj.update.dlm (y,delta=.6+.1*d,m.0=0,C.0=1,n.0=1,S.0=.01)$m,
          lty=d )
}
legend ( 1981, 6, legend=paste(delta=.6+.1*(1:4)), lty=1:4, bty="n" )
text ( 1981, 6.5, "Discount factor", adj=0 )