# R/DLM.R In dlm: Bayesian and Likelihood Analysis of Dynamic Linear Models

```###### 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
m0=rep(0,length(dW)), C0=1e7*diag(nrow=length(dW)))
{
## sanity checks
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)
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)
}

"+.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)
R[[1]] <- mod\$C0
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)
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 )
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 {
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)
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)
}
}
}

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