[R] Setting up a State Space Model in dlm

[Michael Ash]

This question pertains to setting up a model in the package "dlm"
(dynamic linear models,
http://cran.r-project.org/web/packages/dlm/index.html

I have read both the vignette and "An R Package for Dynamic Linear
Models" (http://www.jstatsoft.org/v36/i12/paper), both of which are
very helpful.  There is also some discussion at
https://stat.ethz.ch/pipermail/r-help/2009-May/198463.html

I have what I think is a relatively straightforward state-space model
but am unable to translate it into the terms of dlm.   It would be
very helpful to get a basic dlm setup for the problem and I would
guess that I can then modify it with more lags, etc., etc.

The main equation is
pi[t] = a * pi[t-1] + b*(U[t] - UN[t]) + e[t]

(see http://chart.apis.google.com/chart?cht=tx&chl=%5Cpi_t=a%5Cpi_{t-1}%2bb%28U_t-U^N_{t}%29%2Be_t
for a pretty version)

with pi and U observed, a and b fixed coefficients, and e a
well-behaved error term (gaussian, say, variance unknown).
The object of interest is the unobserved and time-varying component UN
which evolves according to

UN[t] = UN[t-1] + w[t]

(see http://chart.apis.google.com/chart?cht=tx&chl=U%5EN_%7Bt%7D%20=%20U%5EN_%7Bt-1%7D%20%2B%20%5Cepsilon_t
for a pretty version)
that is, a random walk with well-behaved error term with var(w) known
(or assumed).

I'm interested in the estimates of a and b and also in estimating the
time series of UN.

Note that the term b*(U[t] - UN[t]) makes this a nonlinear model.

Below is code that does not work as expected.  I see the model as
having four parameters, a, b, var(e), and UN.  (Or do I have a
parameter UN[t] for every period?)

I do not fully understand the dlm syntax.  Is FF specified properly?
What should X look like?  How does m0 relate to parm()?


I would be grateful if someone would be willing to glance at the code.
 Thanks. Michael

library(quantmod)
library(dlm)
## Get and organize the data
getSymbols("UNRATE",src="FRED")  ## Unemployment rate
getSymbols("GDPDEF",src="FRED")  ## Quarterly GDP Implicit Price Deflator
u <- aggregate(UNRATE,as.yearqtr,mean)
gdpdef <- aggregate(GDPDEF,as.yearqtr,mean)
pi <- diff(log(gdpdef))*400
pilag <- lag(pi,-1)
tvnairu <- cbind(pi,pilag,u)
tvnairu.df <- subset(data.frame(tvnairu), !is.na(pi) & !is.na(u) &
!is.na(pilag))


## First attempt
buildNAIRU <- function(x) {
  modNAIRU <- dlm(FF=t(matrix(c(1,1,1,0))),
                  GG=diag(4),
                  W=matrix(c(0,0,0,0,  0,0,0,0,  0,0,0.04,0,  0,0,0,0),4,4),
                  V=exp(x[4]), m0=rep(0,4), C0=diag(1e07,4),
                  JFF = t(matrix(c(1,1,0,0))),
                  X=cbind( tvnairu.df$pilag, tvnairu.df$u))
  return(modNAIRU)
}

(fitNAIRU <- dlmMLE(tvnairu.df$pi, parm=c(0,0,0,0) , build=buildNAIRU,
hessian=TRUE, control=list(maxit=500)))
(dlmNAIRU <- buildNAIRU(fitNAIRU$par))


## Second attempt
buildNAIRU <- function(x) {
  modNAIRU <- dlm(FF=t(matrix(c(1,1,0,0))),
                  GG=diag(4),
                  W=matrix(c(0,0,0,0, 0,0,0,0, 0,0,0.04,0, 0,0,0,0 ),4,4),
                  V=exp(x[4]), m0=c(0.7,-0.3,5,1), C0=diag(100,4),
                  JFF = t(matrix(c(1,1,0,0))),
                  X=cbind( tvnairu.df$pilag, tvnairu.df$u - x[3] ))
  return(modNAIRU)
}

(fitNAIRU <- dlmMLE(tvnairu.df$pi, parm=c(0,0,0,0) , build=buildNAIRU,
hessian=TRUE, control=list(maxit=500)))
(dlmNAIRU <- buildNAIRU(fitNAIRU$par))