R/DLM.R

Defines functions dlmGibbsDIG dropFirst ergMean print.mcmcMean mcmcMean mcmcSD dlmRandom tsdiag.dlmFiltered residuals.dlmFiltered dlmForecast dlmSvd2var rwishart dlmBSample dlmSmooth.dlmFiltered dlmSmooth.default dlmSmooth dlmFilter dlmLL print.dlm dlmModARMA dlmModTrig dlmModSeas dlmModPoly dlmModReg X.dlm JW.dlm JV.dlm JGG.dlm JFF.dlm X JW JV JGG JFF m0.dlm C0.dlm W.dlm V.dlm GG.dlm FF.dlm m0 C0 W V GG FF is.dlmFiltered as.dlm is.dlm dlm ARtransPars bdiag

Documented in ARtransPars as.dlm bdiag C0 C0.dlm dlm dlmBSample dlmFilter dlmForecast dlmGibbsDIG dlmLL dlmModARMA dlmModPoly dlmModReg dlmModSeas dlmModTrig dlmRandom dlmSmooth dlmSmooth.default dlmSmooth.dlmFiltered dlmSvd2var dropFirst ergMean FF FF.dlm GG GG.dlm is.dlm JFF JFF.dlm JGG JGG.dlm JV JV.dlm JW JW.dlm m0 m0.dlm mcmcMean mcmcSD residuals.dlmFiltered rwishart V V.dlm W W.dlm X X.dlm

###### Function to create block diagonal matrices (from R-help, I suppose)
bdiag <- function(...)
{
    if (nargs() == 1)
        x <- as.list(...)
    else
        x <- list(...)
    n <- length(x)
    if(n==0) return(NULL)
    x <- lapply(x, function(y) if(length(y)) as.matrix(y) else
                stop("Zero-length component in x"))
    d <- array(unlist(lapply(x, dim)), c(2, n))
    rr <- d[1,]
    cc <- d[2,]
    rsum <- sum(rr)
    csum <- sum(cc)
    out <- array(0, c(rsum, csum))
    ind <- array(0, c(4, n))
    rcum <- cumsum(rr)
    ccum <- cumsum(cc)
    ind[1,-1] <- rcum[-n]
    ind[2,] <- rcum
    ind[3,-1] <- ccum[-n]
    ind[4,] <- ccum
    imat <- array(1:(rsum * csum), c(rsum, csum))
    iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2],
                                                 (y[3]+1):y[4]], imat=imat)
    iuse <- as.vector(unlist(iuse))
    out[iuse] <- unlist(x)
    return(out)
}

###### From R^p to the AR parameters of a stationary AR(p) (univariate)
###### For the multivariate analog, see Ansley and Kohn, 1986,
###### J. Statist. Comput. Simul. 24.
ARtransPars <- function(raw) {
    p <- length(raw)
    return(.Call(C_ARtranspar, p, as.double(raw), PACKAGE="dlm"))
}

###### Constructor for dlm objects
dlm <- function(...) {
    if (nargs() == 1)
        x <- as.list(...)
    else
        x <- list(...)
    ## required components
    nm <- c("m0", "C0", "FF", "V", "GG", "W")
    nmInd <- match(nm, names(x))
    if (any(is.na(nmInd)))
        stop(paste("Component(s)", paste(nm[is.na(nmInd)], collapse=", "),
                   "is (are) missing"))
    x[nmInd[-1]] <- lapply(x[nmInd[-1]], as.matrix)
    if (!is.numeric(x$FF))
        stop("Component FF must be numeric")
    m <- nrow(x$FF)
    p <- ncol(x$FF)
    if (!is.numeric(x$V))
        stop("Component V must be numeric")
    if (!( nrow(x$V) == m && ncol(x$V) == m))
        stop("Incompatible dimensions of matrices")
    if (!is.numeric(x$GG))
        stop("Component GG must be numeric")
    if (!( nrow(x$GG) == p && ncol(x$GG) == p))
        stop("Incompatible dimensions of matrices")
    if (!is.numeric(x$W))
        stop("Component W must be numeric")
    if (!( nrow(x$W) == p && ncol(x$W) == p))
        stop("Incompatible dimensions of matrices")
    if (!is.numeric(x$C0))
        stop("Component C0 must be numeric")
    if (!( nrow(x$C0) == p && ncol(x$C0) == p))
        stop("Incompatible dimensions of matrices")
    if (!( is.numeric(x$m0) && NCOL(x$m0) == 1 && NROW(x$m0) == p ))
        stop(paste("Component m0 must be a numeric vector of length",
                   "\n equal to ncol of component FF, or a matrix with one column and",
                   "\n number of rows equal to ncol of component FF"))
    if (!(all.equal(x$C0, t(x$C0)) && all(eigen(x$C0)$values >= 0)))
        stop("C0 is not a valid variance matrix")
    if (any( c(is.na(x$m0), is.na(x$C0))))
        stop("Missing values are not allowed in components m0 and C0")
    ## extra components for time-varying dlm
    nm1 <- c("JFF", "JV", "JGG", "JW")
    nm1Ind <- match(nm1, names(x))
    if (all(is.na(nm1Ind))) {
        if (!(all.equal(x$V, t(x$V)) && all(eigen(x$V)$values >= 0)))
            stop("V is not a valid variance matrix")
        if (!(all.equal(x$W, t(x$W)) && all(eigen(x$W)$values >= 0)))
            stop("W is not a valid variance matrix")
        mod <- x[nmInd]
        class(mod) <- "dlm"
        return(mod)
    }
    x[nm1Ind[!is.na(nm1Ind)]] <-
        lapply(x[nm1Ind[!is.na(nm1Ind)]], function(x) if (!is.null(x)) as.matrix(x))
    if (!is.null(x$JFF)) {
        if (!(is.numeric(x$JFF) &&
              nrow(x$JFF) == m && ncol(x$JFF) == p))
            stop("Invalid component JFF")
        JFF <- round(x$JFF)
        if (all(JFF == 0)) JFF <- NULL
    } else
    JFF <- NULL
    if (!is.null(x$JV)) {
        if (!(is.numeric(x$JV) &&
              nrow(x$JV) == m && ncol(x$JV) == m))
            stop("Invalid component JV")
        JV <- round(x$JV)
        if (all(JV == 0)) JV <- NULL
    } else
    JV <- NULL
    if (!is.null(x$JGG)) {
        if (!(is.numeric(x$JGG) &&
              nrow(x$JGG) == p && ncol(x$JGG) == p))
            stop("Invalid component JGG")
        JGG <- round(x$JGG)
        if (all(JGG == 0)) JGG <- NULL
    } else
    JGG <- NULL
    if (!is.null(x$JW)) {
        if (!(is.numeric(x$JW) &&
              nrow(x$JW) == p && ncol(x$JW) == p))
            stop("Invalid component JW")
        JW <- round(x$JW)
        if (all(JW == 0)) JW <- NULL
    } else
    JW <- NULL
    mx <- max(c(JFF, JV, JGG, JW))
    if ( mx <= 0 ) {
        mod <- x[nmInd]
        class(mod) <- "dlm"
        return(mod)
    }
    if ( is.null(x$X) )
        stop("Component X must be provided for time-varying models")
    x$X <- as.matrix(x$X)
    if (!(is.numeric(x$X) && ncol(x$X) >= mx))
        stop("Invalid component X")
    mod <- c(x[nmInd], list(JFF=JFF, JV=JV, JGG=JGG, JW=JW, X=x$X))
    class(mod) <- "dlm"
    return(mod)
}

is.dlm <- function(obj) inherits(obj, "dlm")

as.dlm <- function(obj)
    if (is.dlm(obj)) obj else dlm(obj)

is.dlmFiltered <- function(obj) inherits(obj, "dlmFiltered")


###### Extractors and replacement functions
FF <- function(x) UseMethod("FF")
"FF<-" <- function(x, value) UseMethod("FF<-")
GG <- function(x) UseMethod("GG")
"GG<-" <- function(x, value) UseMethod("GG<-")
V <- function(x) UseMethod("V")
"V<-" <- function(x, value) UseMethod("V<-")
W <- function(x) UseMethod("W")
"W<-" <- function(x, value) UseMethod("W<-")
C0 <- function(x) UseMethod("C0")
"C0<-" <- function(x, value) UseMethod("C0<-")
m0 <- function(x) UseMethod("m0")
"m0<-" <- function(x, value) UseMethod("m0<-")

FF.dlm <- function(x)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (!is.null(x$JFF))
        warning(paste("Time varying", sQuote("F")))
    return(x$FF)
}

"FF<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("F")))
    if (!(is.numeric(value) && is.matrix(value)))
        stop(paste(deparse(substitute(value)), "is not a valid observation matrix"))
    if ( any(dim(x$FF) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    if (!is.null(x$JFF))
        warning(paste("time varying", sQuote("F"), "in", sQuote("x")))
    x$FF <- value
    return(x)
}

GG.dlm <- function(x)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (!is.null(x$JGG))
        warning(paste("Time varying", sQuote("G")))
    return(x$GG)
}

"GG<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("G")))
    if (!(is.numeric(value) && is.matrix(value)))
        stop(paste(deparse(substitute(value)), "is not a valid evolution matrix"))
    if ( any(dim(x$GG) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    if (!is.null(x$JGG))
        warning(paste("Time varying", sQuote("G"), "in", sQuote("x")))
    x$GG <- value
    return(x)
}

V.dlm <- function(x)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (!is.null(x$JV))
        warning(paste("Time varying", sQuote("V")))
    return(x$V)
}

"V<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("V")))
    if (!(is.numeric(value) && is.matrix(value) && all.equal(value, t(value))
          && all(eigen(value)$values > -.Machine$double.neg.eps)))
        stop(paste(deparse(substitute(value)), "is not a valid variance matrix"))
    if ( any(dim(x$V) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    if (!is.null(x$JV))
        warning(paste("Time varying", sQuote("V"), "in", sQuote("x")))
    x$V <- value
    return(x)
}

W.dlm <- function(x)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (!is.null(x$JW))
        warning(paste("Time varying", sQuote("W")))
    return(x$W)
}

"W<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("W")))
    if (!(is.numeric(value) && is.matrix(value) && all.equal(value, t(value))
          && all(eigen(value)$values > -.Machine$double.neg.eps)))
        stop(paste(deparse(substitute(value)), "is not a valid variance matrix"))
    if ( any(dim(x$W) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    if (!is.null(x$JW))
        warning(paste("Time varying", sQuote("W"), "in", sQuote("x")))
    x$W <- value
    return(x)
}

C0.dlm <- function(x)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    return(x$C0)
}

