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