Setting up a State Space Model in dlm

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Setting up a State Space Model in dlm

Michael Ash-2
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))

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

Gavin Simpson
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))
>
> ______________________________________________
> [hidden email] 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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

Michael Ash-2
Thank you.    The new article linked below is excellent.  BTW, I have
had trouble getting the sample code to run.    In particular, there is
a variable k that does not appear to be set when it is used at line
227.

Here is the line where it is called:
fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS",
  control = list(maxit = 500))

Here is the error message:
$ R --vanilla < kfas.R > kfas.Rlog
Error in rep.int(1, length(par)) : object 'k' not found
Calls: optim -> rep.int
Execution halted

Any advice?  Thanks.

Best,
Michael





On Fri, Jun 10, 2011 at 8:22 AM, Gavin Simpson <[hidden email]> wrote:

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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

quantguy
I get the identical error even when applying the sample code from the dlm vignette on state space models: http://www.jstatsoft.org/v41/i04/paper

The sample code is here:

tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T), start=c(1978,1), frequency=12) * 100
y <- tmp[,1:4] - tmp[,"RKFREE"]
colnames(y) <- colnames(tmp)[1:4]
market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
rm("tmp")
m <- NCOL(y)


Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
dim(Zt) <- c(m, m, length(market))
Rt <- diag(nr = m)
logLik <- function(theta) {
   a <- diag(exp(0.5 * theta[1:m]), nr = m)
   a[upper.tri(a)] <- theta[(m + 1):k]
   Ht <- crossprod(a)
   a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
   a[upper.tri(a)] <- theta[-(1:(k + m))]
   Qt <- crossprod(a)
   lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
   Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
   P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
   return(-lik$lik)
  }

fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS", control = list(maxit = 500))
fit$conv

The parameter K is not defined. If I select an arbitrary value of K (looks like it is just a parameter to initialize optimization) I get a separate error.

Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

Giovanni Petris

Oops..
You need to add the following line, right after the "m <- NCOL(y)"
statement:

k <- m * (m+1) / 2

'k' is the number of independent parameters in an m-by-m covariance
matrix, and two such matrices are estimated (Ht and Qt, both time
invariant). Hence the unknown parameter theta has length 2k.

Hope this clarifies the issue.

Best,
Giovanni Petris

On Thu, 2011-08-25 at 19:30 -0700, quantguy wrote:

> I get the identical error even when applying the sample code from the dlm
> vignette on state space models: http://www.jstatsoft.org/v41/i04/paper
>
> The sample code is here:
>
> tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T),
> start=c(1978,1), frequency=12) * 100
> y <- tmp[,1:4] - tmp[,"RKFREE"]
> colnames(y) <- colnames(tmp)[1:4]
> market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
> rm("tmp")
> m <- NCOL(y)
>
>
> Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
> dim(Zt) <- c(m, m, length(market))
> Rt <- diag(nr = m)
> logLik <- function(theta) {
>    a <- diag(exp(0.5 * theta[1:m]), nr = m)
>    a[upper.tri(a)] <- theta[(m + 1):k]
>    Ht <- crossprod(a)
>    a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
>    a[upper.tri(a)] <- theta[-(1:(k + m))]
>    Qt <- crossprod(a)
>    lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
>    Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
>    P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
>    return(-lik$lik)
>   }
>
> fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS", control =
> list(maxit = 500))
> fit$conv
>
> The parameter K is not defined. If I select an arbitrary value of K (looks
> like it is just a parameter to initialize optimization) I get a separate
> error.
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3769863.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

quantguy
Thanks! However, I added the line of code and receive an error during the optim() procedure:
rror in fn(par, ...) : could not find function "kf"

The update code is here:

tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T), start=c(1978,1), frequency=12) * 100
y <- tmp[,1:4] - tmp[,"RKFREE"]
colnames(y) <- colnames(tmp)[1:4]
market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
rm("tmp")
m <- NCOL(y)
k <- m * (m+1) / 2

# 'k' is the number of independent parameters in an m-by-m covariance 
# matrix, and two such matrices are estimated (Ht and Qt, both time 
# invariant). Hence the unknown parameter theta has length 2k. 

Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
dim(Zt) <- c(m, m, length(market))
Rt <- diag(nr = m)
logLik <- function(theta) {
   a <- diag(exp(0.5 * theta[1:m]), nr = m)
   a[upper.tri(a)] <- theta[(m + 1):k]
   Ht <- crossprod(a)
   a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
   a[upper.tri(a)] <- theta[-(1:(k + m))]
   Qt <- crossprod(a)
   lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
   Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
   P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
   return(-lik$lik)
  }

fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS", control = list(maxit = 500))
fit$conv

# With the MLEs of the unknown parameters one can compute the smoothing estimates of the
# betas (Figure 11), as the following code illustrates

theta <- fit$par
a <- diag(exp(0.5 * theta[1:m]), nr = m)
a[upper.tri(a)] <- theta[(m+1):k]
Ht <- crossprod(a)
a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
a[upper.tri(a)] <- theta[-(1:(k + m))]
Qt <- crossprod(a)
smoothCAPM <- ks(kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt,
Ht = Ht, Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
P1inf = diag(rep(1, m))))
betas <- ts(t(smoothCAPM$ahat), start = start(market),
freq = frequency(market))



On Fri, Aug 26, 2011 at 9:41 AM, Giovanni Petris [via R] <[hidden email]> wrote:

Oops..
You need to add the following line, right after the "m <- NCOL(y)"
statement:

k <- m * (m+1) / 2

'k' is the number of independent parameters in an m-by-m covariance
matrix, and two such matrices are estimated (Ht and Qt, both time
invariant). Hence the unknown parameter theta has length 2k.

Hope this clarifies the issue.

Best,
Giovanni Petris

On Thu, 2011-08-25 at 19:30 -0700, quantguy wrote:

> I get the identical error even when applying the sample code from the dlm
> vignette on state space models: http://www.jstatsoft.org/v41/i04/paper
>
> The sample code is here:
>
> tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T),
> start=c(1978,1), frequency=12) * 100
> y <- tmp[,1:4] - tmp[,"RKFREE"]
> colnames(y) <- colnames(tmp)[1:4]
> market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
> rm("tmp")
> m <- NCOL(y)
>
>
> Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
> dim(Zt) <- c(m, m, length(market))
> Rt <- diag(nr = m)
> logLik <- function(theta) {
>    a <- diag(exp(0.5 * theta[1:m]), nr = m)
>    a[upper.tri(a)] <- theta[(m + 1):k]
>    Ht <- crossprod(a)
>    a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
>    a[upper.tri(a)] <- theta[-(1:(k + m))]
>    Qt <- crossprod(a)
>    lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
>    Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
>    P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
>    return(-lik$lik)
>   }
>
> fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS", control =
> list(maxit = 500))
> fit$conv
>
> The parameter K is not defined. If I select an arbitrary value of K (looks
> like it is just a parameter to initialize optimization) I get a separate
> error.
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3769863.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.



If you reply to this email, your message will be added to the discussion below:
http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3770873.html
To unsubscribe from Setting up a State Space Model in dlm, click here.



--
Ram Ahluwalia


Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

Giovanni Petris
In reply to this post by Michael Ash-2

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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

Giovanni Petris
In reply to this post by quantguy
This was an example using package KFAS, which needs to be loaded before
running the code ('kf' is the function that computes Kalman filter in
package KFAS)

Giovanni

On Fri, 2011-08-26 at 07:11 -0700, quantguy wrote:

