R: mbsts documentaion error – example usage fails.

[edwardstm]

The CRAN documentation for the bsts package contains the following example usage for mbsts:

seed <- 8675309
set.seed(seed)
ntimes <- 250
nseries <- 20
nfactors <- 6
residual.sd <- 1.2
state.innovation.sd <- .75
##---------------------------------------------------------------------------
## simulate latent state for fake data.
##---------------------------------------------------------------------------
state <- matrix(rnorm(ntimes * nfactors, 0, state.innovation.sd), nrow = ntimes)
for (i in 1:ncol(state)) {
    state[, i] <- cumsum(state[, i])
}
##---------------------------------------------------------------------------
## Simulate "observed" data from state.
##---------------------------------------------------------------------------
observation.coefficients <- matrix(rnorm(nseries * nfactors), nrow = nseries)
diag(observation.coefficients) <- 1.0
observation.coefficients[upper.tri(observation.coefficients)] <- 0
errors <- matrix(rnorm(nseries * ntimes, 0, residual.sd), ncol = ntimes)
y <- t(observation.coefficients %*% t(state) + errors)
##---------------------------------------------------------------------------
## Plot the data.
##---------------------------------------------------------------------------
par(mfrow=c(1,2))
plot.ts(y, plot.type="single", col = rainbow(nseries), main = "observed data")
plot.ts(state, plot.type = "single", col = 1:nfactors, main = "latent state")
##---------------------------------------------------------------------------
## Fit the model
##---------------------------------------------------------------------------
ss <- AddSharedLocalLevel(list(), y, nfactors = nfactors)
opts <- list("fixed.state" = t(state),
fixed.residual.sd = rep(residual.sd, nseries),
fixed.regression.coefficients = matrix(rep(0, nseries), ncol = 1))
model <- mbsts(y, shared.state.specification = ss, niter = 100,
data.format = "wide", seed = seed)
##---------------------------------------------------------------------------
## Plot the state
##---------------------------------------------------------------------------
par(mfrow=c(1, nfactors))
ylim <- range(model$shared.state, state)
for (j in 1:nfactors) {
    PlotDynamicDistribution(model$shared.state[, j, ], ylim=ylim)
    lines(state[, j], col = "blue")
}
##---------------------------------------------------------------------------
## Plot the factor loadings.
##---------------------------------------------------------------------------
opar <- par(mfrow=c(nfactors,1), mar=c(0, 4, 0, 4), omi=rep(.25, 4))
burn <- 10
for(j in 1:nfactors) {
    BoxplotTrue(model$shared.local.level.coefficients[-(1:burn), j, ],
        t(observation.coefficients[, j]), axes=F, truth.color="blue")
    abline(h=0, lty=3)
    box()
    axis(2)
}
axis(1)
par(opar)
##---------------------------------------------------------------------------
## Plot the predicted values of the series.
##---------------------------------------------------------------------------
index <- 1:12
nr <- floor(sqrt(length(index)))
nc <- ceiling(length(index) / nr)
opar <- par(mfrow = c(nr, nc), mar = c(2, 4, 1, 2))
for (i in index) {
  PlotDynamicDistribution(
      model$shared.state.contributions[, 1, i, ]
      + model$regression.coefficients[, i, 1]
      , ylim=range(y))
  points(y[, i], col="blue", pch = ".", cex = .2)
}
par(opar)

The "Plot the factor loadings." section does not function as written;

Error in BoxplotTrue(model$shared.local.level.coefficients[-(1:burn),  : 
  truth and x don't conform in boxplot.true
Some(<code style = 'font-size:10p'> Error in BoxplotTrue(model$shared.local.level.coefficients[-(1:burn), : truth and x don't conform in boxplot.true </code>)

The provided x (from model$shared.local.level.coefficients) is a (90, 6) matrix, while the truth (from observation.coefficients) is (1, 20). The full observation.coefficients is (20, 6), and a length 6 truth makes sense, so I believe there's a misplaced parenthesis; changing the BoxplotTrue line to the following results in plots:

    BoxplotTrue(model$shared.local.level.coefficients[-(1:burn), j, ],
        t(observation.coefficients)[, j], axes=F, truth.color="blue")