add.dynamic.regression {bsts}

Dynamic Regression State Component

Description

Add a dynamic regression component to the state specification of a bsts model. A dynamic regression is a regression model where the coefficients change over time according to a random walk.

Usage

AddDynamicRegression(
    state.specification,
    formula,
    data,
    model.options = NULL,
    sigma.mean.prior.DEPRECATED = NULL,
    shrinkage.parameter.prior.DEPRECATED = GammaPrior(a = 10, b = 1),
    sigma.max.DEPRECATED = NULL,
    contrasts = NULL,
    na.action = na.pass)

DynamicRegressionRandomWalkOptions(
    sigma.prior = NULL,
    sdy = NULL,
    sdx = NULL)

DynamicRegressionHierarchicalRandomWalkOptions(
     sdy = NULL,
     sigma.mean.prior = NULL,
     shrinkage.parameter.prior = GammaPrior(a = 10, b = 1),
     sigma.max = NULL)

DynamicRegressionArOptions(lags = 1, sigma.prior = SdPrior(1, 1))

Arguments

state.specification

A list of state components that you wish to add to. If omitted, an empty list will be assumed.

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 vector giving the time series to be modeled.

data

An optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from 'environment(formula)', typically the environment from which AddDynamicRegression is called.

model.options

An object inheriting from DynamicRegressionOptions giving the specific transition model for the dynamic regression coefficients, and the prior distribution for any hyperparameters associated with the transition model.

sigma.mean.prior

An object created by GammaPrior describing the prior distribution of b/a (see details below).

sigma.mean.prior.DEPRECATED

This option should be set using model.options. It will be removed in a future release.

shrinkage.parameter.prior

An object of class GammaPrior describing the shrinkage parameter, a (see details below).

shrinkage.parameter.prior.DEPRECATED

This option should be set using model.options. It will be removed in a future release.

sigma.max

The largest supported value of each sigma[i]. Truncating the support of sigma can keep ill-conditioned models from crashing. This must be a positive number (Inf is okay), or NULL. A NULL value will set sigma.max = sd(y), which is a substantially larger value than one would expect, so in well behaved models this constraint will not affect the analysis.

sigma.max.DEPRECATED

This option should be set using model.options. It will be removed in a future release.

contrasts

An optional list. See the contrasts.arg of model.matrix.default. This argument is only used if a model formula is specified. It can usually be ignored even then.

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.

sdy

The standard deviation of the response variable. This is used to scale default priors and sigma.max if other arguments are left NULL. If all other arguments are non-NULL then sdy is not used.

sdx

The vector of standard deviations of each predictor variable in the dynamic regression. Used only to scale the default prior. This argument is not used if a prior is specified directly.

lags

The number of lags in the autoregressive process for the coefficients.

sigma.prior

Either an object of class SdPrior or a list of such objects. If a single SdPrior is given then it specifies the prior on the innovation variance for all the coefficients. If a list of SdPrior objects is given, then each element gives the prior distribution for the corresponding regression coefficient. The length of such a list must match the number of predictors in the dynamic regression part of the model.

Details

For the standard "random walk" coefficient model, the model is

\beta_{i, t+1} = beta_{i, t} + \epsilon_t \qquad \epsilon_t \sim \mathcal{N}(0, \sigma^2_i / variance_{xi})

\frac{1}{\sigma^2_i} \sim Ga(a, b)

\sqrt{b/a} \sim sigma.mean.prior

a \sim shrinkage.parameter.prior

That is, each coefficient evolves independently, with its own variance term which is scaled by the variance of the i'th column of X. The parameters of the hyperprior are interpretable as: sqrt(b/a) typical amount that a coefficient might change in a single time period, and 'a' is the 'sample size' or 'shrinkage parameter' measuring the degree of similarity in sigma[i] among the arms.

In most cases we hope b/a is small, so that sigma[i]'s will be small and the series will be forecastable. We also hope that 'a' is large because it means that the sigma[i]'s will be similar to one another.

The default prior distribution is a pair of independent Gamma priors for sqrt(b/a) and a. The mean of sigma[i] is set to .01 * sd(y) with shape parameter equal to 1. The mean of the shrinkage parameter is set to 10, but with shape parameter equal to 1.

If the coefficients have AR dynamics, then the model is that each coefficient independently follows an AR(p) process, where p is given by the lags argument. Independent priors are assumed for each coefficient's model, with a uniform prior on AR coefficients (with support restricted to the finite region where the process is stationary), while the sigma.prior argument gives the prior for each coefficient's innovation variance.

Value

Returns a list with the elements necessary to specify a dynamic regression model.

Author(s)

Steven L. Scott

References

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.

See Also

bsts. SdPrior NormalPrior

Examples

## Setting the seed to avoid small sample effects resulting from small
## number of iterations.
set.seed(8675309)
n <- 1000
x <- matrix(rnorm(n))

# beta follows a random walk with sd = .1 starting at -12.
beta <- cumsum(rnorm(n, 0, .1)) - 12

# level is a local level model with sd = 1 starting at 18.
level <- cumsum(rnorm(n)) + 18

# sigma.obs is .1
error <- rnorm(n, 0, .1)

y <- level + x * beta + error
par(mfrow = c(1, 3))
plot(y, main = "Raw Data")
plot(x, y - level, main = "True Regression Effect")
plot(y - x * beta, main = "Local Level Effect")

ss <- list()
ss <- AddLocalLevel(ss, y)
ss <- AddDynamicRegression(ss, y ~ x)
## In a real appliction you'd probably want more than 100
## iterations. See comment above about the random seed.
model <- bsts(y, state.specification = ss, niter = 100, seed = 8675309)
plot(model, "dynamic", burn = 10)

xx <- rnorm(10)
pred <- predict(model, newdata = xx)
plot(pred)