[R] Setting up a State Space Model in dlm

[Giovanni Petris]

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] <
> ml-node+3770873-1409921251-254102 at n4.nabble.com> 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<http://r.789695.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3580664&code=cmFobHV3YWxpYUBnbWFpbC5jb218MzU4MDY2NHwtMjAzMjI4MTk2NQ==>.
> >
> >
> 
> 
>