R help - Jun 2011 - Setting up a State Space Model in dlm

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))

On Tue, 2011-06-07 at 17:24 +0100, Michael Ash wrote:
> 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
The author of the dlm package has just published a paper on state space
models in R including details on setting up dlm:

http://www.jstatsoft.org/v41/i04

That might help with your question - I haven't seen a reply on list, but
am unable to help answer it either.

HTH

G
> 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))
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

I just saw this old post, but it seems that nobody replied, so let me try.

If you can assume that also U[t] evelves as a random walk, I would build a
DLM by taking the state vector to be 

x[t] = (U[t], UN[t], pi[t])'

By plugging in the equation for pi[t] the random walk expressions for U[t]
and UN[t] you get the system equation of the DLM. The observation matrix F
will just be the 2 by 3 matrix that extracts components 2 and 3 from the
3-dimensional state. Set the observation variance V to a tiny multiple of
the 2 by 2 identity matrix, so that it is invertible but practically
negligible. 

I have tried some code to implement this model - here it is:

myMod <- dlm(FF = matrix(c(0, 1, 0, 0, 0, 1), 2, 3, TRUE),
             GG = matrix(c(1, 0, 0, 0, 1, 0, NA, NA, NA), 3, 3, TRUE),
             W = diag(3), # will change this in 'build' function
             V = diag(1e-7, 2),
             m0 = rep(0, 3),
             C0 = diag(1e8, 3))
R <- matrix(c(1, 0, 0, 0, 1, 0, NA, NA, 1), 3, 3, TRUE)
buildFun <- function(theta) {
    ## theta[1] : 'b'
    ## theta[2] : 'a'
    ## theta[3] : log innovation std dev of U
    ## theta[4] : log innovation std dev of UN
    ## theta[5] : log innovation std dev of pi
    GG(myMod)[3, ] <- theta[c(1, 1, 2)] * c(1, -1, 1)
    R[3, 1 : 2] <- theta[1] * c(1, -1)
    dd <- exp(theta[3 : 5])
    W(myMod) <- tcrossprod(R * rep(dd, each = 3))
    return(myMod)
}

outMLE <- dlmMLE(y = tvnairu[, c("pi", "u")], parm = c(1,
1, 0, 0, 0),
                 build = buildFun, lower = c(-Inf, -Inf, rep(-8, 3)),
                 upper = c(Inf, Inf, rep(12, 3)), 
                 control = list(trace = 1, REPORT = 5, maxit = 1000))
outMLE$par

In the estimates I get, the 'b' parameter is tiny, but this may be a
local
optimum - you need to try different starting values for the optimizer.

Best,
Giovanni Petris



Michael Ash-2 wrote:
> 
> 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))
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

--
View this message in context:
http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3771129.html
Sent from the R help mailing list archive at Nabble.com.