"C0<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("C0")))
    if (!(is.numeric(value) && is.matrix(value) && all.equal(value, t(value))
          && all(eigen(value)$values > 0)))
        stop(paste(deparse(substitute(value)), "is not a valid variance matrix"))
    if ( any(dim(x$C0) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    x$C0 <- value
    return(x)
}

m0.dlm <- function(x) {
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    return(x$m0)
}

"m0<-.dlm" <- function(x, value)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("m0")))
    if (!is.numeric(value))
        stop(paste(deparse(substitute(value)), "is not numeric"))
    if (!(length(value) == length(x$m0) && NCOL(value) == 1))
        stop(paste("wrong length/dimension of", deparse(substitute(value))))
    x$m0 <- value
    return(x)
}

### Stuff for time-varying models
JFF <- function(x) UseMethod("JFF")
"JFF<-" <- function(x, value) UseMethod("JFF<-")
JGG <- function(x) UseMethod("JGG")
"JGG<-" <- function(x, value) UseMethod("JGG<-")
JV <- function(x) UseMethod("JV")
"JV<-" <- function(x, value) UseMethod("JV<-")
JW <- function(x) UseMethod("JW")
"JW<-" <- function(x, value) UseMethod("JW<-")
X <- function(x) UseMethod("X")
"X<-" <- function(x, value) UseMethod("X<-")

JFF.dlm <- function(x) {
    return(x$JFF)
}

