auto.ar: Sparse AR(p)

Description

Add a sparse AR(p) process to the state distribution. A sparse AR(p) is an AR(p) process with a spike and slab prior on the autoregression coefficients.

Usage

AddAutoAr(state.specification,
          y,
          lags = 1,
          prior = NULL,
          sdy = NULL,
          ...)

Value

Returns state.specification with an AR(p) state component added to the end.

Arguments

state.specification: A list of state components. If omitted, an empty list is assumed.
y: A numeric vector. The time series to be modeled. This can be omitted if sdy is supplied.
lags: The maximum number of lags ("p") to be considered in the AR(p) process.
prior: An object inheriting from SpikeSlabArPrior, or NULL. If the latter, then a default SpikeSlabArPrior will be created.
sdy: The sample standard deviation of the time series to be modeled. Used to scale the prior distribution. This can be omitted if y is supplied.
...: Extra arguments passed to SpikeSlabArPrior.

Author

Steven L. Scott steve.the.bayesian@gmail.com

Details

The model contributes alpha[t] to the expected value of y[t], where the transition equation is

$α_{t} = ϕ_{1} α_{i, t - 1} + \dots + ϕ_{p} α_{t - p} + ϵ_{t - 1} ϵ_{t} \sim N (0, σ^{2})$

The state consists of the last p lags of alpha. The state transition matrix has phi in its first row, ones along its first subdiagonal, and zeros elsewhere. The state variance matrix has sigma^2 in its upper left corner and is zero elsewhere. The observation matrix has 1 in its first element and is zero otherwise.

This model differs from the one in AddAr only in that some of its coefficients may be set to zero.

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.

Examples

Run this code

n <- 100
residual.sd <- .001

# Actual values of the AR coefficients
true.phi <- c(-.7, .3, .15)
ar <- arima.sim(model = list(ar = true.phi),
                n = n,
                sd = 3)

## Layer some noise on top of the AR process.
y <- ar + rnorm(n, 0, residual.sd)
ss <- AddAutoAr(list(), y, lags = 6)

# Fit the model with knowledge with residual.sd essentially fixed at the
# true value.
model <- bsts(y, state.specification=ss, niter = 500, prior = SdPrior(residual.sd, 100000))

# Now compare the empirical ACF to the true ACF.
acf(y, lag.max = 30)
points(0:30, ARMAacf(ar = true.phi, lag.max = 30), pch = "+")
points(0:30, ARMAacf(ar = colMeans(model$AR6.coefficients), lag.max = 30))
legend("topright", leg = c("empirical", "truth", "MCMC"), pch = c(NA, "+", "o"))

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning