mbsts {bsts}  R Documentation 
Fit a multivariate Bayesian structural time series model, also known as a "dynamic factor model."
** NOTE ** This code is experimental. Please feel free to experiment with it and report any bugs to the maintainer. Expect it to improve substantially in the next release.
mbsts(formula, shared.state.specification, series.state.specification = NULL, data = NULL, timestamps = NULL, series.id = NULL, prior = NULL, # TODO opts = NULL, contrasts = NULL, na.action = na.pass, niter, ping = niter / 10, data.format = c("long", "wide"), seed = NULL, ...)
formula 
A formula describing the regression portion of the relationship between y and X. If no regressors are desired then the formula can be replaced by a numeric matrix giving the multivariate time series to be modeled. 
shared.state.specification 
A list with elements created by
This list defines the components of state which are shared across all time series. These are the "factors" in the dynamic factor model. 
series.state.specification 
This argument specifies state
components needed by a particular series. Not all series need have
the same state components (e.g. some series may require a seasonal
component, while others do not). It can be It can be a list of elements created by In its most general form, this argument can be a list of lists, some of which can be NULL, but with nonNULL lists specifying state components for individual series, as above. 
data 
An optional data frame, list or environment (or object
coercible by 
timestamps 
A vector of timestamps indicating the time of each
observation. If TODO: TEST THIS under wide and long formats in regression and nonregression settings. 
series.id 
A factor (or object coercible to factor) indicating the series to which each observation in "long" format belongs. This argument is ignored for data in "wide" format. 
prior 
A list of 
opts 
A list containing model options. This is currently only
used for debugging, so leave this as 
contrasts 
An optional list containing the names of contrast
functions to use when converting factors numeric variables in a
regression formula. This argument works exactly as it does in

na.action 
What to do about missing values. The default is to allow missing responses, but no missing predictors. Set this to na.omit or na.exclude if you want to omit missing responses altogether. 
niter 
A positive integer giving the desired number of MCMC draws. 
ping 
A scalar giving the desired frequency of status messages.
If ping > 0 then the program will print a status message to the
screen every 
data.format 
Whether the data are store in wide (each row is a
time point, and columns are values from different series) or long
(each row is the value of a particular series at a particular point
in time) format. For 
seed 
An integer to use as the random seed for the underlying
C++ code. If 
... 
Extra arguments to be passed to

An object of class mbsts
which is a list with the
following components
coefficients 
A 
sigma.obs 
A vector of length 
The returned object will also contain named elements holding the MCMC
draws of model parameters belonging to the state models. The names of
each component are supplied by the entries in
state.specification
. If a model parameter is a scalar, then
the list element is a vector with niter
elements. If the
parameter is a vector then the list element is a matrix with
niter
rows. If the parameter is a matrix then the list element
is a 3way array with first dimension niter
.
Finally, if a model formula was supplied, then the returned object will contain the information necessary for the predict method to build the design matrix when a new prediction is made.
Steven L. Scott steve.the.bayesian@gmail.com
Harvey (1990), "Forecasting, structural time series, and the Kalman filter", Cambridge University Press.
Durbin and Koopman (2001), "Time series analysis by state space methods", Oxford University Press.
George and McCulloch (1997) "Approaches for Bayesian variable selection", Statistica Sinica pp 339–374.
bsts
,
AddLocalLevel
,
AddLocalLinearTrend
,
AddSemilocalLinearTrend
,
AddSeasonal
AddDynamicRegression
SpikeSlabPrior
,
SdPrior
.
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) print("observation.coefficients =\n", observation.coefficients) 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) ## Not run: cat("simulated y (first 10 rows): \n") print(y[1:10, ]) ## ## 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") ## End(Not run) ## ## 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. ## ## Not run: 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) ## End(Not run) ## ## 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)