"JFF<-.dlm" <- function(x, value)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    value <- as.matrix(value)
    dm <- dim(value)
    value <- as.integer(value)
    dim(value) <- dm
    if (any(is.na(value)))
        stop(paste("missing values not allowed in", sQuote("JFF")))
    if ( any(dim(x$FF) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    x$JFF <- value
    return(x)
}

JGG.dlm <- function(x) {
    return(x$JGG)
}

"JGG<-.dlm" <- function(x, value)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    value <- as.matrix(value)
    dm <- dim(value)
    value <- as.integer(value)
    dim(value) <- dm
    if (any(is.na(value)))
        stop(paste("missing values not allowed in", sQuote("JGG")))
    if ( any(dim(x$GG) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    x$JGG <- value
    return(x)
}

JV.dlm <- function(x) {
    return(x$JV)
}

"JV<-.dlm" <- function(x, value)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    value <- as.matrix(value)
    dm <- dim(value)
    value <- as.integer(value)
    dim(value) <- dm
    if (any(is.na(value)))
        stop(paste("missing values not allowed in", sQuote("JV")))
    if ( any(dim(x$V) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    x$JV <- value
    return(x)
}

JW.dlm <- function(x) {
    return(x$JW)
}

"JW<-.dlm" <- function(x, value)
{
    if (!is.dlm(x))
        stop(paste(deparse(substitute(x)), "is not a dlm object"))
    value <- as.matrix(value)
    dm <- dim(value)
    value <- as.integer(value)
    dim(value) <- dm
    if (any(is.na(value)))
        stop(paste("missing values not allowed in", sQuote("JW")))
    if ( any(dim(x$W) != dim(value)) )
        stop(paste("wrong dimension of", deparse(substitute(value))))
    x$JW <- value
    return(x)
}

X.dlm <- function(x) {
    return(x$X)
}

"X<-.dlm" <- function(x, value)
{
    value <- as.matrix(value)
    dm <- dim(value)
    value <- as.numeric(value)
    dim(value) <- dm
    if (any(!is.finite(value)))
        stop(paste("finite values needed in", sQuote("X")))
    x$X <- value
    return(x)
}


###### Regression
dlmModReg <- function(X, addInt=TRUE, dV=1, dW=rep(0,NCOL(X)+addInt),
                      m0=rep(0,length(dW)), C0=1e7*diag(nrow=length(dW)))
{
    ## sanity checks
    p <- NCOL(X) + addInt
    if (!( length(dV)==1 && length(dW)==p &&
          length(m0)==p && nrow(C0)==p && ncol(C0)==p))
        stop("Inconsistent dimensions of arguments")
    X <- as.matrix(X)
    JFF <- matrix(1:ncol(X),nrow=1)
    if (addInt)
        JFF <- cbind(0,JFF)
    mod <- list(
                m0 = m0,
                C0 = C0,
                FF = matrix(1,1,p),
                V = as.matrix(dV),
                GG = diag(nrow=p),
                W = diag(x=dW, nrow=p),
                JFF = JFF,
                JV = NULL,
                JGG = NULL,
                JW = NULL,
                X = X)
    class(mod) <- "dlm"
    return(mod)
}

###### Polynomial trend
dlmModPoly <- function(order=2, dV=1,
                       dW=c(rep(0,order-1),1),
                       m0=rep(0,order), C0=1e7*diag(nrow=order))
{
    C0 <- as.matrix(C0)
    ## sanity checks
    if (!( length(dV)==1 && length(dW)==order &&
          length(m0)==order && nrow(C0)==order && ncol(C0)==order))
        stop("Inconsistent dimensions of arguments")
    GG <- diag(order)
    GG[row(GG) == col(GG)-1] <- 1
    mod <- list(
                m0 = m0,
                C0 = C0,
                FF = matrix(c(1,rep(0,order-1)),nrow=1),
                V = as.matrix(dV),
                GG = GG,
                W = diag(dW,nrow=length(dW)),
                JFF = NULL,
                JV = NULL,
                JGG = NULL,
                JW = NULL)
    class(mod) <- "dlm"
    return(mod)
}

###### Seasonal factors
dlmModSeas <- function(frequency, dV=1, dW=c(1,rep(0,frequency-2)),
                       m0=rep(0,frequency-1),
                       C0=1e7*diag(nrow=frequency-1))
{
    frequency <- as.integer(frequency)
    p <- frequency - 1
    if (!( length(dV) == 1 && length(dW) == p && length(m0) == p &&
          nrow(C0) == p && ncol(C0) == p ))
        stop("Inconsistent dimensions of arguments")
    GG <- matrix(0,p,p)
    GG[row(GG) == col(GG)+1] <- 1
    GG[1,] <- -1
    mod <- list(
                m0 = m0,
                C0 = C0,
                FF = matrix(c(1,rep(0,p-1)),nrow=1),
                V = as.matrix(dV),
                GG = GG,
                W = diag(dW,nrow=length(dW)),
                JFF = NULL,
                JV = NULL,
                JGG = NULL,
                JW = NULL)
    class(mod) <- "dlm"
    return(mod)
}

###### Fourier representation
dlmModTrig <- function(s, q, om, tau, dV=1, dW=0, m0, C0)
{
    if ( hasArg(s) ) {
        if ( hasArg(om) )
            stop("Cannot specify both 's' and 'om'")
        if ( hasArg(tau) )
            stop("Cannot specify both 's' and 'tau'")
        s <- as.integer(s)
        sHalf <- s %/% 2
        if ( hasArg(q) ) {
            q <- as.integer(q)
            if ( q > sHalf )
                stop(paste("Can use",sHalf,"components at most (",q,"were asked for )"))
            if ( q <= 0 ) stop("'q' must be positive")
        } else {
            q <- sHalf }
        om <- 2 * base::pi / s
        evenAll <- (q == sHalf) && !(s %% 2)
        if ( sHalf == 1 && evenAll ) {
            FF <- matrix(1,1,1)
            GG <- matrix(-1,1,1)
        } else {
            H1 <- diag(x=cos(om), nrow=2)
            H1[1,2] <- sin(om)
            H1[2,1] <- - H1[1,2]
            if ( q > 1 ) {
                h <- vector("list",q)
                h[[1]] <- H1
                i <- 2
                while ( i < q ) {
                    h[[i]] <- h[[i-1]] %*% H1
                    i <- i + 1
                }
                h[[q]] <- if ( evenAll ) matrix(-1,1,1) else h[[q-1]] %*% H1
                GG <- bdiag(h)
            } else {
                GG <- H1 }
        }
    } else {
        if ( !hasArg(om) )
            if ( !hasArg(tau) )
                stop("One of 's', 'om' or 'tau' must be specified")
            else om <- 2 * base::pi / tau
        if ( !hasArg(q) )
            stop("When 'om' or 'tau' is specified, 'q' must be supplied as well")
        q <- as.integer(q)
        if ( q <= 0 ) stop("'q' must be positive")
        evenAll <- FALSE
        H1 <- diag(x=cos(om), nrow=2)
        H1[1,2] <- sin(om)
        H1[2,1] <- - H1[1,2]
        if ( q > 1 ) {
            h <- vector("list",q)
            h[[1]] <- H1
            for ( i in 2:q )
                h[[i]] <- h[[i-1]] %*% H1
            GG <- bdiag(h)
        } else {
            GG <- H1 }
    }
    if ( hasArg(m0) ) {
        if ( length(m0) != nrow(GG) ) stop("Wrong length of 'm0'")
    } else {
        m0 <- rep(0,nrow(GG)) }
    if ( hasArg(C0) ) {
        if ( (nrow(C0) != nrow(GG)) || (ncol(C0) != nrow(GG)) )
            stop("Wrong dimension of 'C0'")
    } else {
        C0 <- diag(x=1e7, nrow=nrow(GG)) }
    mod <- list(
                m0 = m0,
                C0 = C0,
                FF = matrix(if (evenAll) c(rep(c(1,0),q-1),1)
                else c(1,0),1,nrow(GG)),
                V = matrix(dV,1,1),
                GG = GG,
                W = diag(x=dW,nrow=nrow(GG)),
                JFF = NULL,
                JV = NULL,
                JGG = NULL,
                JW = NULL)
    class(mod) <- "dlm"
    return(mod)
}

###### ARMA
dlmModARMA <- function(ar=NULL, ma=NULL, sigma2=1, dV, m0, C0)
{
    if (is.matrix(sigma2) && (m <- nrow(sigma2)) > 1)
    { ## multivariate
        if (ncol(sigma2) != m)
            stop("sigma2 must be a square matrix")
        r <- max(p <- length(ar), (q <- length(ma)) + 1)
        k <- m * r
        if (hasArg("dV")) {
            if ( length(dV) != m )
                stop("Incompatible dimensions of arguments")
        }
        else
            dV <- rep(0,m)
        if (hasArg("m0")) {
            if ( length(m0) != k )
                stop("Incompatible dimensions of arguments")
        }
        else
            m0 <- rep(0,k)
        if (hasArg("C0")) {
            if ( !( nrow(C0) == k && ncol(C0) == k ))
                stop("Incompatible dimensions of arguments")
        }
        else
            C0 <- 1e7*diag(nrow=k)
        FF <- matrix(0,m,k)
        FF[row(FF) == col(FF)] <- 1
        GG <- matrix(0,k,k)
        for (i in seq(length.out=p))  GG[ (1 + (i-1) * m):(i * m), 1:m ] <- ar[[i]]
        GG[row(GG) == col(GG) - m] <- 1
        R <- matrix(0,k,m)
        R[row(R) == col(R)] <- 1
        for (i in seq(length.out=q)) R[ (1 + i * m):((i+1) * m), ] <- ma[[i]]
        mod <- list(
                    m0 = m0,
                    C0 = C0,
                    FF = FF,
                    V = diag(dV,nrow=m),
                    GG = GG,
                    W = R %*% sigma2 %*% t(R),
                    JFF = NULL,
                    JV = NULL,
                    JGG = NULL,
                    JW = NULL)
    }
    else
    { ## univariate
        r <- max(p <- length(ar), (q <- length(ma)) + 1)
        if (hasArg("dV")) {
            if ( length(dV) != 1 )
                stop("Incompatible dimensions of arguments")
        }
        else
            dV <- 0
        if (hasArg("m0")) {
            if ( length(m0) != r )
                stop("Incompatible dimensions of arguments")
        }
        else
            m0 <- rep(0,r)
        if (hasArg("C0")) {
            if ( !( nrow(C0) == r && ncol(C0) == r ))
                stop("Incompatible dimensions of arguments")
        }
        else
            C0 <- 1e7*diag(nrow=r)
        GG <- matrix(0,r,r)
        if ( p > 0 ) GG[ 1:p, 1 ] <- ar
        GG[row(GG) == col(GG) - 1] <- 1
        R <- rep(0,r)
        R[1] <- 1
        if ( q > 0 ) R[ 1 + 1:q ] <- ma
        mod <- list(
                    m0 = m0,
                    C0 = C0,
                    FF = matrix(c(1,rep(0,r-1)),nrow=1),
                    V = matrix(dV),
                    GG = GG,
                    W = sigma2 * crossprod(matrix(R,nrow=1)),
                    JFF = NULL,
                    JV = NULL,
                    JGG = NULL,
                    JW = NULL)
    }
    class(mod) <- "dlm"
    return(mod)
}


###### Addition for "dlm" objects
"+.dlm" <- function(mod1, mod2)
{
    if ( (m <- nrow(mod1$FF)) != nrow(mod2$FF) )
        stop("Incompatible models")
    if ( (x1 <- !is.null(mod1$X)) | (x2 <- !is.null(mod2$X)) ) {
        if ( x1 && x2 && ( nrow(mod1$X) != nrow(mod2$X) ) )
            stop("Number of rows of mod1$X and mod2$X must match")
        plus <- function(u,v) ifelse( u==0, 0, u+v )
        p1 <- ncol(mod1$FF)
        p2 <- ncol(mod2$FF)
        r <- if ( x1 ) ncol(mod1$X) else 0
        if ( (one <- !is.null(mod1$JFF)) | (two <- !is.null(mod2$JFF)) ) {
            JFF1 <- if ( one ) mod1$JFF else matrix(0,m,p1)
            JFF2 <- if ( two ) plus(mod2$JFF,r) else matrix(0,m,p2)
            JFF <- cbind(JFF1, JFF2)
        } else
        JFF <- NULL
        if ( (one <- !is.null(mod1$JGG)) | (two <- !is.null(mod2$JGG)) ) {
            JGG1 <- if ( one ) mod1$JGG else matrix(0,p1,p1)
            JGG2 <- if ( two ) plus(mod2$JGG,r) else matrix(0,p2,p2)
            JGG <- bdiag(list(JGG1, JGG2))
        } else
        JGG <- NULL
        if ( (one <- !is.null(mod1$JW)) | (two <- !is.null(mod2$JW)) ) {
            JW1 <- if ( one ) mod1$JW else matrix(0,p1,p1)
            JW2 <- if ( two ) plus(mod2$JW,r) else matrix(0,p2,p2)
            JW <- bdiag(list(JW1, JW2))
        } else
        JW <- NULL
        if ( (one <- !is.null(mod1$JV)) | (two <- !is.null(mod2$JV)) ) {
            if ( one && two )
                stop("mod1$V and mod2$V cannot be both time-varying")
            if ( one ) {
                if ( any(mod2$V != 0) ) {
                    mod2$V[] <- 0
                    warning("the value of mod2$V has been discarded")
                }
                JV <- mod1$JV
            }
            if ( two ) {
                if ( any(mod1$V != 0) ) {
                    mod1$V[] <- 0
                    warning("the value of mod1$V has been discarded")
                }
                JV <- plus(mod2$JV,r)
            }
        } else
        JV <- NULL
        mod <- list(
                    m0 = c(mod1$m0, mod2$m0),
                    C0 = bdiag(list(mod1$C0, mod2$C0)),
                    FF = cbind(mod1$FF, mod2$FF),
                    V = mod1$V + mod2$V,
                    GG = bdiag(list(mod1$GG, mod2$GG)),
                    W = bdiag(list(mod1$W, mod2$W)),
                    JFF = JFF,
                    JV = JV,
                    JGG = JGG,
                    JW = JW,
                    X = switch( x1 + 2 * x2,
                    mod1$X,
                    mod2$X,
                    cbind( mod1$X, mod2$X) ))
    } else
    mod <- list(
                m0 = c(mod1$m0, mod2$m0),
                C0 = bdiag(list(mod1$C0, mod2$C0)),
                FF = cbind(mod1$FF, mod2$FF),
                V = mod1$V + mod2$V,
                GG = bdiag(list(mod1$GG, mod2$GG)),
                W = bdiag(list(mod1$W, mod2$W)),
                JFF = NULL,
                JV = NULL,
                JGG = NULL,
                JW = NULL)
    class(mod) <- "dlm"
    return(mod)
}


###### Outer sum of DLMs
"dlmSum" <- function(...)
{
    if ( (narg <- nargs()) == 1  ) {
        if ( is.dlm(...) )
            return(...)
        else
            stop("Argument must be a \"dlm\" object")
    }
    else
    {
        args <- list(...)
        if ( !all(sapply(args, is.dlm)) )
            stop("Arguments must be \"dlm\" objects")
        if ( narg > 2 )
            return( dlmSum(args[[1]], do.call("dlmSum", args[-1])))
        ## now we are in the case narg == 2
        if ((x1 <- !is.null(args[[1]]$X)) | (x2 <- !is.null(args[[2]]$X)))
            stop("Sum of dlm's is only implemented for constant models")
        else
            mod <- list(m0 = c(args[[1]]$m0, args[[2]]$m0),
                        C0 = bdiag(list(args[[1]]$C0, args[[2]]$C0)),
                        FF = bdiag(list(args[[1]]$FF, args[[2]]$FF)),
                        V = bdiag(list(args[[1]]$V, args[[2]]$V)),
                        GG = bdiag(list(args[[1]]$GG, args[[2]]$GG)),
                        W = bdiag(list(args[[1]]$W, args[[2]]$W)),
                        JFF = NULL, JV = NULL, JGG = NULL, JW = NULL)
        class(mod) <- "dlm"
        return(mod)
    }
}

"%+%" <- function(x, y)
    dlmSum(x, y)


###### Print method for "dlm" objects
print.dlm <- function(x, ...) {
    nmRef <- c("FF","V","GG","W","JFF","JV","JGG","JW","X","m0","C0")
    nm <- names(x)
    what <- !sapply(x, is.null)
    ind <- match(nmRef,nm)
    ind <- ind[!is.na(ind)]
    good <- nm[ind][what[ind]]
    if ( !is.na(match("X",good)) && nrow(x$X) > 2 ) {
        x$X <- rbind(formatC(x$X[1:2,,drop=FALSE],digits=4,format="fg"),
                       c("...",rep("",ncol(x$X)-1)))
    }
    print(x[match(good,nm)], quote=FALSE)
    invisible(x)
}

###### Negative loglikelihood
dlmLL <- function(y, mod, debug=FALSE)
{
    ## calculations based on singular value decomposition
    ## Note: V must be nonsingular
    ## The C code relies on the order of the elements in 'mod'
    ## mod = list(m0, C0, FF, V, GG, W)
     if (!debug) {
        storage.mode(y) <- "double"
        matchnames <- match(c("m0", "C0", "FF", "V", "GG", "W"), names(mod))
        for (i in matchnames)
            storage.mode(mod[[i]]) <- "double"
        ## define flags for time-varying components
        if (is.null(mod$JFF))
            tvFF <- FALSE
        else {
            tvFF <- TRUE
            nz <- mod$JFF != 0
            mod$JFF <- cbind(row(mod$JFF)[nz], col(mod$JFF)[nz], mod$JFF[nz]) - 1
            storage.mode(mod$JFF) <- "integer"
        }
        if (is.null(mod$JV))
            tvV <- FALSE
        else {
            tvV <- TRUE
            nz <- mod$JV != 0
            mod$JV <- cbind(row(mod$JV)[nz], col(mod$JV)[nz], mod$JV[nz]) - 1
            storage.mode(mod$JV) <- "integer"
        }
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz]) - 1
            storage.mode(mod$JGG) <- "integer"
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz]) - 1
            storage.mode(mod$JW) <- "integer"
        }
        if (any(c(tvFF,tvV,tvGG,tvW))) {
            mod <- mod[match(c("m0", "C0", "FF", "V", "GG", "W",
                               "JFF", "JV", "JGG", "JW", "X"), names(mod))]
            return(.Call(C_dlmLL, y, mod, tvFF, tvV, tvGG, tvW, PACKAGE="dlm"))
        } else {
            mod <- mod[match(c("m0", "C0", "FF", "V", "GG", "W"), names(mod))]
            return(.Call(C_dlmLL0, y, mod, PACKAGE="dlm"))
        }
    }
    else {
        eps <- .Machine$double.eps^.3
        y <- as.matrix(y)
        n <- nrow(y)
        ll <- 0 # negative loglikelihood
        ## define flags for time-varying components
        if (is.null(mod$JFF))
            tvFF <- FALSE
        else {
            tvFF <- TRUE
            nz <- mod$JFF != 0
            mod$JFF <- cbind(row(mod$JFF)[nz], col(mod$JFF)[nz], mod$JFF[nz])
        }
        if (is.null(mod$JV))
            tvV <- FALSE
        else {
            tvV <- TRUE
            nz <- mod$JV != 0
            mod$JV <- cbind(row(mod$JV)[nz], col(mod$JV)[nz], mod$JV[nz])
        }
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz])
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz])
        }
        tvFV <- tvFF || tvV
        ## preliminary calculations, if possible (non time-varying case)
        if ( !tvV ) {
            tmp <- La.svd(mod$V,nu=0)
            Dv <- sqrt(tmp$d)
            if (any(Dv < eps)) {
              Dv <- pmax(Dv, eps)
              warning("a numerically singular 'V' has been slightly perturbed to make it nonsingular")
            }
            Dv.inv <- 1/Dv
            sqrtVinv <- Dv.inv * tmp$vt # t() %*% () = V^{-1}
            sqrtV <- Dv * tmp$vt # t()%*%() = V
            if ( !tvFF )
                tF.Vinv <- t(mod$FF) %*% crossprod(sqrtVinv)
        }
        if ( !tvW ) {
            svdW <- La.svd(mod$W,nu=0)
            sqrtW <- sqrt(svdW$d) * svdW$vt # t()%*%() = W
        }
        tmp <- La.svd(mod$C0,nu=0)
        Ux <- t(tmp$vt); Dx <- sqrt(tmp$d)
        for (i in seq(length = n)) {
            ## set time-varying matrices
            if ( tvFF )
                mod$FF[mod$JFF[,-3,drop=FALSE]] <- mod$X[i,mod$JFF[,3]]
            if ( tvV ) {
                mod$V[mod$JV[,-3,drop=FALSE]] <- mod$X[i,mod$JV[,3]]
                tmp <- La.svd(mod$V,nu=0)
                Dv <- sqrt(tmp$d)
                Dv.inv <- 1/Dv; Dv.inv[abs(Dv.inv)==Inf] <- 0
                sqrtVinv <- Dv.inv * tmp$vt
                sqrtV <- sqrt(tmp$d) * tmp$vt # t()%*%() = V
            }
            if ( tvGG )
                mod$GG[mod$JGG[,-3,drop=FALSE]] <- mod$X[i,mod$JGG[,3]]
            if ( tvW ) {
                mod$W[mod$JW[,-3,drop=FALSE]] <- mod$X[i,mod$JW[,3]]
                svdW <- La.svd(mod$W,nu=0)
                sqrtW <- sqrt(svdW$d) * svdW$vt # t()%*%() = W
            }
            if ( tvFV )
                tF.Vinv <- t(mod$FF) %*% crossprod(sqrtVinv)

            if (!any(whereNA <- is.na(y[i, ]))) { ## No missing values

                ## prior
                a <- mod$GG %*% mod$m0
                tmp <- La.svd(rbind( Dx*t(mod$GG%*%Ux), sqrtW ), nu=0)
                Ux.prior <- t(tmp$vt)
                Dx.prior <- tmp$d
                ## one-step forecast
                f <- mod$FF %*% a
                tmp <- La.svd(rbind( Dx.prior*t(mod$FF%*%Ux.prior), sqrtV ), nu=0)
                Uy <- t(tmp$vt)
                Dy <- tmp$d
                ## posterior
                D.inv <- 1/Dx.prior
                D.inv[abs(D.inv)==Inf] <- 0
                tmp <- La.svd(rbind(sqrtVinv%*%mod$FF%*%Ux.prior,
                                    diag(x=D.inv,nrow=length(D.inv))), nu=0)
                Ux <- Ux.prior %*% t(tmp$vt)
                Dx <- 1/tmp$d
                Dx[abs(Dx)==Inf] <- 0
                e <- as.matrix(y[i,]-f)
                mod$m0 <- a + crossprod(Dx*t(Ux)) %*%
                    tF.Vinv %*% e
                ## update scaled negative loglikelihood
                ll <- ll + 2*sum(log(Dy)) + crossprod(crossprod(Uy,e)/Dy)

            } else {
                if (all(whereNA)) { ## All components missing

                    ## prior & posterior
                    mod$m0 <- mod$GG %*% mod$m0
                    tmp <- La.svd(rbind( Dx*t(mod$GG%*%Ux), sqrtW ), nu=0)
                    Ux <- t(tmp$vt)
                    Dx <- tmp$d

                } else { ## Some components missing

                    good <- !whereNA
                    tmp <- La.svd(mod$V[good, good], nu=0)
                    Dv <- sqrt(tmp$d)
                    Dv.inv <- 1/Dv; Dv.inv[abs(Dv.inv)==Inf] <- 0
                    sqrtVinvTMP <- Dv.inv * tmp$vt
                    tF.VinvTMP <- t(mod$FF[good,,drop=FALSE]) %*% crossprod(sqrtVinvTMP)
                    sqrtVTMP <- Dv * tmp$vt
                    ## prior
                    a <- mod$GG %*% mod$m0
                    tmp <- La.svd(rbind( Dx*t(mod$GG%*%Ux), sqrtW ), nu=0)
                    Ux.prior <- t(tmp$vt)
                    Dx.prior <- tmp$d
                    ## one-step forecast
                    f <- mod$FF[good,,drop=FALSE] %*% a
                    tmp <- La.svd(rbind( Dx.prior*t(mod$FF[good,,drop=FALSE]%*%Ux.prior), sqrtVTMP ), nu=0)
                    Uy <- t(tmp$vt)
                    Dy <- tmp$d
                    ## posterior
                    D.inv <- 1/Dx.prior
                    D.inv[abs(D.inv)==Inf] <- 0
                    tmp <- La.svd(rbind(sqrtVinvTMP%*%mod$FF[good,,drop=FALSE]%*%Ux.prior,
                                        diag(x=D.inv,nrow=length(D.inv))), nu=0)
                    Ux <- Ux.prior %*% t(tmp$vt)
                    Dx <- 1/tmp$d
                    Dx[abs(Dx)==Inf] <- 0
                    e <- as.matrix(y[i,good]-f)
                    mod$m0 <- a + crossprod(Dx*t(Ux)) %*%
                        tF.VinvTMP %*% e
                    ## update scaled negative loglikelihood
                    ll <- ll + 2*sum(log(Dy)) + crossprod(crossprod(Uy,e)/Dy)
                }
            }
        }
    }
    return(drop(ll)*0.5) # negative loglikelihood
}


"dlmMLE" <- function(y, parm, build, method = "L-BFGS-B",
                     ..., debug = FALSE)
{
    logLik <- function(parm, ...)
    {
        mod <- build(parm, ...)
        return(dlmLL(y=y, mod=mod, debug=debug))
    }
    out <- optim(parm, logLik, method=method, ...)
    return(out)
}


dlmFilter <- function(y, mod, debug = FALSE, simplify = FALSE)
{
    ## Note: V must be nonsingular. It will be forced to be so,
    ## with a warning, otherwise (time-invariant case only).
    eps <- .Machine$double.eps^.3
    mod1 <- mod
    yAttr <- attributes(y)
    ytsp <- tsp(y)
    y <- as.matrix(y)
    timeNames <- dimnames(y)[[1]]
    stateNames <- names(mod$m0)
    if (!debug) {
        storage.mode(y) <- "double"
        matchnames <- match(c("m0", "C0", "FF", "V", "GG", "W"), names(mod))
        for (i in matchnames)
            storage.mode(mod[[i]]) <- "double"
        if ("X" %in% names(mod))
          storage.mode(mod$X) <- "double"
        ## define flags for time-varying components
        if (is.null(mod$JFF))
            tvFF <- FALSE
        else {
            tvFF <- TRUE
            nz <- mod$JFF != 0
            mod$JFF <- cbind(row(mod$JFF)[nz], col(mod$JFF)[nz], mod$JFF[nz]) - 1
            storage.mode(mod$JFF) <- "integer"
        }
        if (is.null(mod$JV))
            tvV <- FALSE
        else {
            tvV <- TRUE
            nz <- mod$JV != 0
            mod$JV <- cbind(row(mod$JV)[nz], col(mod$JV)[nz], mod$JV[nz]) - 1
            storage.mode(mod$JV) <- "integer"
        }
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz]) - 1
            storage.mode(mod$JGG) <- "integer"
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz]) - 1
            storage.mode(mod$JW) <- "integer"
        }
        if (any(c(tvFF,tvV,tvGG,tvW))) {
            mod <- mod[match(c("m0", "C0", "FF", "V", "GG", "W",
                               "JFF", "JV", "JGG", "JW", "X"), names(mod))]
            ans <- .Call(C_dlmFilter, y, mod, tvFF, tvV, tvGG, tvW, PACKAGE = "dlm")
        } else {
            mod <- mod[match(c("m0", "C0", "FF", "V", "GG", "W"), names(mod))]
            ans <- .Call(C_dlmFilter0, y, mod, PACKAGE = "dlm")
        }
        names(ans) <- c("m", "U.C", "D.C", "a", "U.R", "D.R", "f")
    }
    else {
        ## define flags for time-varying components
        if (is.null(mod$JFF))
            tvFF <- FALSE
        else {
            tvFF <- TRUE
            nz <- mod$JFF != 0
            mod$JFF <- cbind(row(mod$JFF)[nz], col(mod$JFF)[nz], mod$JFF[nz])
        }
        if (is.null(mod$JV))
            tvV <- FALSE
        else {
            tvV <- TRUE
            nz <- mod$JV != 0
            mod$JV <- cbind(row(mod$JV)[nz], col(mod$JV)[nz], mod$JV[nz])
        }
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz])
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz])
        }
        tvFV <- tvFF || tvV
        m <- rbind(mod$m0,matrix(0,nrow=nrow(y),ncol=length(mod$m0))) # filtered values
        a <- matrix(0,nrow=nrow(y),ncol=length(mod$m0))
        f <- matrix(0,nrow=nrow(y),ncol=ncol(y))
        U.C <- vector(1+nrow(y),mode="list")
        D.C <- matrix(0,1+nrow(y),length(mod$m0))
        U.R <- vector(nrow(y),mode="list")
        D.R <- matrix(0,nrow(y),length(mod$m0))
        ## preliminary calculations, if possible (non time-varying case)
        if ( !tvV ) {
            tmp <- La.svd(mod$V,nu=0)
            Uv <- t(tmp$vt); Dv <- sqrt(tmp$d)
            if (any(Dv < eps)) {
              Dv <- pmax(Dv, eps)
              warning("a numerically singular 'V' has been slightly perturbed to make it nonsingular")
            }
            Dv.inv <- 1/Dv
            sqrtVinv <- Dv.inv * tmp$vt # t() %*% () = V^{-1}
            if ( !tvFF )
                tF.Vinv <- t(mod$FF) %*% crossprod(sqrtVinv)
        }
        if ( !tvW ) {
            svdW <- La.svd(mod$W,nu=0)
            sqrtW <- sqrt(svdW$d) * svdW$vt # t()%*%() = W
        }
        tmp <- La.svd(mod$C0,nu=0)
        U.C[[1]] <- t(tmp$vt)
        D.C[1,] <- sqrt(tmp$d)
        for (i in seq(length=nrow(y))) {
            ## set time-varying matrices
            if ( tvFF )
                mod$FF[mod$JFF[,-3,drop=FALSE]] <- mod$X[i,mod$JFF[,3]]
            if ( tvV ) {
                mod$V[mod$JV[,-3,drop=FALSE]] <- mod$X[i,mod$JV[,3]]
                tmp <- La.svd(mod$V,nu=0)
                Uv <- t(tmp$vt); Dv <- sqrt(tmp$d)
                Dv.inv <- 1/Dv; Dv.inv[abs(Dv.inv)==Inf] <- 0
                sqrtVinv <- Dv.inv * tmp$vt
            }
            if ( tvGG )
                mod$GG[mod$JGG[,-3,drop=FALSE]] <- mod$X[i,mod$JGG[,3]]
            if ( tvW ) {
                mod$W[mod$JW[,-3,drop=FALSE]] <- mod$X[i,mod$JW[,3]]
                svdW <- La.svd(mod$W,nu=0)
                sqrtW <- sqrt(svdW$d) * svdW$vt # t()%*%() = W
            }
            if ( tvFV )
                tF.Vinv <- t(mod$FF) %*% crossprod(sqrtVinv)

            if (!any(whereNA <- is.na(y[i, ]))) { ## No missing values

                ## prior
                a[i,] <- mod$GG %*% m[i,]
                tmp <- La.svd(rbind( D.C[i,]*t(mod$GG%*%U.C[[i]]), sqrtW ), nu=0)
                U.R[[i]] <- t(tmp$vt)
                D.R[i,] <- tmp$d
                ## one-step forecast
                f[i,] <- mod$FF %*% a[i,]
                ## posterior
                D.Rinv <- 1/D.R[i,]
                D.Rinv[abs(D.Rinv)==Inf] <- 0
                tmp <- La.svd(rbind(sqrtVinv %*% mod$FF %*% U.R[[i]],
                                    diag(x=D.Rinv,nrow=length(D.Rinv))), nu=0)
                U.C[[i+1]] <- U.R[[i]] %*% t(tmp$vt)
                foo <- 1/tmp$d; foo[abs(foo)==Inf] <- 0
                D.C[i+1,] <- foo
                m[i+1,] <- a[i,] + crossprod(D.C[i+1,]*t(U.C[[i+1]])) %*%
                    tF.Vinv %*% as.matrix(y[i,]-f[i,])

            } else {
                if (all(whereNA)) { ## All components missing

                    ## prior & posterior
                    m[i+1,] <- a[i,] <- mod$GG %*% m[i,]
                    tmp <- La.svd(rbind( D.C[i,]*t(mod$GG%*%U.C[[i]]), sqrtW ), nu=0)
                    U.C[[i+1]] <- U.R[[i]] <- t(tmp$vt)
                    D.C[i+1,] <- D.R[i,] <- tmp$d
                    ## one-step forecast
                    f[i,] <- mod$FF %*% a[i,]

                } else { ## Some components missing

                    good <- !whereNA
                    tmp <- La.svd(mod$V[good, good], nu=0)
                    Dv <- sqrt(tmp$d)
                    Dv.inv <- 1/Dv; Dv.inv[abs(Dv.inv)==Inf] <- 0
                    sqrtVinvTMP <- Dv.inv * tmp$vt
                    tF.VinvTMP <- t(mod$FF[good,,drop=FALSE]) %*% crossprod(sqrtVinvTMP)
                    ## prior
                    a[i,] <- mod$GG %*% m[i,]
                    tmp <- La.svd(rbind( D.C[i,]*t(mod$GG%*%U.C[[i]]), sqrtW ), nu=0)
                    U.R[[i]] <- t(tmp$vt)
                    D.R[i,] <- tmp$d
                    ## one-step forecast
                    f[i,] <- mod$FF %*% a[i,]
                    ## posterior
                    D.Rinv <- 1/D.R[i,]
                    D.Rinv[abs(D.Rinv)==Inf] <- 0
                    tmp <- La.svd(rbind(sqrtVinvTMP %*% mod$FF[good,,drop=FALSE] %*% U.R[[i]],
                                        diag(x=D.Rinv,nrow=length(D.Rinv))), nu=0)
                    U.C[[i+1]] <- U.R[[i]] %*% t(tmp$vt)
                    foo <- 1/tmp$d; foo[abs(foo)==Inf] <- 0
                    D.C[i+1,] <- foo
                    m[i+1,] <- a[i,] + crossprod(D.C[i+1,]*t(U.C[[i+1]])) %*%
                        tF.VinvTMP %*% as.matrix(y[i,good]-f[i,good])
                }
            }
        }

        ans <- list(m=m,U.C=U.C,D.C=D.C,a=a,U.R=U.R,D.R=D.R,f=f)
    }

    ans$m <- drop(ans$m); ans$a <- drop(ans$a); ans$f <- drop(ans$f)
    attributes(ans$f) <- yAttr
    if (!is.null(ytsp)) {
        tsp(ans$a) <- ytsp
        tsp(ans$m) <- c(ytsp[1] - 1/ytsp[3], ytsp[2:3])
        class(ans$a) <- class(ans$m) <- if (length(mod$m0) > 1) c("mts","ts") else "ts"
    }
    if (!(is.null(timeNames) && is.null(stateNames)))
        if (is.matrix(ans$a))
        {
            dimnames(ans$a) <- list(timeNames, stateNames)
            dimnames(ans$m) <- list(if(is.null(timeNames)) NULL else c("",timeNames),
                                    stateNames)
        }
        else
            if (!is.null(timeNames))
            {
                names(ans$a) <- timeNames
                names(ans$m) <- c("", timeNames)
            }
    if (simplify)
        ans <- c(mod=list(mod1), ans)
    else
    {
        attributes(y) <- yAttr
        ans <- c(y=list(y), mod=list(mod1), ans)
    }
    class(ans) <- "dlmFiltered"
    return(ans)
}

dlmSmooth <- function(y, ...)
    UseMethod("dlmSmooth", y)

dlmSmooth.default <- function(y, mod, ...)
{
    dlmSmooth(dlmFilter(y, mod, ...), ...)
}

dlmSmooth.dlmFiltered <- function(y, ..., debug = FALSE)
{
    big <- 1 / sqrt(.Machine$double.eps)
    mod <- c(y[match(c("m", "U.C", "D.C", "a", "U.R", "D.R"),names(y))],
             y$mod[match(c("GG", "W", "JGG", "JW", "X"), names(y$mod))])
    mAttr <- attributes(mod$m)
    mod$m <- as.matrix(mod$m)
    mod$a <- as.matrix(mod$a)
    if (!debug) {
        ## define flags for time-varying components
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz]) - 1
            storage.mode(mod$JGG) <- "integer"
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz]) - 1
            storage.mode(mod$JW) <- "integer"
        }
        if (tvGG || tvW) {
            ans <- .Call(C_dlmSmooth, mod, tvGG, tvW, big, PACKAGE="dlm")
            names(ans) <- c("s", "U.S", "D.S")
        }
        else {
            ans <- .Call(C_dlmSmooth0, mod, big, PACKAGE="dlm")
            names(ans) <- c("s", "U.S", "D.S")
        }
        } else {
        if (is.null(mod$JGG))
            tvGG <- FALSE
        else {
            tvGG <- TRUE
            nz <- mod$JGG != 0
            mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz])
        }
        if (is.null(mod$JW))
            tvW <- FALSE
        else {
            tvW <- TRUE
            nz <- mod$JW != 0
            mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz])
        }
        n <- length(mod$U.R) # number of obs
        p <- NCOL(mod$m) # dimension of state vector
        s <- rbind(matrix(0,n,p), mod$m[n+1,])
        U.S <- vector("list", length=n+1)
        U.S[[n+1]] <- mod$U.C[[n+1]]
        D.S <- rbind(matrix(0,n,p), mod$D.C[n+1,])
        ## preliminary calculations, if possible (time-invariant case)
        if ( !tvW ) {
            tmp <- La.svd(mod$W,nu=0)
            Dw <- sqrt(tmp$d)
            Dw.inv <- pmin(1/Dw, big)
            sqrtWinv <- Dw.inv * tmp$vt # t()%*%() = W^(-1)
        }
        if (n > 0)
            for (i in n:1)
            {
                ## set relevant time-varying matrices
                if ( tvGG )
                    mod$GG[mod$JGG[,-3,drop=FALSE]] <- mod$X[i,mod$JGG[,3]]
                if ( tvW ) {
                    mod$W[mod$JW[,-3,drop=FALSE]] <- mod$X[i,mod$JW[,3]]
                    tmp <- La.svd(mod$W,nu=0)
                    Dw <- sqrt(tmp$d)
                    Dw.inv <- pmin(1/Dw, big)
                    sqrtWinv <- Dw.inv * tmp$vt # t()%*%() = W^(-1)
                }
                Dinv <- 1/mod$D.R[i,]; Dinv[abs(Dinv)==Inf] <- 0
                H <- crossprod(mod$D.C[i,]*t(mod$U.C[[i]])) %*%
                    t(mod$GG) %*% crossprod(Dinv*t(mod$U.R[[i]]))
                Dinv <- 1/mod$D.C[i,]
                Dinv[abs(Dinv)==Inf] <- 0
                tmp <- La.svd(rbind( sqrtWinv%*%mod$GG, Dinv*t(mod$U.C[[i]])), nu=0)
                Dinv <- 1/tmp$d
                Dinv[abs(Dinv)==Inf] <- 0
                tmp <- La.svd(rbind(Dinv*tmp$vt, D.S[i+1,]*t(H%*%U.S[[i+1]])))
                U.S[[i]] <- t(tmp$vt)
                D.S[i,] <- tmp$d
                s[i,] <- mod$m[i,] + H %*% (s[i+1,]-mod$a[i,])
            }
        ans <- list(s=s, U.S=U.S, D.S=D.S)
    }
    attributes(ans$s) <- mAttr
                                        #     if (!is.null(tsp(mod$m))) {
#         tsp(ans$s) <- tsp(mod$m)
#         class(ans$s) <- if (NCOL(ans$s) > 1) c("mts","ts") else "ts"
#     }
#     dimnames(ans$s) <- dimnames(mod$m)
    return(ans)
}


