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