[R] Setting up a State Space Model in dlm

[Giovanni Petris]

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