dlmBSample <- function(modFilt)
{
    eps <- .Machine$double.eps^.4
    mod <- c(modFilt[match(c("m", "U.C", "D.C", "a", "U.R", "D.R"),names(modFilt))],
             modFilt$mod[match(c("GG", "W", "JGG", "JW", "X"), names(modFilt$mod))])
    n <- length(mod$U.R) # number of obs
    p <- NCOL(mod$m) # dimension of state vector
    mtsp <- tsp(mod$m)
    if (p==1) {
        dim(mod$m) <- c(n+1, 1)
        dim(mod$a) <- c(n, 1)
    }
    if (is.null(mod$JGG))
        tvGG <- FALSE
    else {
        tvGG <- TRUE
        nz <- mod$JGG != 0
        mod$JGG <- cbind(row(mod$JGG)[nz], col(mod$JGG)[nz], mod$JGG[nz])
    }
    if (is.null(mod$JW))
        tvW <- FALSE
    else {
        tvW <- TRUE
        nz <- mod$JW != 0
        mod$JW <- cbind(row(mod$JW)[nz], col(mod$JW)[nz], mod$JW[nz])
    }
    tvGW <- tvGG || tvW
    theta <- matrix(0,n+1,p)
    ## preliminary calculations, if possible (time-invariant case)
    if ( !tvW ) {
        tmp <- La.svd(mod$W,nu=0)
        Dw <- sqrt(tmp$d)
        Dw <- pmax(Dw, eps)
        Dw.inv <- 1/Dw
        sqrtWinv <- Dw.inv * tmp$vt # t()%*%() = W^(-1)
        if ( !tvGG ) tG.Winv <- t(mod$GG) %*% crossprod(sqrtWinv)
    }
    ## generate last theta
    theta[n+1,] <- mod$m[n+1,] + mod$U.C[[n+1]] %*% matrix(mod$D.C[n+1,]*rnorm(p))
    ## generate all the other theta's
    for (i in (n:1))
    {
        if ( tvGG )
            mod$GG[mod$JGG[,-3,drop=FALSE]] <- mod$X[i,mod$JGG[,3]]
        if ( tvW ) {
            mod$W[mod$JW[,-3,drop=FALSE]] <- mod$X[i,mod$JW[,3]]
            tmp <- La.svd(mod$W,nu=0)
            Dw <- sqrt(tmp$d)
            Dw <- pmax(Dw, eps)
            Dw.inv <- 1/Dw
            sqrtWinv <- Dw.inv * tmp$vt # t()%*%() = W^(-1)
        }
        if ( tvGW )
            tG.Winv <- t(mod$GG) %*% crossprod(sqrtWinv)
        D.inv <- 1/mod$D.C[i,]; D.inv[abs(D.inv)==Inf] <- 0
        tmp <- La.svd(rbind(sqrtWinv %*% mod$GG %*% mod$U.C[[i]],
                            diag(x=D.inv,nrow=length(D.inv))), nu=0)
        U.H <- mod$U.C[[i]] %*% t(tmp$vt)
        D.H <- 1/tmp$d; D.H[abs(D.H)==Inf] <- 0
        h <- mod$m[i,] + crossprod(D.H*t(U.H)) %*%
            tG.Winv %*% (t(theta[i+1,,drop=F])-mod$a[i,])
        theta[i,] <- h + U.H %*% matrix(D.H*rnorm(p))
    }
    if (!is.null(mtsp))
    {
        theta <- drop(theta)
        tsp(theta) <- mtsp
        class(theta) <- if (p > 1) c("mts","ts") else "ts"
    }
    return(theta=theta)
}


rwishart <- function(df, p = nrow(SqrtSigma), Sigma,
                     SqrtSigma = diag(p))
{
    ## generate a Wishart-distributed matrix - from S-news, due to B. Venables
    ## note: Sigma = crossprod(SqrtSigma), and df must be integer
    if (!missing(Sigma))
    {
        ## compute SqrtSigma
        tmp <- svd(Sigma)
        SqrtSigma <- sqrt(tmp$d) * t(tmp$u)
    }
    if((Ident <- missing(SqrtSigma)) && missing(p))
        stop("either p, Sigma or SqrtSigma must be specified")
    Z <- matrix(0, p, p)
    diag(Z) <- sqrt(rchisq(p, df:(df-p+1)))
    if(p > 1)
    {
        pseq <- 1:(p-1)
        Z[rep(p*pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p*(p-1)/2)
    }
    if(Ident)
        crossprod(Z)
    else
        crossprod(Z %*% SqrtSigma)
}


dlmSvd2var <- function(u, d)
{
    if (is.matrix(u)) {
        if ( nrow(u) != length(d) || ncol(u) != (n <- length(d)) )
            stop("inconsistent dimensions")
        return(tcrossprod(rep(d, each = n) * u))
    }
    if (is.list(u)) {
        if ( length(u) != NROW(d) )
            stop("length of 'u' must be equal to the number of rows of 'd'")
        n <- NCOL(d)
        return(lapply(seq(along = u), function(i)
                      tcrossprod(rep(d[i,], each = n) * u[[i]])))

    }
    stop("wrong argument 'u'")
}

dlmForecast <- function(mod, nAhead=1, method=c("plain","svd"), sampleNew=FALSE) {
    method <- match.arg(method)

    if (is.dlmFiltered(mod)) {
        modFuture <- mod$mod
        lastObsIndex <- NROW(mod$m)
        modFuture$C0 <- with(mod, dlmSvd2var(U.C[[lastObsIndex]], D.C[lastObsIndex,]))
        if (is.ts(mod$m))
            modFuture$m0 <- window(mod$m, start=end(mod$m))
        else {
            modFuture$m0 <- window(mod$m, start=lastObsIndex)
            tsp(modFuture$m0) <- NULL
        }
        mod <- modFuture
    }

    if (! (is.null(mod$JFF) && is.null(mod$JV) && is.null(mod$JGG) && is.null(mod$JW)))
        stop("dlmForecast only works with constant models")

    ytsp <- tsp(mod$m0)
    p <- length(mod$m0)
    m <- nrow(mod$FF)
    a <- rbind(mod$m0, matrix(0,nAhead,p))
    R <- vector("list",nAhead+1)
    R[[1]] <- mod$C0
    f <- matrix(0,nAhead,m)
    Q <- vector("list", nAhead)
    for (it in 1:nAhead) {
        a[it+1,] <- mod$GG %*% a[it,]
        R[[it+1]] <- mod$GG %*% R[[it]] %*% t(mod$GG) + mod$W
        f[it,] <- mod$FF %*% a[it+1,]
        Q[[it]] <- mod$FF %*% R[[it+1]] %*% t(mod$FF) + mod$V
    }
    a <- a[-1,,drop=FALSE]
    R <- R[-1]
    if ( sampleNew ) {
        newStates <- vector("list", sampleNew)
        newObs <- vector("list", sampleNew)
        newS <- matrix(0, nAhead, p)
        newO <- matrix(0, nAhead, m)
        tmp <- La.svd(mod$V,nu=0)
        Ut.V <- tmp$vt; D.V <- sqrt(tmp$d)
        tmp <- La.svd(mod$W,nu=0)
        Ut.W <- tmp$vt; D.W <- sqrt(tmp$d)
        for (i in 1:sampleNew) {
            tmp <- La.svd(R[[1]],nu=0)
            newS[1,] <- crossprod(tmp$vt, rnorm(p, sd=sqrt(tmp$d))) + a[1,]
            newO[1,] <- crossprod(Ut.V, rnorm(m, sd=D.V)) + mod$FF %*% newS[1,]
            if ( nAhead > 1 )
                for (it in 2:nAhead) {
                    newS[it,] <- crossprod(Ut.W, rnorm(p, sd=D.W)) + mod$GG %*% newS[it-1,]
                    newO[it,] <- crossprod(Ut.V, rnorm(m, sd=D.V)) + mod$FF %*% newS[it,]
                }
            newStates[[i]] <- newS
            newObs[[i]] <- newO
        }
        if (!is.null(ytsp)) {
            a <- ts(a, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3])
            f <- ts(f, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3])
            newStates <- lapply(newStates, function(x)
                                ts(x, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3]))
            newObs <- lapply(newObs, function(x)
                             ts(x, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3]))
        }
        ans <- list(a=a, R=R, f=f, Q=Q, newStates=newStates, newObs=newObs)
    } else {
        if (!is.null(ytsp)) {
            a <- ts(a, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3])
            f <- ts(f, start = ytsp[2] + 1/ytsp[3], frequency = ytsp[3])
        }
        ans <- list(a=a, R=R, f=f, Q=Q)
    }


    return(ans)
}

###### One-step forecast errors
residuals.dlmFiltered <- function(object, ...,
                                  type=c("standardized", "raw"), sd=TRUE) {
    if (is.null(object$y))
        stop("\'object\' argument has no \'y\' component")
    type <- match.arg(type)
    if (is.null(object$mod$JFF)) tvFF <- FALSE else tvFF <- TRUE
    if (is.null(object$mod$JV)) tvV <- FALSE else tvV <- TRUE
    FF <- object$mod$FF
    if (!( tvFF || tvV )) { ## constant model
        f <- object$a %*% t(FF)
        res <- drop(object$y - f) # one-step forecasting errors
        if (sd || (type == "standardized")) {
            V <- object$mod$V
            SD <- drop(t(sqrt(sapply(seq(along=object$U.R),
                                     function(i)
                                     diag(crossprod(object$D.R[i,] *
                                                    t(FF%*%object$U.R[[i]])) + V)))))
        }
    } else
    if ( !tvFF ) { ## only V time-varying
        f <- object$a %*% t(FF)
        res <- drop(object$y - f) # one-step forecasting errors
        if (sd || (type == "standardized")) {
            nz <- object$mod$JV != 0
            JV <- cbind(row(object$mod$JV)[nz], col(object$mod$JV)[nz],
                        object$mod$JV[nz])
            V <- object$mod$V
            getSD <- function(i) {
                V[JV[,-3,drop=FALSE]] <- object$mod$X[i,JV[,3]]
                diag(crossprod(object$D.R[i,] * t(FF%*%object$U.R[[i]])) + V)
            }
            SD <- drop(t(sqrt(sapply(seq(along=object$U.R), getSD))))
        }
    } else
    if ( !tvV ) { ## only FF time-varying
        if (!(sd || (type == "standardized"))) {
            nz <- object$mod$JFF != 0
            JFF <- cbind(row(object$mod$JFF)[nz], col(object$mod$JFF)[nz],
                         object$mod$JFF[nz])
            getFore <- function(i) {
                FF[JFF[,-3,drop=FALSE]] <- object$mod$X[i,JFF[,3]]
                FF %*% object$a[i,]
            }
            f <- drop(t(sapply(seq(along=object$U.R), getFore)))
            res <- drop(object$y - f) # one-step forecasting errors
        } else {
            nz <- object$mod$JFF != 0
            JFF <- cbind(row(object$mod$JFF)[nz], col(object$mod$JFF)[nz],
                         object$mod$JFF[nz])
            V <- object$mod$V
            getBoth <- function(i) {
                FF[JFF[,-3,drop=FALSE]] <- object$mod$X[i,JFF[,3]]
                c(FF %*% object$a[i,],
                  diag(crossprod(object$D.R[i,] * t(FF%*%object$U.R[[i]])) + V))
            }
            tmp <- t(sapply(seq(along=object$U.R), getBoth))
            m <- ncol(tmp) / 2
            res <- drop(object$y - tmp[,1:m])
            SD <- drop(sqrt(tmp[,-(1:m)]))
        }
    } else { ## both FF and V time-varying
        if (!(sd || (type == "standardized"))) {
            nz <- object$mod$JFF != 0
            JFF <- cbind(row(object$mod$JFF)[nz], col(object$mod$JFF)[nz],
                         object$mod$JFF[nz])
            getFore <- function(i) {
                FF[JFF[,-3,drop=FALSE]] <- object$mod$X[i,JFF[,3]]
                FF %*% object$a[i,]
            }
            f <- drop(t(sapply(seq(along=object$U.R), getFore)))
            res <- drop(object$y - f) # one-step forecasting errors
        } else {
            nz <- object$mod$JFF != 0
            JFF <- cbind(row(object$mod$JFF)[nz], col(object$mod$JFF)[nz],
                         object$mod$JFF[nz])
            nz <- object$mod$JV != 0
            JV <- cbind(row(object$mod$JV)[nz], col(object$mod$JV)[nz],
                        object$mod$JV[nz])
            V <- object$mod$V
            getBoth <- function(i) {
                FF[JFF[,-3,drop=FALSE]] <- object$mod$X[i,JFF[,3]]
                V[JV[,-3,drop=FALSE]] <- object$mod$X[i,JV[,3]]
                c(FF %*% object$a[i,],
                  diag(crossprod(object$D.R[i,] * t(FF%*%object$U.R[[i]])) + V))
            }
            tmp <- t(sapply(seq(along=object$U.R), getBoth))
            m <- ncol(tmp) / 2
            res <- drop(object$y - tmp[,1:m])
            SD <- drop(sqrt(tmp[,-(1:m)]))
        }
    }

    if ( type == "standardized" )
        res <- res / SD
    if (sd) {
        if (is.ts(res)) attributes(SD) <- attributes(res) # makes a time series of SD
        return(list(res=res, sd=SD))
    } else
    return(res)
}

