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-254...@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==>.
>
>
--
Ram Ahluwalia
--
View this message in context:
http://r.789695.n4.nabble.com/Setting-up-a-State-Space-Model-in-dlm-tp3580664p3770976.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
______________________________________________
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.
