# Setting up a State Space Model in dlm Classic List Threaded 9 messages Open this post in threaded view
|

## 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.htmlI 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.htmlI 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_tfor 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), 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), 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 ))   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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Setting up a State Space Model in dlm

Open this post in threaded view
|

## Re: Setting up a State Space Model in dlm

 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Setting up a State Space Model in dlm

 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/paperThe 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.
Open this post in threaded view
|

## Re: Setting up a State Space Model in dlm

 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Setting up a State Space Model in dlm

 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] 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-helpPLEASE 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
Open this post in threaded view
|