###### Diagnostic plots
tsdiag.dlmFiltered <- function (object, gof.lag = 10, ...) {
    stdres <- residuals(object, sd=FALSE)
    if ((ns <- NCOL(stdres)) == 1) {
        oldpar <- par(mfrow = c(3, 1))
        on.exit(par(oldpar))
        plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
        abline(h = 0)
        acf(stdres, plot = TRUE, main = "ACF of Residuals",
            na.action = na.pass)
        nlag <- gof.lag
        pval <- numeric(nlag)
        for (i in 1:nlag) pval[i] <- Box.test(stdres, i, type = "Ljung-Box")$p.value
        plot(1:nlag, pval, xlab = "lag", ylab = "p value",
             ylim = c(0, 1), main = "p values for Ljung-Box statistic")
        abline(h = 0.05, lty = 2, col = "blue")
    } else {
        ask <- dev.interactive()
        oldpar <- par(mfrow = c(3, 1), oma=c(0, 0, 2, 0), "ask")
        on.exit(par(oldpar))
        hasNames <- !is.null(nm <- attr(stdres,"dimnames")[[2]])
        for (j in 1:ns) {
            plot(stdres[,j], type = "h", main = "Standardized Residuals", ylab = "")
            abline(h = 0)
            if (ask) par(ask=FALSE)
            acf(stdres[,j], plot = TRUE, main = "ACF of Residuals",
                na.action = na.pass)
            nlag <- gof.lag
            pval <- numeric(nlag)
            for (i in 1:nlag) pval[i] <- Box.test(stdres[,j], i, type = "Ljung-Box")$p.value
            plot(1:nlag, pval, xlab = "lag", ylab = "p value",
                 ylim = c(0, 1), main = "p values for Ljung-Box statistic")
            abline(h = 0.05, lty = 2, col = "blue")
            mtext(if (hasNames) nm[j] else paste("Series",j), line=1, outer=TRUE)
            if (ask) par(ask=TRUE)
        }
    }
}

###### Generating a random DLM
dlmRandom <- function(m, p, nobs = 0, JFF, JV, JGG, JW)
{
    ### Assume for now that each of FF, V, GG, W is either fixed or
    ### time-varying in every entry
    FF <- matrix(rnorm(m*p),m,p)
    V <- rwishart(2*m, m)
    GGtrial <- matrix(rnorm(p*p),p,p)
    e <- eigen(GGtrial)
    if ((ab <- max(abs(e$values))) > 1) {
        r <- runif(1)
        GG <- with(e, Re(vectors %*% (r * values / ab * solve(vectors))))
    }
    else
        GG <- GGtrial
    W <- rwishart(2*p, p)
    m0 <- rep(0,p)
    C0 <- diag(nrow = p) * 100
    if (nobs > 0) {
        count <- 0
        if (hasArg(JGG) && JGG == TRUE) {
            tvGG <- TRUE
            JGG <- structure(1:(p*p), dim=dim(GG))
            count <- p * p
        }
        else {
            tvGG <- FALSE
            JGG <- NULL
        }
        if (hasArg(JW) && JW == TRUE) {
            tvW <- TRUE
            JW <- structure(count + 1:(p*p), dim=dim(W))
            count <- count + p * p
        }
        else {
            tvW <- FALSE
            JW <- NULL
            tmp <- La.svd(W,nu=0)
            Ut.W <- tmp$vt; D.W <- sqrt(tmp$d)
        }
        if (hasArg(JFF) && JFF == TRUE) {
            tvFF <- TRUE
            JFF <- structure(count + 1:(m*p), dim=dim(FF))
            count <- count + m * p
        }
        else {
            tvFF <- FALSE
            JFF <- NULL
        }
        if (hasArg(JV) && JV == TRUE) {
            tvV <- TRUE
            JV <- structure(count + 1:(m*m), dim=dim(V))
        }
        else {
            tvV <- FALSE
            JV <- NULL
            tmp <- La.svd(V,nu=0)
            Ut.V <- tmp$vt; D.V <- sqrt(tmp$d)
        }
        if (any(c(tvFF, tvV, tvGG, tvW))) {
            newStates <- matrix(0, nrow = nobs + 1, ncol = p)
            newObs <- matrix(0, nrow = nobs, ncol = m)
            X <- matrix(0, nrow = nobs,
                        ncol = m * p * tvFF + m * m * tvV +
                        p * p * tvGG + p * p * tvW)
            tmp <- La.svd(C0,nu=0)
            newStates[1,] <- crossprod(tmp$vt, rnorm(p, sd=sqrt(tmp$d)))
            for (it in 1:nobs) {
                count <- 0
                if (tvGG) {
                    GGtrial <- matrix(rnorm(p*p),p,p)
                    e <- eigen(GGtrial)
                    if ((ab <- max(abs(e$values))) > 1) {
                        r <- runif(1)
                        GG <- with(e, Re(vectors %*% (r * values / ab * solve(vectors))))
                    }
                    else
                        GG <- GGtrial
                    X[it, 1:(p*p)] <- GG
                    count <- p * p
                }
                if (tvW) {
                    W <- rwishart(2*p, p)
                    X[it, count + 1:(p*p)] <- W
                    count <- count + p * p
                    tmp <- La.svd(W,nu=0)
                    Ut.W <- tmp$vt; D.W <- sqrt(tmp$d)
                }
                if (tvFF) {
                    FF <- matrix(rnorm(m*p),m,p)
                    X[it, count + 1:(m*p)] <- FF
                    count <- count + m * p
                }
                if (tvV) {
                    V <- rwishart(2*m, m)
                    X[it, count + 1:(m*m)] <- V
                    tmp <- La.svd(V,nu=0)
                    Ut.V <- tmp$vt; D.V <- sqrt(tmp$d)
                }
                newStates[it + 1,] <- crossprod(Ut.W, rnorm(p, sd=D.W)) + GG %*% newStates[it,]
                newObs[it,] <- crossprod(Ut.V, rnorm(m, sd=D.V)) + FF %*% newStates[it+1,]
            }
            mod <- dlm(list(m0=m0, C0=C0, FF=FF, V=V, GG=GG, W=W, JFF=JFF, JV=JV,
                            JGG=JGG, JW=JW, X=X))
            ans <- list(mod = mod, theta = newStates[-1,,drop=FALSE], y = newObs)
        }
        else {
            mod <- dlm(list(m0=m0, C0=C0, FF=FF, V=V, GG=GG, W=W))
            tmp <- dlmForecast(mod, nAhead = nobs, sampleNew = 1)
            ans <- list(mod = mod, theta = tmp$newStates[[1]], y = tmp$newObs[[1]])
        }
    }
    else
        ans <- dlm(list(m0=m0, C0=C0, FF=FF, V=V, GG=GG, W=W))
    return(ans)
}

###### Utility functions for MCMC output analysis
mcmcSD <- function(x) {
    ## Monte Carlo standard deviation using Sokal's method
    msg <- "Numeric vector or matrix argument needed"
    if (!is.numeric(x))
        stop(msg)
    univariate <- function(x) {
        l <- floor(30*log10(length(x)))
        ac <- drop(acf(x, lag.max=l, plot=FALSE)$acf)
        tau <- cumsum(c(1, 2 * ac[-1]))
        k <- 0
        while ( (k < 3 * tau[k+1]) &&  (k < l)) k <- k+1
        tau <- tau[k+1]
        if (tau < 0)
            {
                tau <- 1
                warning("estimated s.d. may not be reliable")
            }
        return(sqrt(var(x) * tau / length(x)))
    }
    if (is.null(dim(x)))
        ans <- univariate(x)
    else
        if (is.matrix(x))
            ans <- apply(x, 2, univariate)
        else
            stop(msg)
    return(ans)
}

mcmcMean <- function(x, sd = TRUE)
{
    ## Ergodic means of mcmc output
    ## Estimated standard deviations of means using Sokal's method
    msg <- "Numeric vector or matrix argument needed"
    if (!is.numeric(x))
        stop(msg)
    if (v <- is.null(dim(x))) # univariate input
    {
        mn <- mean(x)
        nm <- deparse(substitute(x))
    }
    else
        if (is.matrix(x)) # multivariate input
        {
            mn <- colMeans(x)
            nm <- colnames(x)
            if ( is.null(nm) )
                nm <- paste(deparse(substitute(x)), 1:NCOL(x), sep='.')
        }
        else
            stop(msg)
    if (sd) {
        univariateSD <- function(x) {
            l <- floor(30*log10(length(x)))
            ac <- acf(x, lag.max=l, plot=FALSE)$acf
            tau <- cumsum(c(1, 2 * ac[-1]))
            k <- 0
            while ( (k < 3 * tau[k+1]) &&  (k < l)) k <- k+1
            tau <- tau[k+1]
            return(sqrt(var(x) * tau / length(x)))
        }
        if ( v )
            sd <- univariateSD(x)
        else
            sd <- apply(x, 2, univariateSD)
        ans <- rbind(mn, sd)
        dimnames(ans) <- list(c("mean", "sd"), nm)
        class(ans) <- "mcmcMean"
        return(ans)
    } else {
        if ( !v ) names(mn) <- nm
        return(mn)
    }
}

