Setting up a State Space Model in dlm
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.