> Thanks! However, I added the line of code and receive an error during the
> optim() procedure:
>
> rror in fn(par, ...) : could not find function "kf"
>
>
> The update code is here:
>
>
> tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt",
> header=T), start=c(1978,1), frequency=12) * 100
> y <- tmp[,1:4] - tmp[,"RKFREE"]
> colnames(y) <- colnames(tmp)[1:4]
> market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
> rm("tmp")
> m <- NCOL(y)
> k <- m * (m+1) / 2
>
> # 'k' is the number of independent parameters in an m-by-m covariance
> # matrix, and two such matrices are estimated (Ht and Qt, both time
> # invariant). Hence the unknown parameter theta has length 2k.
>
> Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
> dim(Zt) <- c(m, m, length(market))
> Rt <- diag(nr = m)
> logLik <- function(theta) {
>    a <- diag(exp(0.5 * theta[1:m]), nr = m)
>    a[upper.tri(a)] <- theta[(m + 1):k]
>    Ht <- crossprod(a)
>    a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
>    a[upper.tri(a)] <- theta[-(1:(k + m))]
>    Qt <- crossprod(a)
>    lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
>    Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
>    P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
>    return(-lik$lik)
>   }
>
> fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS",
> control = list(maxit = 500))
> fit$conv
>
> # With the MLEs of the unknown parameters one can compute the
> smoothing estimates of the
> # betas (Figure 11), as the following code illustrates
>
> theta <- fit$par
> a <- diag(exp(0.5 * theta[1:m]), nr = m)
> a[upper.tri(a)] <- theta[(m+1):k]
> Ht <- crossprod(a)
> a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
> a[upper.tri(a)] <- theta[-(1:(k + m))]
> Qt <- crossprod(a)
> smoothCAPM <- ks(kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt,
> Ht = Ht, Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
> P1inf = diag(rep(1, m))))
> betas <- ts(t(smoothCAPM$ahat), start = start(market),
> freq = frequency(market))
>
>
>
>
> On Fri, Aug 26, 2011 at 9:41 AM, Giovanni Petris [via R] <
> [hidden email]> wrote:
>
> >
> > Oops..
> > You need to add the following line, right after the "m <- NCOL(y)"
> > statement:
> >
> > k <- m * (m+1) / 2
> >
> > 'k' is the number of independent parameters in an m-by-m covariance
> > matrix, and two such matrices are estimated (Ht and Qt, both time
> > invariant). Hence the unknown parameter theta has length 2k.
> >
> > Hope this clarifies the issue.
> >
> > Best,
> > Giovanni Petris
> >
> > On Thu, 2011-08-25 at 19:30 -0700, quantguy wrote:
> >
> > > I get the identical error even when applying the sample code from the dlm
> >
> > > vignette on state space models: http://www.jstatsoft.org/v41/i04/paper
> > >
> > > The sample code is here:
> > >
> > > tmp <- ts(read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T),
> >
> > > start=c(1978,1), frequency=12) * 100
> > > y <- tmp[,1:4] - tmp[,"RKFREE"]
> > > colnames(y) <- colnames(tmp)[1:4]
> > > market <- tmp[,"MARKET"] - tmp[,"RKFREE"]
> > > rm("tmp")
> > > m <- NCOL(y)
> > >
> > >
> > > Zt <- sapply(seq_along(market), function(i) market[i] %x% diag(m))
> > > dim(Zt) <- c(m, m, length(market))
> > > Rt <- diag(nr = m)
> > > logLik <- function(theta) {
> > >    a <- diag(exp(0.5 * theta[1:m]), nr = m)
> > >    a[upper.tri(a)] <- theta[(m + 1):k]
> > >    Ht <- crossprod(a)
> > >    a <- diag(exp(0.5 * theta[1:m + k]), nr = m)
> > >    a[upper.tri(a)] <- theta[-(1:(k + m))]
> > >    Qt <- crossprod(a)
> > >    lik <- kf(yt = t(y), Zt = Zt, Tt = diag(nr = m), Rt = Rt, Ht = Ht,
> > >    Qt = Qt, a1 = rep(0, m), P1 = matrix(0, m, m),
> > >    P1inf = diag(rep(1, m)), optcal = c(FALSE, FALSE, FALSE, FALSE))
> > >    return(-lik$lik)
> > >   }
> > >
> > > fit <- optim(par = rep(0, 2 * k), fn = logLik, method = "BFGS", control =
> >
> > > list(maxit = 500))
> > > fit$conv
> > >
> > > The parameter K is not defined. If I select an arbitrary value of K
> > (looks
> > > like it is just a parameter to initialize optimization) I get a separate
> > > error.
> > >
> > >
> > >
> > > --
> > > View this message in context:
> > http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3769863.html
> > > Sent from the R help mailing list archive at Nabble.com.
> > >
> > > ______________________________________________
> > > [hidden email] <http://user/SendEmail.jtp?type=node&node=3770873&i=0>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.
> >
> > ______________________________________________
> > [hidden email] <http://user/SendEmail.jtp?type=node&node=3770873&i=1>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.
> >
> >
> > ------------------------------
> >  If you reply to this email, your message will be added to the discussion
> > below:
> >
> > http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3770873.html
> >  To unsubscribe from Setting up a State Space Model in dlm, click here<
> >
> >
>
>
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Setting up a State Space Model in dlm

kch382001@yahoo.com
This post has NOT been accepted by the mailing list yet.
In reply to this post by Michael Ash-2
Hi, Michael:

   I am working on similar model. Have you get your DLM model setup to run the model? Do you mind sharing the setup code for DLM?
   Thanks.