mcmcMeans <- mcmcMean

print.mcmcMean <- function(x, digits = getOption("digits"), ...)
{
    ans <- format(x, digits = 3)
    ans[1, ] <- sapply(ans[1, ], function(x) paste("", x))
    ans[2, ] <- sapply(ans[2, ], function(x) paste("(", x, ")", sep = ""))
    dn <- dimnames(ans)
    dn[[1]] <- rep("", 2)
    if (is.null(dn[[2]]))
        dn[[2]] <- paste('V', 1:dim(ans)[2], sep='.')
    dn[[2]] <- paste(substring("      ", 1, (nchar(ans[2, ]) -
        nchar(dn[[2]]))%/%2), dn[[2]])
    dn[[2]] <- paste(dn[[2]], substring("      ", 1, (nchar(ans[2,
        ]) - nchar(dn[[2]]))%/%2))
    dimnames(ans) <- dn
    print(ans, quote = FALSE)
    return(x)
}

ergMean <- function(x, m = 1)
{
    ### ergodic means
    if (hasArg(m)) m <- max(1, round(m))
    n <- NROW(x)
    if ( m > n )
        stop("Need m <= n")
    if ( is.null(dm <- dim(x)) || dm[2] == 1 )
    {
        ## univariate
        if ( m == 1 )
            ans <- cumsum(x) / 1:n
        else
            if ( m == n )
                ans <- mean(x)
            else
                ans <- cumsum(c(sum(x[1:m]), x[(m+1):n])) / m:n
    } else {
        if ( length(dm) == 2 )
        {
            ## matrix "x" - multivariate
            nm <- colnames(x)
            if ( is.null(nm) )
                nm <- paste(deparse(substitute(x)), 1:NCOL(x), sep='.')
            if ( m == 1 )
                ans <- apply(x, 2, cumsum) / 1:n
            else
                if ( m == n )
                    ans <- matrix(colMeans(x), nrow = 1)
                else
                    ans <- apply(rbind(colSums(x[1:m,]), x[(m+1):n,]), 2, cumsum) / m:n
            colnames(ans) <- nm
        }
        else
            stop("\'x\' must be a vector or a matrix")
    }
    return(ans)
}


###### drop first element of a ts or mts retaining the class
dropFirst <- function(x)
{
    st <- start(x) + c(0, 1)
    freq <- frequency(x)
    newStart <- c(st[1], st[2] %% freq) + c(st[2] %/% freq, 0)
    y <- window(x, start = newStart)
    if (is.null(tsp(x)))
        attr(y, "tsp") <- NULL
    return(y)
}

######
###### Gibbs sampler for "d-inverse-gamma" model
######
dlmGibbsDIG <- function(y, mod, a.y, b.y, a.theta, b.theta, shape.y, rate.y,
                        shape.theta, rate.theta, n.sample = 1,
                        thin = 0, ind, save.states = TRUE,
                        progressBar = interactive())
##################################################################################
##################################################################################
### Gibbs sampler for the 'd-inverse-gamma' model                              ###
### Constant DLMs and univariate observations only                             ###
###                                                                            ###
### y       : data (vector or univariate time series).                         ###
### mod     : a dlm model for the data, with a diagonal 'W' component.         ###
### a.y     : prior mean of the observation precision.                         ###
### b.y     : prior variance of the observation precision.                     ###
### a.theta : vector of prior mean(s) of the system precision(s);              ###
###           recycled if needed.                                              ###
### b.theta : vector of prior variance(s) of the system precision(s);          ###
###           recycled if needed.                                              ###
### shape.y, rate.y : shape and rate parameters of the prior gamma             ###
###           distribution of the observation precision. Can be specified      ###
###           in alternative to 'a' and 'b'.                                   ###
### shape.theta, rate.theta : vectors of shape and rate parameters of the      ###
###           prior distributions of the system precision(s). Can be           ###
###           specified in alternative to 'alpha' and 'beta'.                  ###
### n.sample : number of simulated values in the output.                       ###
### thin    : number of sweeps to discard between each pair of returned        ###
###           draws.                                                           ###
### ind     : vector of indices. If specified, the sampler will only draw      ###
###           the diagonal elements of 'W' having the specified indices.       ###
###           Useful when some of the system variances are zero.               ###
### save.states : if TRUE, the generated states will be returned together      ###
###           with the generated parameter values.                             ###
### progressBar : if TRUE, a text progress bar will be displayed               ###
###           during execution                                                 ###
###                                                                            ###
### Value:                                                                     ###
### A list with components 'dV', 'dW', 'theta' (only if 'save.states' is       ###
### TRUE). 'dV' contains the generated observation variances, 'dW' the         ###
### generated system variances, 'theta' the generated states.                  ###
##################################################################################
##################################################################################
{
    msg1 <- "Either \"a.y\" and \"b.y\" or \"shape.y\" and \"rate.y\" must be specified"
    msg2 <- "Unexpected length of \"shape.y\" and/or \"rate.y\""
    msg3 <- "Unexpected length of \"a.y\" and/or \"b.y\""
    msg4 <- paste("Either \"a.theta\" and \"b.theta\" or \"shape.theta\"",
                  "and \"rate.theta\" must be specified")
    msg5 <- "Unexpected length of \"shape.theta\" and/or \"rate.theta\""
    msg6 <- "Unexpected length of \"a.theta\" and/or \"b.theta\""
    msg7 <- "\"thin\" must be a nonnegative integer"
    msg8 <- "multivariate observations are not allowed"
    msg9 <- "inadmissible value of \"ind\""
    if ( NCOL(y) > 1 )
        stop(msg8)
    r <- ncol(mod$FF)
    if ( hasArg(ind) ) {
        ind <- unique(as.integer(ind))
        s <- 1:r
        if ( !all(ind %in% s) )
            stop(msg9)
        perm <- s[c(ind, s[ !(s %in% ind)])]
        FF(mod) <- mod$FF[, perm, drop = FALSE]
        GG(mod) <- mod$GG[perm, perm, drop = FALSE]
        W(mod) <- mod$W[perm, perm, drop = FALSE]
        p <- length(ind)
    }
    else {
        perm <- ind <- 1 : r
        p <- r
    }
    nobs <- NROW(y)
    if ( is.numeric(thin) && (thin <- as.integer(thin)) >= 0 )
    {
        every <- thin + 1
        mcmc <- n.sample * every
    }
    else
        stop(msg7)
    ## check hyperpriors for precision of 'y'
    if ( !hasArg(a.y) )
        if ( !hasArg(shape.y) ) stop(msg1)
        else
            if ( !hasArg(rate.y) ) stop(msg1)
            else
            {
                ## check length of shape.y and rate.y
                if (!all(c(length(shape.y), length(rate.y)) == 1))
                    warning(msg2)
            }
    else
        if ( !hasArg(b.y) ) stop(msg1)
        else
        {
            if (!all(c(length(a.y), length(b.y)) == 1))
                warning(msg3)
            shape.y <- a.y^2 / b.y
            rate.y <- a.y / b.y
        }
    ## check hyperpriors for precision(s) of 'theta'
    if ( !hasArg(a.theta) )
        if ( !hasArg(shape.theta) ) stop(msg4)
        else
            if ( !hasArg(rate.theta) ) stop(msg4)
            else
            {
                ## check length of shape.theta and rate.theta
                if (!all(c(length(shape.theta), length(rate.theta)) %in% c(1,p)))
                    warning(msg5)
            }
    else
        if ( !hasArg(b.theta) ) stop(msg4)
        else
        {
            if (!all(c(length(a.theta), length(b.theta)) %in% c(1,p)))
                warning(msg6)
            shape.theta <- a.theta^2 / b.theta
            rate.theta <- a.theta / b.theta
        }
    shape.y <- shape.y + 0.5 * nobs
    shape.theta <- shape.theta + 0.5 * nobs
    shape.theta <- rep(shape.theta, length.out = p)
    rate.theta <- rep(rate.theta, length.out = p)
    theta <- matrix(0, nobs + 1, r)
    if ( save.states )
        gibbsTheta <- array(0, dim = c(nobs + 1, r, n.sample))
    gibbsV <- vector("numeric", n.sample)
    gibbsW <- matrix(0, nrow = n.sample, ncol = p)
    it.save <- 0
    if (progressBar) pb <- txtProgressBar(0, mcmc, style = 3)
    for (it in 1:mcmc)
    {
        if (progressBar) setTxtProgressBar(pb, it)
        ## generate states - FFBS
        modFilt <- dlmFilter(y, mod, simplify=TRUE)
        theta[] <- dlmBSample(modFilt)
        ## generate V
        y.center <- y - tcrossprod(theta[-1,,drop=FALSE], mod$FF)
        SSy <- drop(crossprod(y.center))
        mod$V[] <- 1 / rgamma(1, shape = shape.y,
                              rate = rate.y + 0.5 * SSy)
        ## generate W
        theta.center <- theta[-1,,drop=FALSE] -
            tcrossprod(theta[-(nobs + 1),,drop=FALSE], mod$GG)
        SStheta <- drop(sapply( 1 : p, function(i) crossprod(theta.center[,i])))
        SStheta <- colSums((theta[-1,1:p,drop=FALSE] -
                            tcrossprod(theta[-(nobs + 1),,drop=FALSE],mod$GG)[,1:p])^2)
        diag(mod$W)[1:p] <-
            1 / rgamma(p, shape = shape.theta, rate = rate.theta + 0.5 * SStheta)
        ## save
        if ( !(it %% every) )
        {
            it.save <- it.save + 1
            if ( save.states )
                gibbsTheta[,,it.save] <- theta
            gibbsV[it.save] <- diag(mod$V)
            gibbsW[it.save,] <- diag(mod$W)[1:p]
        }
    }
    colnames(gibbsW) <- paste("W", ind, sep = '.')
    if (progressBar) close(pb)
    if ( save.states )
        return(list(dV = gibbsV, dW = gibbsW,
                    theta = gibbsTheta[, order(perm), , drop = FALSE]))
    else
        return(list(dV = gibbsV, dW = gibbsW))
}

Try the dlm package in your browser

Any scripts or data that you put into this service are public.

dlm documentation built on June 14, 2018, 1:03 a.m.