bsts

Description

Time series regression using dynamic linear models fit using MCMC. See Scott and Varian (2014) <DOI:10.1504/IJMMNO.2014.059942>, among many other sources.

Details

Installation note for Linux users

If you are installing bsts using install.packages on a Linux machine (and thus compiling yourself) you will almost certainly want to set the Ncpus argument to a large number. Windows and Mac users can ignore this advice.

Author(s)

Author: Steven L. Scott <steve.the.bayesian@gmail.com> Maintainer: Steven L. Scott <steve.the.bayesian@gmail.com>

References

Please see the references in the help page for the bsts function.

See Also

See the examples in the bsts function.

HarveyCumulator

Description

Given a state space model on a fine scale, the Harvey cumulator aggregates the model to a coarser scale (e.g. from days to weeks, or weeks to months).

Usage

HarveyCumulator(fine.series,
                contains.end,
                membership.fraction)

Arguments

Value

Returns a vector containing the course scale partial aggregates of fine.series.

Author(s)

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

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.mixed,

Examples

data(goog)
days <- factor(weekdays(index(goog)),
               levels = c("Monday", "Tuesday", "Wednesday",
                          "Thursday", "Friday"),
               ordered = TRUE)

## Because of holidays, etc the days do not always go in sequence.
## (Sorry, Rebecca Black! https://www.youtube.com/watch?v=kfVsfOSbJY0)
## diff.days[i] is the number of days between days[i-1] and days[i].
## We know that days[i] is the end of a week if diff.days[i] < 0.
diff.days <- tail(as.numeric(days), -1) - head(as.numeric(days), -1)
contains.end <- c(FALSE, diff.days < 0)

goog.weekly <- HarveyCumulator(goog, contains.end, 1)

Match Numeric Timestamps

Description

S3 generic method for MATCH function supplied in the zoo package.

Usage

  ## S3 method for class 'NumericTimestamps'
MATCH(x, table, nomatch = NA, ...)

Arguments

Details

Numeric timestamps match if they agree to 8 significant digits.

Value

Returns the index of the entry in table matched by each argument in x. If an entry has no match then nomatch is returned at that position.

See Also

MATCH

Add a state component to a Bayesian structural time series model

Description

Add a state component to the state.specification argument in a bsts model.

Author(s)

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

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 Ar1CoefficientPrior

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

Suggested burn-in size

Description

Suggest the size of an MCMC burn in sample as a proportion of the total run.

Usage

SuggestBurn(proportion, bsts.object)

Arguments

Value

An integer number of iterations to discard.

See Also

bsts

AR(p) state component

Description

Add an AR(p) state component to the state specification.

Usage

AddAr(state.specification,
      y,
      lags = 1,
      sigma.prior,
      initial.state.prior = NULL,
      sdy)

Arguments

Details

The model is

\alpha_{t} = \phi_1\alpha_{i, t-1} + \cdots + \phi_p \alpha_{t-p} + \epsilon_{t-1} \qquad \epsilon_t \sim \mathcal{N}(0, \sigma^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.

Value

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

Author(s)

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

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

Examples

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 <- AddAr(list(), lags = 3, sigma.prior = SdPrior(3.0, 1.0))

# 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$AR3.coefficients), lag.max = 30))
legend("topright", leg = c("empirical", "truth", "MCMC"), pch = c(NA, "+", "o"))

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

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)

Local level trend state component

Description

Add a local level model to a state specification. The local level model assumes the trend is a random walk:

\alpha_{t+1} = \alpha_t + \epsilon_t \qquad \epsilon_t \sim \mathcal{N}(0,\sigma).

The prior is on the \sigma parameter.

Usage

  AddLocalLevel(
     state.specification,
     y,
     sigma.prior,
     initial.state.prior,
     sdy,
     initial.y)

Arguments

Value

Returns a list with the elements necessary to specify a local linear trend state model.

Author(s)

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

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

Local linear trend state component

Description

Add a local linear trend model to a state specification. The local linear trend model assumes that both the mean and the slope of the trend follow random walks. The equation for the mean is

\mu_{t+1} = \mu_t + \delta_t + \epsilon_t \qquad \epsilon_t \sim \mathcal{N}(0, \sigma_\mu).

The equation for the slope is

\delta_{t+1} = \delta_t + \eta_t \qquad \eta_t \sim \mathcal{N}(0, \sigma_\delta).

The prior distribution is on the level standard deviation \sigma_\mu and the slope standard deviation \sigma_\delta.

Usage

  AddLocalLinearTrend(
     state.specification = NULL,
     y,
     level.sigma.prior = NULL,
     slope.sigma.prior = NULL,
     initial.level.prior = NULL,
     initial.slope.prior = NULL,
     sdy,
     initial.y)

Arguments

Value

Returns a list with the elements necessary to specify a local linear trend state model.

Author(s)

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

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

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

Monthly Annual Cycle State Component

Description

A seasonal state component for daily data, representing the contribution of each month to the annual seasonal cycle. I.e. this is the "January, February, March, ..." effect, with 12 seasons. There is a step change at the start of each month, and then the contribution of that month is constant over the course of the month.

Note that if you have anything other than daily data, then you're probably looking for AddSeasonal instead.

The state of this model is an 11-vector \gamma_t where the first element is the contribution to the mean for the current month, and the remaining elements are the values for the 10 most recent months. When t is the first day in the month then

\gamma_{t+1} = -\sum_{i = 2}^11 \gamma_{t, i} + % \epsilon_t \qquad \epsilon_t \sim \mathcal{N}(0, \sigma)

And the remaining elements are \gamma_t shifted down one. When t is any other day then \gamma_{t+1} = % \gamma_t.

Usage

AddMonthlyAnnualCycle(state.specification,
                      y,
                      date.of.first.observation = NULL,
                      sigma.prior = NULL,
                      initial.state.prior = NULL,
                      sdy)

Arguments

Examples

  ## Let's simulate some fake daily data with a monthly cycle.
## Not run: 
  residuals <- rnorm(365 * 5)

## End(Not run)

  n <- length(residuals)
  dates <- seq.Date(from = as.Date("2014-01-01"),
                    len = n,
                    by = 1)
  monthly.cycle <- rnorm(12)
  monthly.cycle <- monthly.cycle - mean(monthly.cycle)
  timestamps <- as.POSIXlt(dates)
  month <- timestamps$mon + 1

  new.month <- c(TRUE, diff(timestamps$mon) != 0)
  month.effect <- cumsum(new.month) 
  month.effect[month.effect == 0] <- 12

  response <- monthly.cycle[month] + residuals
  response <- zoo(response, timestamps)

  ## Now let's fit a bsts model to the daily data with a monthly annual
  ## cycle.
  ss <- AddLocalLevel(list(), response)
  ss <- AddMonthlyAnnualCycle(ss, response)

  ## In real life you'll probably want more iterations.
  model <- bsts(response, state.specification = ss, niter = 200)
  plot(model)
  plot(model, "monthly")

Random Walk Holiday State Model

Description

Adds a random walk holiday state model to the state specification. This model says

% y_t = \alpha_{d(t), t} + \epsilon_t

where there is one element in \alpha_t for each day in the holiday influence window. The transition equation is

% \alpha_{d(t+1), t+1} = \alpha_{d(t+1), t} + \epsilon_{t+1}

if t+1 occurs on day d(t+1) of the influence window, and

\alpha_{d(t+1), t+1} = \alpha_{d(t+1), t} %

otherwise.

Usage

AddRandomWalkHoliday(state.specification = NULL,
                     y,
                     holiday,
                     time0 = NULL, 
                     sigma.prior = NULL,
                     initial.state.prior = NULL,
                     sdy = sd(as.numeric(y), na.rm = TRUE))

Arguments

Value

A list describing the specification of the random walk holiday state model, formatted as expected by the underlying C++ code.

Author(s)

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

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. RegressionHolidayStateModel HierarchicalRegressionHolidayStateModel

Examples

trend <- cumsum(rnorm(730, 0, .1))
dates <- seq.Date(from = as.Date("2014-01-01"), length = length(trend),
  by = "day")
y <- zoo(trend + rnorm(length(trend), 0, .2), dates)

AddHolidayEffect <- function(y, dates, effect) {
  ## Adds a holiday effect to simulated data.
  ## Args:
  ##   y: A zoo time series, with Dates for indices.
  ##   dates: The dates of the holidays.
  ##   effect: A vector of holiday effects of odd length.  The central effect is
  ##     the main holiday, with a symmetric influence window on either side.
  ## Returns:
  ##   y, with the holiday effects added.
  time <- dates - (length(effect) - 1) / 2
  for (i in 1:length(effect)) {
    y[time] <- y[time] + effect[i]
    time <- time + 1
  }
  return(y)
}

## Define some holidays.
memorial.day <- NamedHoliday("MemorialDay")
memorial.day.effect <- c(.3, 3, .5)
memorial.day.dates <- as.Date(c("2014-05-26", "2015-05-25"))
y <- AddHolidayEffect(y, memorial.day.dates, memorial.day.effect)

presidents.day <- NamedHoliday("PresidentsDay")
presidents.day.effect <- c(.5, 2, .25)
presidents.day.dates <- as.Date(c("2014-02-17", "2015-02-16"))
y <- AddHolidayEffect(y, presidents.day.dates, presidents.day.effect)

labor.day <- NamedHoliday("LaborDay")
labor.day.effect <- c(1, 2, 1)
labor.day.dates <- as.Date(c("2014-09-01", "2015-09-07"))
y <- AddHolidayEffect(y, labor.day.dates, labor.day.effect)

## The holidays can be in any order.
holiday.list <- list(memorial.day, labor.day, presidents.day)
number.of.holidays <- length(holiday.list)

## In a real example you'd want more than 100 MCMC iterations.
niter <- 100
ss <- AddLocalLevel(list(), y)
ss <- AddRandomWalkHoliday(ss, y, memorial.day)
ss <- AddRandomWalkHoliday(ss, y, labor.day)
ss <- AddRandomWalkHoliday(ss, y, presidents.day)
model <- bsts(y, state.specification = ss, niter = niter, seed = 8675309)

## Plot model components.
plot(model, "comp")

## Plot the effect of the specific state component.
plot(ss[[2]], model)

Seasonal State Component

Description

Add a seasonal model to a state specification.

The seasonal model can be thought of as a regression on nseasons dummy variables with coefficients constrained to sum to 1 (in expectation). If there are S seasons then the state vector \gamma is of dimension S-1. The first element of the state vector obeys

\gamma_{t+1, 1} = -\sum_{i = 2}^S \gamma_{t, i} + \epsilon_t \qquad \epsilon_t \sim \mathcal{N}(0, \sigma)

Usage

  AddSeasonal(
     state.specification,
     y,
     nseasons,
     season.duration = 1,
     sigma.prior,
     initial.state.prior,
     sdy)

Arguments

Value

Returns a list with the elements necessary to specify a seasonal state model.

Author(s)

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

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

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

Semilocal Linear Trend

Description

The semi-local linear trend model is similar to the local linear trend, but more useful for long-term forecasting. It assumes the level component moves according to a random walk, but the slope component moves according to an AR1 process centered on a potentially nonzero value D. The equation for the level is

\mu_{t+1} = \mu_t + \delta_t + \epsilon_t \qquad \epsilon_t \sim \mathcal{N(0, \sigma_\mu)}.

The equation for the slope is

\delta_{t+1} = D + \phi (\delta_t - D) + \eta_t \qquad \eta_t \sim \mathcal{N(0, \sigma_\delta)}.

This model differs from the local linear trend model in that the latter assumes the slope \delta_t follows a random walk. A stationary AR(1) process is less variable than a random walk when making projections far into the future, so this model often gives more reasonable uncertainty estimates when making long term forecasts.

The prior distribution for the semi-local linear trend has four independent components. These are:

• an inverse gamma prior on the level standard deviation \sigma_\mu,

• an inverse gamma prior on the slope standard deviation \sigma_\delta,

• a Gaussian prior on the long run slope parameter D,

• and a potentially truncated Gaussian prior on the AR1 coefficient \phi. If the prior on \phi is truncated to (-1, 1), then the slope will exhibit short term stationary variation around the long run slope D.

Usage

  AddSemilocalLinearTrend(
     state.specification = list(),
     y = NULL,
     level.sigma.prior = NULL,
     slope.mean.prior = NULL,
     slope.ar1.prior = NULL,
     slope.sigma.prior = NULL,
     initial.level.prior = NULL,
     initial.slope.prior = NULL,
     sdy = NULL,
     initial.y = NULL)

Arguments

Value

Returns a list with the elements necessary to specify a generalized local linear trend state model.

Author(s)

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

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

Local level trend state component

Description

Add a shared local level model to a state specification. The shared local level model assumes the trend is a multivariate random walk:

\alpha_{t+1} = \alpha_t + \eta_t \qquad \eta_{tj} \sim \mathcal{N}(0,\sigma_j).

The contribution to the mean of the observed series obeys

y_{t} = B \alpha_t + \epsilon_t.

plus observation error. Identifiability constraints imply that the observation coefficients B form a rectangular lower triangular matrix with diagonal 1.0.

If there are m time series and p factors, then B has m rows and p columns. Having B be lower triangular means that the first factor affects all series. The second affects all but the first, the third excludes the first two, etc.

Usage

  AddSharedLocalLevel(
     state.specification,
     response,
     nfactors,
     coefficient.prior = NULL,
     initial.state.prior = NULL,
     timestamps = NULL,
     series.id = NULL,
     sdy,
     ...)

Arguments

Value

Returns a list with the elements necessary to specify a local linear trend state model.

Author(s)

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

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

Static Intercept State Component

Description

Adds a static intercept term to a state space model. If the model includes a traditional trend component (e.g. local level, local linear trend, etc) then a separate intercept is not needed (and will probably cause trouble, as it will be confounded with the initial state of the trend model). However, if there is no trend, or the trend is an AR process centered around zero, then adding a static intercept will shift the center to a data-determined value.

Usage

AddStaticIntercept(
    state.specification,
    y,
    initial.state.prior = NormalPrior(y[1], sd(y, na.rm = TRUE)))

Arguments

Value

Returns a list with the information required to specify the state component. If initial.state.prior is specified then y is unused.

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

Robust local linear trend

Description

Add a local level model to a state specification. The local linear trend model assumes that both the mean and the slope of the trend follow random walks. The equation for the mean is

\mu_{t+1} = \mu_t + \delta_t + \epsilon_t \qquad \epsilon_t \sim \mathcal{T}_{\nu_\mu}(0, \sigma_\mu).

The equation for the slope is

\delta_{t+1} = \delta_t + \eta_t \qquad \eta_t \sim \mathcal{T}_{\nu_\delta}(0, \sigma_\delta).

Independent prior distributions are assumed on the level standard deviation, \sigma_\mu the slope standard deviation \sigma_\delta, the level tail thickness \nu_\mu, and the slope tail thickness \nu_\delta.

Usage

  AddStudentLocalLinearTrend(
     state.specification = NULL,
     y,
     save.weights = FALSE,
     level.sigma.prior = NULL,
     level.nu.prior = NULL,
     slope.sigma.prior = NULL,
     slope.nu.prior = NULL,
     initial.level.prior = NULL,
     initial.slope.prior = NULL,
     sdy,
     initial.y)

Arguments

Value

Returns a list with the elements necessary to specify a local linear trend state model.

Author(s)

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

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

  data(rsxfs)
  ss <- AddStudentLocalLinearTrend(list(), rsxfs)
  model <- bsts(rsxfs, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

Trigonometric Seasonal State Component

Description

Add a trigonometric seasonal model to a state specification.

Usage

  AddTrig(
     state.specification = NULL,
     y,
     period,
     frequencies,
     sigma.prior = NULL,
     initial.state.prior = NULL,
     sdy = sd(y, na.rm = TRUE),
     method = c("harmonic", "direct"))

Arguments

Details

Harmonic Method

Each frequency lambda_j = 2\pi j / S where S is the period (number of time points in a full cycle) is associated with two time-varying random components: \gamma_{jt}, and gamma^*_{jt}. They evolve through time as

% \gamma_{j, t + 1} = \gamma_{jt} \cos(\lambda_j) + \gamma^*_{j, t} % \sin(\lambda_j) + \epsilon_{0t}

% \gamma^*_{j, t + 1} = \gamma^*[j, t] \cos(\lambda_j) - \gamma_{jt} * \sin(\lambda_j) + \epsilon_1

where \epsilon_0 and \epsilon_1 are independent with the same variance. This is the real-valued version of a harmonic function: \gamma \exp(i\theta).

The transition matrix multiplies the function by \exp(i \lambda_j, so that after 't' steps the harmonic's value is \gamma \exp(i \lambda_j t).

The model dynamics allows gamma to drift over time in a random walk.

The state of the model is (\gamma_{jt}, \gamma^*_{jt}), for j = 1, ... number of frequencies.

The state transition matrix is a block diagonal matrix, where block 'j' is

\cos(\lambda_j) \sin(\lambda_j)

-\sin(\lambda_j) \cos(\lambda_j)

The error variance matrix is sigma^2 * I. There is a common sigma^2 parameter shared by all frequencies.

The model is full rank, so the state error expander matrix R_t is the identity.

The observation_matrix is (1, 0, 1, 0, ...), where the 1's pick out the 'real' part of the state contributions.

Direct Method

Under the 'direct' method the trig component adds a collection of sine and cosine terms with randomly varying coefficients to the state model. The coefficients are the states, while the sine and cosine values are part of the "observation matrix".

This state component adds the sum of its terms to the observation equation.

y_t = \sum_j \beta_{jt} sin(f_j t) + \gamma_{jt} cos(f_j t)

The evolution equation is that each of the sinusoid coefficients follows a random walk with standard deviation sigma[j].

\beta_{jt} = \beta_{jt-1} + N(0, sigma_{sj}^2) \gamma_{jt} = \gamma_{j-1} + N(0, sigma_{cj}^2)

The direct method is generally inferior to the harmonic method. It may be removed in the future.

Value

Returns a list with the elements necessary to specify a seasonal state model.

Author(s)

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

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 MvnPrior

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddTrig(ss, y, period = 12, frequencies = 1:3)
  model <- bsts(y, state.specification = ss, niter = 200)
  plot(model)

  ## The "harmonic" method is much more stable than the "direct" method.
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddTrig(ss, y, period = 12, frequencies = 1:3, method = "direct")
  model2 <- bsts(y, state.specification = ss, niter = 200)
  plot(model2)

Aggregate a fine time series to a coarse summary

Description

Aggregate measurements from a fine scaled time series into a coarse time series. This is similar to functions from the xts package, but it can handle aggregation from weeks to months.

Usage

   AggregateTimeSeries(fine.series,
                       contains.end,
                       membership.fraction,
                       trim.left = any(membership.fraction < 1),
                       trim.right = NULL,
                       byrow = TRUE)

Arguments

Value

A matrix (if fine.series is a matrix) or vector (otherwise) containing the aggregated values of fine.series.

Author(s)

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

Examples

  week.ending <- as.Date(c("2011-11-05",
                           "2011-11-12",
                           "2011-11-19",
                           "2011-11-26",
                           "2011-12-03",
                           "2011-12-10",
                           "2011-12-17",
                           "2011-12-24",
                           "2011-12-31"))
  membership.fraction <- GetFractionOfDaysInInitialMonth(week.ending)
  which.month <- MatchWeekToMonth(week.ending, as.Date("2011-11-01"))
  contains.end <- WeekEndsMonth(week.ending)

  weekly.values <- rnorm(length(week.ending))
  monthly.values <- AggregateTimeSeries(weekly.values, contains.end, membership.fraction)

Aggregate a weekly time series to monthly

Description

Aggregate measurements from a weekly time series into a monthly time series.

Usage

   AggregateWeeksToMonths(weekly.series,
                          membership.fraction = NULL,
                          trim.left = TRUE,
                          trim.right = NULL)

Arguments

Value

A zoo-matrix (if weekly.series is a matrix) or vector (otherwise) containing the aggregated values of weekly.series.

Author(s)

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

See Also

AggregateTimeSeries

Examples

  week.ending <- as.Date(c("2011-11-05",
                           "2011-11-12",
                           "2011-11-19",
                           "2011-11-26",
                           "2011-12-03",
                           "2011-12-10",
                           "2011-12-17",
                           "2011-12-24",
                           "2011-12-31"))

  weekly.values <- zoo(rnorm(length(week.ending)), week.ending)
  monthly.values <- AggregateWeeksToMonths(weekly.values)

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,
          ...)

Arguments

Details

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

\alpha_{t} = \phi_1\alpha_{i, t-1} + \cdots + \phi_p \alpha_{t-p} + \epsilon_{t-1} \qquad \epsilon_t \sim \mathcal{N}(0, \sigma^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.

Value

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

Author(s)

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

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

Examples

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"))

Bayesian Structural Time Series

Description

Uses MCMC to sample from the posterior distribution of a Bayesian structural time series model. This function can be used either with or without contemporaneous predictor variables (in a time series regression).

If predictor variables are present, the regression coefficients are fixed (as opposed to time varying, though time varying coefficients might be added as state component). The predictors and response in the formula are contemporaneous, so if you want lags and differences you need to put them in the predictor matrix yourself.

If no predictor variables are used, then the model is an ordinary state space time series model.

The model allows for several useful extensions beyond standard Bayesian dynamic linear models.

• A spike-and-slab prior is used for the (static) regression component of models that include predictor variables. This is especially useful with large numbers of regressor series.

• Both the spike-and-slab component (for static regressors) and the Kalman filter (for components of time series state) require observations and state variables to be Gaussian. The bsts package allows for non-Gaussian error families in the observation equation (as well as some state components) by using data augmentation to express these families as conditionally Gaussian.

• As of version 0.7.0, bsts supports having multiple observations at the same time point. In this case the basic model is taken to be

  y_{t,j} = Z_t^T \alpha_t + \beta^Tx_{t, j} + \epsilon_{t,j}.

  This is a regression model where all observations with the same time point share a common time series effect.

Usage

bsts(formula,
     state.specification,
     family = c("gaussian", "logit", "poisson", "student"),
     data,
     prior,
     contrasts = NULL,
     na.action = na.pass,
     niter,
     ping = niter / 10,
     model.options = BstsOptions(),
     timestamps = NULL,
     seed = NULL,
     ...)

Arguments

Details

If the model family is logit, then there are two ways one can format the response variable. If the response is 0/1, TRUE/FALSE, or 1/-1, then the response variable can be passed as with any other model family. If the response is a set of counts out of a specified number of trials then it can be passed as a two-column matrix, where the first column contains the counts of successes and the second contains the count of failures.

Likewise, if the model family is Poisson, the response can be passed as a single vector of counts, under the assumption that each observation has unit exposure. If the exposures differ across observations, then the resopnse can be a two column matrix, with the first column containing the event counts and the second containing exposure times.

Value

An object of class bsts which is a list with the following components

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 3-way 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.

Author(s)

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

References

Scott and Varian (2014) "Predicting the Present with Bayesian Structural Time Series", International Journal of Mathematical Modelling and Numerical Optimization. 4–23.

Scott and Varian (2015) "Bayesian Variable Selection for Nowcasting Economic Time Series", Economic Analysis of the Digital Economy, pp 119-135.

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.

See Also

bsts, AddLocalLevel, AddLocalLinearTrend, AddSemilocalLinearTrend, AddSeasonal AddDynamicRegression SpikeSlabPrior, SdPrior.

Examples

  ## Example 1:  Time series (ts) data
  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  par(mfrow = c(1,2))
  plot(model)
  plot(pred)

## Not run: 

  MakePlots <- function(model, ask = TRUE) {
    ## Make all the plots callable by plot.bsts.
    opar <- par(ask = ask)
    on.exit(par(opar))
    plot.types <- c("state", "components", "residuals",
                    "prediction.errors", "forecast.distribution")
    for (plot.type in plot.types) {
      plot(model, plot.type)
    }
    if (model$has.regression) {
      regression.plot.types <- c("coefficients", "predictors", "size")
      for (plot.type in regression.plot.types) {
        plot(model, plot.type)
      }
    }
  }

  ## Example 2: GOOG is the Google stock price, an xts series of daily
  ##            data.
  data(goog)
  ss <- AddSemilocalLinearTrend(list(), goog)
  model <- bsts(goog, state.specification = ss, niter = 500)
  MakePlots(model)

  ## Example 3:  Change GOOG to be zoo, and not xts.
  goog <- zoo(goog, index(goog))
  ss <- AddSemilocalLinearTrend(list(), goog)
  model <- bsts(goog, state.specification = ss, niter = 500)
  MakePlots(model)

  ## Example 4:  Naked numeric data works too
  y <- rnorm(100)
  ss <- AddLocalLinearTrend(list(), y)
  model <- bsts(y, state.specification = ss, niter = 500)
  MakePlots(model)

  ## Example 5:  zoo data with intra-day measurements
  y <- zoo(rnorm(100),
           seq(from = as.POSIXct("2012-01-01 7:00 EST"), len = 100, by = 100))
  ss <- AddLocalLinearTrend(list(), y)
  model <- bsts(y, state.specification = ss, niter = 500)
  MakePlots(model)

\dontrun{
  ## Example 6:  Including regressors
  data(iclaims)
  ss <- AddLocalLinearTrend(list(), initial.claims$iclaimsNSA)
  ss <- AddSeasonal(ss, initial.claims$iclaimsNSA, nseasons = 52)
  model <- bsts(iclaimsNSA ~ ., state.specification = ss, data =
                initial.claims, niter = 1000)
  plot(model)
  plot(model, "components")
  plot(model, "coefficients")
  plot(model, "predictors")
}

## End(Not run)

## Not run: 
  ## Example 7:  Regressors with multiple time stamps.
  number.of.time.points <- 50
  sample.size.per.time.point <- 10
  total.sample.size <- number.of.time.points * sample.size.per.time.point
  sigma.level <- .1
  sigma.obs <- 1

  ## Simulate some fake data with a local level state component.
  trend <- cumsum(rnorm(number.of.time.points, 0, sigma.level))
  predictors <- matrix(rnorm(total.sample.size * 2), ncol = 2)
  colnames(predictors) <- c("X1", "X2")
  coefficients <- c(-10, 10)
  regression <- as.numeric(predictors %*% coefficients)
  y.hat <- rep(trend, sample.size.per.time.point) + regression
  y <- rnorm(length(y.hat), y.hat, sigma.obs)

  ## Create some time stamps, with multiple observations per time stamp.
  first <- as.POSIXct("2013-03-24")
  dates <- seq(from = first, length = number.of.time.points, by = "month")
  timestamps <- rep(dates, sample.size.per.time.point)

  ## Run the model with a local level trend, and an unnecessary seasonal component.
  ss <- AddLocalLevel(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 7)
  model <- bsts(y ~ predictors, ss, niter = 250, timestamps = timestamps,
                seed = 8675309)
  plot(model)
  plot(model, "components")

## End(Not run)

## Example 8: Non-Gaussian data
## Poisson counts of shark attacks in Florida.
data(shark)
logshark <- log1p(shark$Attacks)
ss.level <- AddLocalLevel(list(), y = logshark)
model <- bsts(shark$Attacks, ss.level, niter = 500,
              ping = 250, family = "poisson", seed = 8675309)

## Poisson data can have an 'exposure' as the second column of a
## two-column matrix.
model <- bsts(cbind(shark$Attacks, shark$Population / 1000),
              state.specification = ss.level, niter = 500,
              family = "poisson", ping = 250, seed = 8675309)

Bsts Model Options

Description

Rarely used modeling options for bsts models.

Usage

BstsOptions(save.state.contributions = TRUE,
            save.prediction.errors = TRUE,
            bma.method = c("SSVS", "ODA"),
            oda.options = list(
                fallback.probability = 0.0,
                eigenvalue.fudge.factor = 0.01),
            timeout.seconds = Inf,
            save.full.state = FALSE)

Arguments

Value

The arguments are checked to make sure they have legal types and values, then a list is returned containing the arguments.

Author(s)

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

Compare bsts models

Description

Produce a set of line plots showing the cumulative absolute one step ahead prediction errors for different models. This plot not only shows which model is doing the best job predicting the data, it highlights regions of the data where the predictions are particularly good or bad.

Usage

CompareBstsModels(model.list,
                  burn = SuggestBurn(.1, model.list[[1]]),
                  filename = "",
                  colors = NULL,
                  lwd = 2,
                  xlab = "Time",
                  main = "",
                  grid = TRUE,
                  cutpoint = NULL)

Arguments

Value

Invisibly returns the matrix of cumulative one-step ahead prediction errors (the lines in the top panel of the plot). Each row in the matrix corresponds to a model in model.list.

Author(s)

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

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  trend.only <- bsts(y, ss, niter = 250)

  ss <- AddSeasonal(ss, y, nseasons = 12)
  trend.and.seasonal <- bsts(y, ss, niter = 250)

  CompareBstsModels(list(trend = trend.only,
                         "trend and seasonal" = trend.and.seasonal))

  CompareBstsModels(list(trend = trend.only,
                         "trend and seasonal" = trend.and.seasonal),
                          cutpoint = 100)

Date Range

Description

Returns the first and last dates of the influence window for the given holiday, among the given timestamps.

Usage

  DateRange(holiday, timestamps)

Arguments

Value

Returns a two-column data frame giving the first and last dates of the influence window for the holiday in the period covered by timestamps.

Author(s)

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

Examples

holiday <- NamedHoliday("MemorialDay", days.before = 2, days.after = 2)
timestamps <- seq.Date(from = as.Date("2001-01-01"), by = "day",
   length.out = 365 * 10)
influence <- DateRange(holiday, timestamps)

Descriptive Plots

Description

Plots for describing time series data.

Usage

  DayPlot(y, colors = NULL, ylab = NULL, ...)
  MonthPlot(y, seasonal.identifier = months, colors = NULL, ylab = NULL, ...)
  YearPlot(y, colors = NULL, ylab = NULL, ylim = NULL, legend = TRUE, ...)

Arguments

Details

DayPlot and MonthPlot plot the time series one season at a time, on the same set of axes. The intent is to use DayPlot for daily data and MonthPlot for monthly or quarterly data.

YearPlot plots each year of the time series as a separate line on the same set of axes.

Both sets of plots help visualize seasonal patterns.

Value

Returns invisible{NULL}.

See Also

monthplot is a base R function for plotting time series of type ts.

Examples

## Plot a 'ts' time series.
data(AirPassengers)
par(mfrow = c(1,2))
MonthPlot(AirPassengers)
YearPlot(AirPassengers)

## Plot a 'zoo' time series.
data(turkish)
par(mfrow = c(1,2))
YearPlot(turkish)
DayPlot(turkish)

Diagnostic Plots

Description

Diagnostic plots for distributions of residuals.

Usage

  qqdist(draws, ...)
  AcfDist(draws, lag.max = NULL, xlab = "Lag", ylab = "Autocorrelation", ...)

Arguments

Details

qqdist sorts the columns of draws by their mean, and plots the resulting set of curves against the quantiles of the standard normal distribution. A reference line is added, and the mean of each column of draws is represented by a blue dot. The dots and the line are the transpose of what you get with qqnorm and qqline.

AcfDist plots the posterior distribution of the autocorrelation function using a set of side-by-side boxplots.

Examples

data(AirPassengers)
y <- log(AirPassengers)

ss <- AddLocalLinearTrend(list(), y)
ss <- AddSeasonal(ss, y, nseasons = 12)
model <- bsts(y, ss, niter = 500)

r <- residuals(model)
par(mfrow = c(1,2))
qqdist(r)   ## A bit of departure in the upper tail
AcfDist(r)  

Dynamic intercept regression model

Description

A dynamic intercept regression is a regression model where the intercept term is a state space model. This model differs from bsts in that there can be multiple observations per time point.

Usage

dirm(formula,
     state.specification,
     data,
     prior = NULL,
     contrasts = NULL,
     na.action = na.pass,
     niter,
     ping = niter / 10,
     model.options = DirmModelOptions(),
     timestamps = NULL,
     seed = NULL,
     ...)

Arguments

Details

The fitted model is a regression model with an intercept term given by a structural time series model. This is similar to the model fit by bsts, but it allows for multiple observations per time period.

Currently dirm only supports Gaussian observation errors, but look for that to change in future releases.

Value

An object of class bsts which is a list with the following components

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 3-way 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.

Author(s)

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

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.

George and McCulloch (1997) "Approaches for Bayesian variable selection", Statistica Sinica pp 339–374.

See Also

bsts, AddLocalLevel, AddLocalLinearTrend, AddSemilocalLinearTrend, AddSeasonal AddDynamicRegression SpikeSlabPrior, SdPrior.

Examples

SimulateDirmData <- function(observation.sd = 1, trend.sd = .1,
                             time.dimension = 100, nobs.per.period = 3,
                             xdim = 4) {
  trend <- cumsum(rnorm(time.dimension, 0, trend.sd))
  total.sample.size <- nobs.per.period * time.dimension
  predictors <- matrix(rnorm(total.sample.size * xdim),
    nrow = total.sample.size)
  coefficients <- rnorm(xdim)
  expanded.trend <- rep(trend, each = nobs.per.period)
  response <- expanded.trend + predictors %*% coefficients + rnorm(
    total.sample.size, 0, observation.sd)
  timestamps <- seq.Date(from = as.Date("2008-01-01"),
                         len = time.dimension, by = "day")
  extended.timestamps <- rep(timestamps, each = nobs.per.period)
  return(list(response = response,
    predictors = predictors,
    timestamps = extended.timestamps,
    trend = trend,
    coefficients = coefficients))
}


data <- SimulateDirmData(time.dimension = 20)
ss <- AddLocalLevel(list(), data$response)

# In real life you'd want more than 50 MCMC iterations.
model <- dirm(data$response ~ data$predictors, ss, niter = 50,
  timestamps = data$timestamps)

Specify Options for a Dynamic Intercept Regression Model

Description

Specify modeling options for a dynamic intercept regression model.

Usage

DirmModelOptions(timeout.seconds = Inf,
                 high.dimensional.threshold.factor = 1.0)

Arguments

Value

An object of class DirmModelOptions, which is simply a list containing values of the function arguments.

The value of using this function instead of making a list "by hand" is that argument types are properly checked, and list names are sure to be correct.

Intervals between dates

Description

Estimate the time scale used in time series data.

Usage

EstimateTimeScale(dates)

Arguments

Value

A character string. Either "daily", "weekly", "yearly", "monthly", "quarterly", or "other". The value is determined based on counting the number of days between successive observations in dates.

Author(s)

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

Examples

weekly.data <- as.Date(c("2011-10-01",
                         "2011-10-08",
                         "2011-10-15",
                         "2011-10-22",
                         "2011-10-29",
                         "2011-11-05"))

EstimateTimeScale(weekly.data) # "weekly"

almost.weekly.data <- as.Date(c("2011-10-01",
                                "2011-10-08",
                                "2011-10-15",
                                "2011-10-22",
                                "2011-10-29",
                                "2011-11-06"))  # last day is one later

EstimateTimeScale(weekly.data) # "other"

Extends a vector of dates to a given length

Description

Pads a vector of dates to a specified length.

Usage

   ExtendTime(dates, number.of.periods, dt = NULL)

Arguments

Value

If number.of.periods is longer than length(dates), then dates will be padded to the desired length. Extra dates are added at time intervals matching the average interval in dates. Thus they may not be

Author(s)

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

See Also

bsts.mixed.

Examples

  origin.month <- as.Date("2011-09-01")
  week.ending <- as.Date(c("2011-10-01",   ## 1
                           "2011-10-08",   ## 2
                           "2011-12-03",   ## 3
                           "2011-12-31"))  ## 4
  MatchWeekToMonth(week.ending, origin.month) == 1:4

Checking for Regularity

Description

Tools for checking if a series of timestamps is 'regular' meaning that it has no duplicates, and no gaps. Checking for regularity can be tricky. For example, if you have monthly observations with Date or POSIXt timestamps then gaps between timestamps can be 28, 29, 30, or 31 days, but the series is still "regular".

Usage

  NoDuplicates(timestamps)
  NoGaps(timestamps)
  IsRegular(timestamps)

  HasDuplicateTimestamps(bsts.object)

Arguments

Value

All four functions return scalar logical values. NoDuplicates returns TRUE if all elements of timestamps are unique.

NoGaps examines the smallest nonzero gap between time points. As long as no gaps between time points are more than twice as wide as the smallest gap, it returns TRUE, indicating that there are no missing timestamps. Otherwise it returns FALSE.

IsRegular returns TRUE if NoDuplicates and NoGaps both return TRUE.

HasDuplicateTimestamps returns FALSE if the data used to fit bsts.model either has NULL timestamps, or if the timestamps contain no duplicate values.

Author(s)

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

Examples

  first <- as.POSIXct("2015-04-19 08:00:04")
  monthly <- seq(from = first, length.out = 24, by = "month")
  IsRegular(monthly) ## TRUE

  skip.one <- monthly[-8]
  IsRegular(skip.one) ## FALSE

  has.duplicates <- monthly
  has.duplicates[1] <- has.duplicates[2]
  IsRegular(has.duplicates) ## FALSE

Gross Domestic Product for 57 Countries

Description

Annual gross domestic product for 57 countries, as produced by the OECD.

Fields:

• LOCATION: Three letter country code.

• MEASURE: MLN_USD signifies a total GDP number in millions of US dollars. USD_CAP is per capita GDP in US dollars.

• TIME: The year of the measurement.

• Value: The measured value.

• Flag.Codes: P for provisional data, B for a break in the series, and E for an estimated value.

Usage

data(gdp)

Format

data frame

Source

OECD website: See https://data.oecd.org/gdp/gross-domestic-product-gdp.htm

Create a Geometric Sequence

Description

Create a geometric sequence.

Usage

  GeometricSequence(length, initial.value = 1, discount.factor = .5)

Arguments

Value

A numeric vector containing the desired sequence.

Author(s)

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

Examples

GeometricSequence(4, .8, .6)
# [1] 0.8000 0.4800 0.2880 0.1728

GeometricSequence(5, 2, 3)
# [1]   2   6  18  54 162

## Not run: 
GeometricSequence(0, -1, -2)
# Error: length > 0 is not TRUE

## End(Not run)

Compute membership fractions

Description

Returns the fraction of days in a week that occur in the ear

Usage

   GetFractionOfDaysInInitialMonth(week.ending)
   GetFractionOfDaysInInitialQuarter(week.ending)

Arguments

Value

Returns a numeric vector of the same length as week.ending. Each entry gives the fraction of days in the week that occur in the coarse time interval (month or quarter) containing the start of the week (i.e the date 6 days before).

Author(s)

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

See Also

bsts.mixed.

Examples

    dates <- as.Date(c("2003-03-31",
                       "2003-04-01",
                       "2003-04-02",
                       "2003-04-03",
                       "2003-04-04",
                       "2003-04-05",
                       "2003-04-06",
                       "2003-04-07"))
    fraction <- GetFractionOfDaysInInitialMonth(dates)
    fraction == c(1, 6/7, 5/7, 4/7, 3/7, 2/7, 1/7, 1)

Google stock price

Description

Daily closing price of Google stock.

Usage

data(goog)

Format

xts time series

Source

The Internets

Specifying Holidays

Description

Specify holidays for use with holiday state models.

Usage

FixedDateHoliday(holiday.name,
                 month = base::month.name,
                 day,
                 days.before = 1,
                 days.after = 1)

NthWeekdayInMonthHoliday(holiday.name,
                         month = base::month.name,
                         day.of.week = weekday.names,
                         week.number = 1,
                         days.before = 1,
                         days.after = 1)

LastWeekdayInMonthHoliday(holiday.name,
                          month = base::month.name,
                          day.of.week = weekday.names,
                          days.before = 1,
                          days.after = 1)

NamedHoliday(holiday.name = named.holidays,
             days.before = 1,
             days.after = 1)

DateRangeHoliday(holiday.name,
                 start.date,
                 end.date)

Arguments

Value

Each function returns a list containing the information from the function arguments, formatted as expected by the underlying C++ code. State models that focus on holidays, such as AddRandomWalkHoliday, AddRegressionHoliday, and AddHierarchicalRegressionHoliday, will expect one or more holiday objects as arguments.

• FixedDateHoliday describes a holiday that occurs on the same date each year, like US independence day (July 4).

• NthWeekdayInMonthHoliday describes a holiday that occurs a particular weekday of a particular week of a particular month. For example, US Labor Day is the first Monday in September.

• LastWeekdayInMonthHoliday describes a holiday that occurs on the last instance of a particular weekday in a particular month. For example, US Memorial Day is the last Monday in May.

• DateRangeHoliday describes an irregular holiday that might not follow a particular pattern. You can handle this type of holiday by manually specifying a range of dates for each instance of the holiday in your data set. NOTE: If you plan on using the model to forecast, be sure to include date ranges in the forecast period as well as the period covered by the training data.

• NamedHoliday is a convenience class for describing several important holidays in the US.

Author(s)

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

See Also

AddRandomWalkHoliday, AddRegressionHoliday, AddHierarchicalRegressionHoliday

Examples

july4 <- FixedDateHoliday("July4", "July", 4)
memorial.day <- LastWeekdayInMonthHoliday("MemorialDay", "May", "Monday")
labor.day <- NthWeekdayInMonthHoliday("LaborDay", "September", "Monday", 1)
another.way.to.get.memorial.day <- NamedHoliday("MemorialDay")
easter <- NamedHoliday("Easter")
winter.olympics <- DateRangeHoliday("WinterOlympicsSince2000",
                     start = as.Date(c("2002-02-08",
                                        "2006-02-10",
                                        "2010-02-12",
                                        "2014-02-07",
                                        "2018-02-07")),
                     end = as.Date(c("2002-02-24",
                                     "2006-02-26",
                                     "2010-02-28",
                                     "2014-02-23",
                                     "2018-02-25")))

Initial Claims Data

Description

A weekly time series of US initial claims for unemployment. The first column contains the initial claims numbers from FRED. The others contain a measure of the relative popularity of various search queries identified by Google Correlate.

Usage

data(iclaims)

Format

zoo time series

Source

FRED. http://research.stlouisfed.org/fred2/series/ICNSA,
Google correlate. http://www.google.com/trends/correlate

See Also

bsts

Examples

data(iclaims)
plot(initial.claims)

Find the last day in a month

Description

Finds the last day in the month containing a specefied date.

Usage

   LastDayInMonth(dates)

Arguments

Value

A vector of class Date where each entry is the last day in the month containing the corresponding entry in dates.

Author(s)

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

Examples

 inputs <- as.Date(c("2007-01-01",
                     "2007-01-31",
                     "2008-02-01",
                     "2008-02-29",
                     "2008-03-14",
                     "2008-12-01",
                     "2008-12-31"))
 expected.outputs <- as.Date(c("2007-01-31",
                              "2007-01-31",
                              "2008-02-29",
                              "2008-02-29",
                              "2008-03-31",
                              "2008-12-31",
                              "2008-12-31"))
 LastDayInMonth(inputs) == expected.outputs

Find the month containing a week

Description

Returns the index of a month, in a sequence of months, that contains a given week.

Usage

   MatchWeekToMonth(week.ending, origin.month)

Arguments

Value

The index of the month matching the month containing the first day in week.ending. The origin is month 1. It is the caller's responsibility to ensure that these indices correspond to legal values in a particular vector of months.

Author(s)

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

See Also

bsts.mixed.

Examples

  origin.month <- as.Date("2011-09-01")
  week.ending <- as.Date(c("2011-10-01",   ## 1
                           "2011-10-08",   ## 2
                           "2011-12-03",   ## 3
                           "2011-12-31"))  ## 4
  MatchWeekToMonth(week.ending, origin.month) == 1:4

Maximum Window Width for a Holiday

Description

The maximum width of a holiday's influence window

Usage

## Default S3 method:
MaxWindowWidth(holiday, ...)
## S3 method for class 'DateRangeHoliday'
MaxWindowWidth(holiday, ...)

Arguments

Value

Returns the number of days in a holiday's influence window.

Author(s)

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

See Also

Holiday. AddRegressionHoliday. AddRandomWalkHoliday. AddHierarchicalRegressionHoliday.

Examples

easter <- NamedHoliday("Easter", days.before = 2, days.after = 1)
if (MaxWindowWidth(easter) == 4) {
  print("That's the right answer!\n")
}

## This holiday lasts two days longer in 2005 than in 2004.
may18 <- DateRangeHoliday("May18",
     start = as.Date(c("2004-05-17",
                       "2005-05-16")),
     end   = as.Date(c("2004-05-19",
                       "2005-05-20")))

if (MaxWindowWidth(may18) == 5) {
   print("Right again!\n")
}

Multivariate Bayesian Structural Time Series

Description

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.

Usage

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,
      ...)

Arguments

Value

An object of class mbsts which is a list with the following components

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 3-way 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.

Author(s)

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

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.

George and McCulloch (1997) "Approaches for Bayesian variable selection", Statistica Sinica pp 339–374.

See Also

bsts, AddLocalLevel, AddLocalLinearTrend, AddSemilocalLinearTrend, AddSeasonal AddDynamicRegression SpikeSlabPrior, SdPrior.

Examples

## Not run: 

# This example takes 12s on Windows, which is longer than CRAN's 10s
# limit.  Marking code as 'dontrun' to prevent CRAN auto checks from
# timing out.

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)
# next line closes 'dontrun'

## End(Not run)
# next line closes 'examples'

Models for mixed frequency time series

Description

Fit a structured time series to mixed frequncy data.

Usage

   bsts.mixed(target.series,
              predictors,
              which.coarse.interval,
              membership.fraction,
              contains.end,
              state.specification,
              regression.prior,
              niter,
              ping = niter / 10,
              seed = NULL,
              truth = NULL,
              ...)

Arguments

Value

An object of class bsts.mixed, which is a list with the following elements. Many of these are arrays, in which case the first index of the array corresponds to the MCMC iteration number.

The returned object also contains MCMC draws for the parameters of the state models supplied as part of state.specification, relevant information passed to the function call, and other supplemental information.

Author(s)

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

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, AddLocalLevel, AddLocalLinearTrend, AddSemilocalLinearTrend, SpikeSlabPrior, SdPrior.

Examples

## Not run: 
  data <- SimulateFakeMixedFrequencyData(nweeks = 104, xdim = 20)

  ## Setting an upper limit on the standard deviations can help keep the
  ## MCMC from flying off to infinity.
  sd.limit <- sd(data$coarse.target)
  state.specification <-
       AddLocalLinearTrend(list(),
                     data$coarse.target,
                     level.sigma.prior = SdPrior(1.0, 5, upper.limit = sd.limit),
                     slope.sigma.prior = SdPrior(.5, 5, upper.limit = sd.limit))
  weeks <- index(data$predictor)
  months <- index(data$coarse.target)
  which.month <- MatchWeekToMonth(weeks, months[1])
  membership.fraction <- GetFractionOfDaysInInitialMonth(weeks)
  contains.end <- WeekEndsMonth(weeks)

  model <- bsts.mixed(target.series = data$coarse.target,
                      predictors = data$predictors,
                      membership.fraction = membership.fraction,
                      contains.end = contains.end,
                      which.coarse = which.month,
                      state.specification = state.specification,
                      niter = 500,
                      expected.r2 = .999,
                      prior.df = 1)

  plot(model, "state")
  plot(model, "components")

## End(Not run)

Elapsed time in months

Description

The (integer) number of months between dates.

Usage

   MonthDistance(dates, origin)

Arguments

Value

Returns a numeric vector giving the integer number of months that have elapsed between origin and each element in dates. The daily component of each date is ignored, so two dates that are in the same month will have the same measured distance. Distances are signed, so months that occur before origin will have negative values.

Author(s)

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

Examples

  dates <- as.Date(c("2008-04-17",
                     "2008-05-01",
                     "2008-05-31",
                     "2008-06-01"))
  origin <- as.Date("2008-05-15")
  MonthDistance(dates, origin) ==  c(-1, 0, 0, 1)

Holidays Recognized by Name

Description

A character vector listing the names of pre-specified holidays.

Usage

named.holidays

Value

"NewYearsDay" "SuperBowlSunday" "MartinLutherKingDay" "PresidentsDay" "ValentinesDay" "SaintPatricksDay" "USDaylightSavingsTimeBegins" "USDaylightSavingsTimeEnds" "EasterSunday" "USMothersDay" "IndependenceDay" "LaborDay" "ColumbusDay" "Halloween" "Thanksgiving" "MemorialDay" "VeteransDay" "Christmas"

New home sales and Google trends

Description

The first column, HSN1FNSA is a time series of new home sales in the US, obtained from the FRED online data base. The series has been manually deseasonalized. The remaining columns contain search terms from Google trends (obtained from http://trends.google.com/correlate). These show the relative popularity of each search term among all serach terms typed into Google. All series in this data set have been standardized by subtracting off their mean and dividing by their standard deviation.

Usage

data(new.home.sales)

Format

zoo time series

Source

FRED and trends.google.com

Prediction Errors

Description

Computes the one-step-ahead prediction errors for a bsts model.

Usage

bsts.prediction.errors(bsts.object,
                       cutpoints = NULL,
                       burn = SuggestBurn(.1, bsts.object),
                       standardize = FALSE)

Arguments

Details

Returns the posterior distribution of the one-step-ahead prediction errors from the bsts.object. The errors are computing using the Kalman filter, and are of two types.

Purely in-sample errors are computed as a by-product of the Kalman filter as a result of fitting the model. These are stored in the bsts.object assuming the save.prediction.errors option is TRUE, which is the default (See BstsOptions). The in-sample errors are 'in-sample' in the sense that the parameter values used to run the Kalman filter are drawn from their posterior distribution given complete data. Conditional on the parameters in that MCMC iteration, each 'error' is the difference between the observed y[t] and its expectation given data to t-1.

Purely out-of-sample errors can be computed by specifying the 'cutpoints' argument. If cutpoints are supplied then a separate MCMC is run using just data up to the cutpoint. The Kalman filter is then run on the remaining data, again finding the difference between y[t] and its expectation given data to t-1, but conditional on parameters estimated using data up to the cutpoint.

Value

A matrix of draws of the one-step-ahead prediction errors. Rows of the matrix correspond to MCMC draws. Columns correspond to time.

Author(s)

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

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, AddLocalLevel, AddLocalLinearTrend, AddSemilocalLinearTrend, SpikeSlabPrior, SdPrior.

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)

## Not run: 
  model <- bsts(y, state.specification = ss, niter = 500)

## End(Not run)

  errors <- bsts.prediction.errors(model, burn = 100)
  PlotDynamicDistribution(errors$in.sample)

  ## Compute out of sample prediction errors beyond times 80 and 120.
  errors <- bsts.prediction.errors(model, cutpoints = c(80, 120))
  standardized.errors <- bsts.prediction.errors(
      model, cutpoints = c(80, 120), standardize = TRUE)
  plot(model, "prediction.errors", cutpoints = c(80, 120))
  str(errors)     ## three matrices, with 400 ( = 500 - 100) rows
                  ## and length(y) columns

Plotting functions for Bayesian structural time series

Description

Functions to plot the results of a model fit using bsts.

Usage

  ## S3 method for class 'bsts'
plot(x, y = c("state", "components", "residuals",
              "coefficients", "prediction.errors",
              "forecast.distribution",
              "predictors", "size", "dynamic", "seasonal", "monthly",
              "help"),
      ...)

  PlotBstsCoefficients(bsts.object, burn = SuggestBurn(.1, bsts.object),
                       inclusion.threshold = 0, number.of.variables = NULL, ...)

  PlotBstsComponents(bsts.object,
                     burn = SuggestBurn(.1, bsts.object),
                     time,
                     same.scale = TRUE,
                     layout = c("square", "horizontal", "vertical"),
                     style = c("dynamic", "boxplot"),
                     ylim = NULL,
                     components = 1:length(bsts.object$state.specification),
                     ...)

  PlotDynamicRegression(bsts.object,
                        burn = SuggestBurn(.1, bsts.object),
                        time = NULL,
                        same.scale = FALSE,
                        style = c("dynamic", "boxplot"),
                        layout = c("square", "horizontal", "vertical"),
                        ylim = NULL,
                        zero.width = 2,
                        zero.color = "green",
                        ...)

  PlotBstsState(bsts.object, burn = SuggestBurn(.1, bsts.object),
                       time, show.actuals = TRUE,
                       style = c("dynamic", "boxplot"),
                       scale = c("linear", "mean"),
                       ylim = NULL,
                       ...)

  PlotBstsResiduals(bsts.object, burn = SuggestBurn(.1, bsts.object),
                       time, style = c("dynamic", "boxplot"), means =
                       TRUE, ...)

  PlotBstsPredictionErrors(bsts.object, cutpoints = NULL,
                           burn = SuggestBurn(.1, bsts.object),
                           style = c("dynamic", "boxplot"),
                           xlab = "Time", ylab = "", main = "",
                           ...)

  PlotBstsForecastDistribution(bsts.object, cutpoints = NULL,
                               burn = SuggestBurn(.1, bsts.object),
                               style = c("dynamic", "boxplot"),
                               xlab = "Time",
                               ylab = "",
                               main = "",
                               show.actuals = TRUE,
                               col.actuals = "blue",
                               ...)

  PlotBstsSize(bsts.object, burn = SuggestBurn(.1, bsts.object), style =
                       c("histogram", "ts"), ...)

  PlotSeasonalEffect(bsts.object, nseasons = 7, season.duration = 1,
                     same.scale = TRUE, ylim = NULL, get.season.name = NULL,
                     burn = SuggestBurn(.1, bsts.object),  ...)

  PlotMonthlyAnnualCycle(bsts.object, ylim = NULL, same.scale = TRUE,
                     burn = SuggestBurn(.1, bsts.object),  ...)

Arguments

Details

PlotBstsState, PlotBstsComponents, and PlotBstsResiduals all produce dynamic distribution plots. PlotBstsState plots the aggregate state contribution (including regression effects) to the mean, while PlotBstsComponents plots the contribution of each state component. PlotBstsResiduals plots the posterior distribution of the residuals given complete data (i.e. looking forward and backward in time). PlotBstsPredictionErrors plots filtering errors (i.e. the one-step-ahead prediction errors given data up to the previous time point). PlotBstsForecastDistribution plots the one-step-ahead forecasts instead of the prediction errors.

PlotBstsCoefficients creates a significance plot for the predictors used in the state space regression model. It is obviously not useful for models with no regressors.

PlotBstsSize plots the distribution of the number of predictors included in the model.

PlotSeasonalEffect generates an array of plots showing how the distibution of the seasonal effect changes, for each season, for models that include a seasonal state component.

PlotMonthlyAnnualCycle produces an array of plots much like PlotSeasonalEffect, for models that include a MonthlyAnnualCycle state component.

Value

These functions are called for their side effect, which is to produce a plot on the current graphics device.

PlotBstsState invisibly returns the state object being plotted.

See Also

bsts PlotDynamicDistribution plot.lm.spike

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  plot(model, burn = 100)
  plot(model, "residuals", burn = 100)
  plot(model, "components", burn = 100)
  plot(model, "forecast.distribution", burn = 100)

Plotting functions for mixed frequency Bayesian structural time series

Description

Functions for plotting the output of a mixed frequency time series regression.

Usage

  ## S3 method for class 'bsts.mixed'
plot(x,
                            y = c("state", "components",
                                  "coefficients", "predictors", "size"),
                            ...)

   PlotBstsMixedState(bsts.mixed.object,
                      burn = SuggestBurn(.1, bsts.mixed.object),
                      time = NULL,
                      fine.scale = FALSE,
                      style = c("dynamic", "boxplot"),
                      trim.left = NULL,
                      trim.right = NULL,
                      ...)

   PlotBstsMixedComponents(bsts.mixed.object,
                           burn = SuggestBurn(.1, bsts.mixed.object),
                           time = NULL,
                           same.scale = TRUE,
                           fine.scale = FALSE,
                           style = c("dynamic", "boxplot"),
                           layout = c("square", "horizontal", "vertical"),
                           ylim = NULL,
                           trim.left = NULL,
                           trim.right = NULL,
                           ...)

Arguments

Details

PlotBstsMixedState plots the aggregate state contribution (including regression effects) to the mean, while PlotBstsComponents plots the contribution of each state component separately. PlotBstsCoefficients creates a significance plot for the predictors used in the state space regression model.

Value

These functions are called for their side effect, which is to produce a plot on the current graphics device.

See Also

bsts.mixed PlotDynamicDistribution plot.lm.spike PlotBstsSize

Examples

## Not run: 
## This example is flaky and needs to be fixed
  data <- SimulateFakeMixedFrequencyData(nweeks = 104, xdim = 20)
  state.specification <- AddLocalLinearTrend(list(), data$coarse.target)
  weeks <- index(data$predictor)
  months <- index(data$coarse.target)
  which.month <- MatchWeekToMonth(weeks, months[1])
  membership.fraction <- GetFractionOfDaysInInitialMonth(weeks)
  contains.end <- WeekEndsMonth(weeks)

  model <- bsts.mixed(target.series = data$coarse.target,
                      predictors = data$predictors,
                      membership.fraction = membership.fraction,
                      contains.end = contains.end,
                      which.coarse = which.month,
                      state.specification = state.specification,
                      niter = 500)

  plot(model, "state")
  plot(model, "components")

## End(Not run)

Plot predictions from Bayesian structural time series

Description

Plot the posterior predictive distribution from a bsts prediction object.

Usage

  ## S3 method for class 'bsts.prediction'
plot(x,
     y = NULL,
     burn = 0,
     plot.original = TRUE,
     median.color = "blue",
     median.type = 1,
     median.width = 3,
     interval.quantiles = c(.025, .975),
     interval.color = "green",
     interval.type = 2,
     interval.width = 2,
     style = c("dynamic", "boxplot"),
     ylim = NULL,
     ...)

Arguments

Details

Plots the posterior predictive distribution described by x using a dynamic distribution plot generated by PlotDynamicDistribution. Overlays the posterior median and 95% prediction limits for the predictive distribution.

Value

Returns NULL.

See Also

bsts PlotDynamicDistribution plot.lm.spike

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 500)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

Plot the most likely predictors

Description

Creates a time series plot showing the most likely predictors of a time series used to fit a bsts object.

Usage

  PlotBstsPredictors(bsts.object,
                     burn = SuggestBurn(.1, bsts.object),
                     inclusion.threshold = .1,
                     ylim = NULL,
                     flip.signs = TRUE,
                     show.legend = TRUE,
                     grayscale = TRUE,
                     short.names = TRUE,
                     ...)

Arguments

See Also

bsts PlotDynamicDistribution plot.lm.spike

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  lag.y <- c(NA, head(y, -1))
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  ## Call bsts with na.action = na.omit to omit the leading NA in lag.y
  model <- bsts(y ~ lag.y, state.specification = ss, niter = 500,
                na.action = na.omit)
  plot(model, "predictors")

Plot Holiday Effects

Description

Plot the estimated effect of the given holiday.

Usage

  PlotHoliday(holiday, model, show.raw.data = TRUE, ylim = NULL, ...)

Arguments

Value

Returns invisible{NULL}.

See Also

bsts AddRandomWalkHoliday

Examples

  trend <- cumsum(rnorm(730, 0, .1))
  dates <- seq.Date(from = as.Date("2014-01-01"), length = length(trend),
    by = "day")
  y <- zoo(trend + rnorm(length(trend), 0, .2), dates)

  AddHolidayEffect <- function(y, dates, effect) {
    ## Adds a holiday effect to simulated data.
    ## Args:
    ##   y: A zoo time series, with Dates for indices.
    ##   dates: The dates of the holidays.
    ##   effect: A vector of holiday effects of odd length.  The central effect is
    ##     the main holiday, with a symmetric influence window on either side.
    ## Returns:
    ##   y, with the holiday effects added.
    time <- dates - (length(effect) - 1) / 2
    for (i in 1:length(effect)) {
      y[time] <- y[time] + effect[i]
      time <- time + 1
    }
    return(y)
  }

  ## Define some holidays.
  memorial.day <- NamedHoliday("MemorialDay")
  memorial.day.effect <- c(.3, 3, .5)
  memorial.day.dates <- as.Date(c("2014-05-26", "2015-05-25"))
  y <- AddHolidayEffect(y, memorial.day.dates, memorial.day.effect)
  
  presidents.day <- NamedHoliday("PresidentsDay")
  presidents.day.effect <- c(.5, 2, .25)
  presidents.day.dates <- as.Date(c("2014-02-17", "2015-02-16"))
  y <- AddHolidayEffect(y, presidents.day.dates, presidents.day.effect)
  
  labor.day <- NamedHoliday("LaborDay")
  labor.day.effect <- c(1, 2, 1)
  labor.day.dates <- as.Date(c("2014-09-01", "2015-09-07"))
  y <- AddHolidayEffect(y, labor.day.dates, labor.day.effect)
  
  ## The holidays can be in any order.
  holiday.list <- list(memorial.day, labor.day, presidents.day)
  number.of.holidays <- length(holiday.list)
  
  ## In a real example you'd want more than 100 MCMC iterations.
  niter <- 100
  ss <- AddLocalLevel(list(), y)
  ss <- AddRegressionHoliday(ss, y, holiday.list = holiday.list)
  model <- bsts(y, state.specification = ss, niter = niter)

  PlotHoliday(memorial.day, model)

Plotting Functions for Multivariate Bayesian Structural Time Series

Description

Functions to plot the results of a model fit using mbsts.

Usage

## S3 method for class 'mbsts'
plot(x, y = c("means", "help"), ...)

PlotMbstsSeriesMeans(mbsts.object,
                     series.id = NULL,
                     same.scale = TRUE,
                     burn = SuggestBurn(.1, mbsts.object),
                     time,
                     show.actuals = TRUE,
                     ylim = NULL,
                     gap = 0,
                     cex.actuals = 0.2,
                     ...)

Arguments

See Also

plot.bsts

Plot Multivariate Bsts Predictions

Description

Plot the posterior predictive distribution from an mbsts prediction object.

Usage

  ## S3 method for class 'mbsts.prediction'
plot(x,
     y = NULL,
     burn = 0,
     plot.original = TRUE,
     median.color = "blue",
     median.type = 1,
     median.width = 3,
     interval.quantiles = c(.025, .975),
     interval.color = "green",
     interval.type = 2,
     interval.width = 2,
     style = c("dynamic", "boxplot"),
     ylim = NULL,
     series.id = NULL,
     same.scale = TRUE,
     gap = 0,
     ...)

Arguments

Details

Plots the posterior predictive distribution described by x using a dynamic distribution plot generated by PlotDynamicDistribution. Overlays the posterior median and 95% prediction limits for the predictive distribution.

Value

Returns NULL.

Prediction for Bayesian Structural Time Series

Description

Generate draws from the posterior predictive distribution of a bsts object.

Usage

## S3 method for class 'bsts'
predict(object,
        horizon = 1,
        newdata = NULL,
        timestamps = NULL,
        burn = SuggestBurn(.1, object),
        na.action = na.exclude,
        olddata = NULL,
        olddata.timestamps = NULL,
        trials.or.exposure = 1,
        quantiles = c(.025, .975),
        seed = NULL,
        ...)

Arguments

Details

Samples from the posterior distribution of a Bayesian structural time series model. This function can be used either with or without contemporaneous predictor variables (in a time series regression).

If predictor variables are present, the regression coefficients are fixed (as opposed to time varying, though time varying coefficients might be added as state component). The predictors and response in the formula are contemporaneous, so if you want lags and differences you need to put them in the predictor matrix yourself.

If no predictor variables are used, then the model is an ordinary state space time series model.

Value

Returns an object of class bsts.prediction, which is a list with the following components.

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. AddLocalLevel. AddLocalLinearTrend. AddSemilocalLinearTrend.

Examples

# The number of MCMC draws in the following examples is artificially low.


  ## Making predictions when there is no regression component.
  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 250)
  pred <- predict(model, horizon = 12, burn = 100)
  plot(pred)

  ## An example using the olddata argument.
  full.pred <- pred
  training <- window(y, end = c(1959, 12))
  model <- bsts(training, state.specification = ss, niter = 250)
  ## Predict the next 12 months.
  pred <- predict(model, horizon = 12)
  ## Compare the predictions to the actual data.
  plot(pred)
  lines(as.numeric(y, col = "red", lty = 2, lwd = 2))

  ## Predict the 12 months of 1961 based on the posterior distribution
  ## of the model fit to data through 1959, but with state filtered
  ## through 1960.
  updated.pred <- predict(model, horizon = 12, olddata = y)
  par(mfrow = c(1, 2))
  plot(full.pred, ylim = c(4, 7))
  plot(updated.pred, ylim = c(4, 7))

  ## Examples including a regression component.
  ##
  data(iclaims)
  training <- initial.claims[1:402, ]
  holdout1 <- initial.claims[403:450, ]
  holdout2 <- initial.claims[451:456, ]

## Not run: 

## This example puts the total run time over 5 seconds, which is a CRAN
## violation.

  ss <- AddLocalLinearTrend(list(), training$iclaimsNSA)
  ss <- AddSeasonal(ss, training$iclaimsNSA, nseasons = 52)
  ## In real life you'd want more iterations...
  model <- bsts(iclaimsNSA ~ ., state.specification = ss, data =
                training, niter = 100)

  ## Predict the holdout set given the training set.
  ## This is really fast, because we can use saved state from the MCMC
  ## algorithm.
  pred.full <- predict(model, newdata = rbind(holdout1, holdout2))

  ## Predict holdout 2, given training and holdout1.
  ## This is much slower because we need to re-filter the 'olddata' before
  ## simulating the predictions.
  pred.update <- predict(model, newdata = holdout2,
    olddata = rbind(training, holdout1))

## End(Not run)

Prediction for Multivariate Bayesian Structural Time Series

Description

Generate draws from the posterior predictive distribution of an mbsts object.

Usage

## S3 method for class 'mbsts'
predict(object,
        horizon = 1,
        newdata = NULL,
        timestamps = NULL,
        burn = SuggestBurn(.1, object),
        na.action = na.exclude,
        quantiles = c(.025, .975),
        seed = NULL,
        ...)

Arguments

Details

The prediction is based off of samples taken from the posterior distribution of a multivariate Bayesian structural time series model.

As an added convenience, means and interval estimates are produced from the posterior predictive distribution.

Value

Returns an object of class mbsts.prediction, which is a list.

Author(s)

Steven L. Scott

See Also

mbsts. predict.bsts plot.mbsts.prediction

Find the quarter in which a date occurs

Description

Returns the quarter and year in which a date occurs.

Usage

   Quarter(date)

Arguments

Value

A numeric vector identifying the quarter that each element of date corresponds to, expressed as a number of years since 1900. Thus Q1-2000 is 100.00, and Q3-2007 is 107.50.

Author(s)

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

Examples

  Quarter(c("2008-02-29", "2008-04-29"))
  # [1] 108.00 108.25

Regression Based Holiday Models

Description

Add a regression-based holiday model to the state specification.

Usage

AddRegressionHoliday(
    state.specification = NULL,
    y,
    holiday.list,
    time0 = NULL,
    prior = NULL,
    sdy = sd(as.numeric(y), na.rm = TRUE))

AddHierarchicalRegressionHoliday(
    state.specification = NULL,
    y,
    holiday.list,
    coefficient.mean.prior = NULL,
    coefficient.variance.prior = NULL,
    time0 = NULL,
    sdy = sd(as.numeric(y), na.rm = TRUE))

Arguments

Details

The model assumes that

y_t = \beta_{d(t)} + \epsilon_t %

The regression state model assumes vector of regression coefficients \beta contains elements \beta_d \sim N(0, \sigma).

The HierarchicalRegressionHolidayModel assumes \beta is composed of holiday-specific sub-vectors \beta_h \sim N(b_0, V), where each \beta_h contains coefficients describing the days in the influence window of holiday h. The hierarchical version of the model treats b_0 and V as parameters to be learned, with prior distributions

b_0 \sim N(\bar b, \Omega)

and

V \sim IW(\nu, S)

where IW represents the inverse Wishart distribution.

Value

Returns a list with the elements necessary to specify a local linear trend state model.

Author(s)

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

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. RandomWalkHolidayStateModel. SdPrior NormalPrior

Examples

trend <- cumsum(rnorm(730, 0, .1))
  dates <- seq.Date(from = as.Date("2014-01-01"), length = length(trend), by = "day")
  y <- zoo(trend + rnorm(length(trend), 0, .2), dates)

AddHolidayEffect <- function(y, dates, effect) {
  ## Adds a holiday effect to simulated data.
  ## Args:
  ##   y: A zoo time series, with Dates for indices.
  ##   dates: The dates of the holidays.
  ##   effect: A vector of holiday effects of odd length.  The central effect is
  ##     the main holiday, with a symmetric influence window on either side.
  ## Returns:
  ##   y, with the holiday effects added.
  time <- dates - (length(effect) - 1) / 2
  for (i in 1:length(effect)) {
    y[time] <- y[time] + effect[i]
    time <- time + 1
  }
  return(y)
}

## Define some holidays.
memorial.day <- NamedHoliday("MemorialDay")
memorial.day.effect <- c(.3, 3, .5)
memorial.day.dates <- as.Date(c("2014-05-26", "2015-05-25"))
y <- AddHolidayEffect(y, memorial.day.dates, memorial.day.effect)

presidents.day <- NamedHoliday("PresidentsDay")
presidents.day.effect <- c(.5, 2, .25)
presidents.day.dates <- as.Date(c("2014-02-17", "2015-02-16"))
y <- AddHolidayEffect(y, presidents.day.dates, presidents.day.effect)

labor.day <- NamedHoliday("LaborDay")
labor.day.effect <- c(1, 2, 1)
labor.day.dates <- as.Date(c("2014-09-01", "2015-09-07"))
y <- AddHolidayEffect(y, labor.day.dates, labor.day.effect)

## The holidays can be in any order.
holiday.list <- list(memorial.day, labor.day, presidents.day)

## In a real example you'd want more than 100 MCMC iterations.
niter <- 100

## Fit the model
ss <- AddLocalLevel(list(), y)
ss <- AddRegressionHoliday(ss, y, holiday.list = holiday.list)
model <- bsts(y, state.specification = ss, niter = niter)

## Plot all model state components.
plot(model, "comp")

## Plot the specific holiday state component.
plot(ss[[2]], model)

## Try again with some shrinkage.  With only 3 holidays there won't be much
## shrinkage.  
ss2 <- AddLocalLevel(list(), y)

## Plot the specific holiday state component.
ss2 <- AddHierarchicalRegressionHoliday(ss2, y, holiday.list = holiday.list)
model2 <- bsts(y, state.specification = ss2, niter = niter)

plot(model2, "comp")
plot(ss2[[2]], model2)

Produce a Regular Series of Time Stamps

Description

Given an set of timestamps that might contain duplicates and gaps, produce a set of timestamps that has no duplicates and no gaps.

Usage

  RegularizeTimestamps(timestamps)

  ## Default S3 method:
RegularizeTimestamps(timestamps)

  ## S3 method for class 'numeric'
RegularizeTimestamps(timestamps)

  ## S3 method for class 'Date'
RegularizeTimestamps(timestamps)

  ## S3 method for class 'POSIXt'
RegularizeTimestamps(timestamps)

Arguments

If the argument is NULL a NULL will be returned.

Value

A set of regularly spaced timestamps of the same class as the argument (which might be NULL).

Author(s)

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

Examples

  first <- as.POSIXct("2015-04-19 08:00:04")
  monthly <- seq(from = first, length.out = 24, by = "month")
  skip.one <- monthly[-8]
  has.duplicates <- monthly
  has.duplicates[2] <- has.duplicates[3]

  reg1 <- RegularizeTimestamps(skip.one)
  all.equal(reg1, monthly) ## TRUE

  reg2 <- RegularizeTimestamps(has.duplicates)
  all.equal(reg2, monthly)  ## TRUE

Residuals from a bsts Object

Description

Residuals (or posterior distribution of residuals) from a bsts object.

Usage

  ## S3 method for class 'bsts'
residuals(object,
    burn = SuggestBurn(.1, object),
    mean.only = FALSE,
    ...)

Arguments

Value

If mean.only is TRUE then this function returns a vector of residuals with the same "time stamp" as the original series. If mean.only is FALSE then the posterior distribution of the residuals is returned instead, as a matrix of draws. Each row of the matrix is an MCMC draw, and each column is a time point. The colnames of the returned matrix will be the timestamps of the original series, as text.

See Also

bsts, plot.bsts.

Retail sales, excluding food services

Description

A monthly time series of retail sales in the US, excluding food services. In millions of dollars. Seasonally adjusted.

Usage

data(rsxfs)

Format

zoo time series

Source

FRED. See http://research.stlouisfed.org/fred2/series/RSXFS

Examples

data(rsxfs)
plot(rsxfs)

Shark Attacks in Florida.

Description

An annual time series of shark attacks and fatalities in Florida.

Usage

data(shark)

Format

zoo time series

Source

From Jeffrey Simonoff "Analysis of Categorical Data". http://people.stern.nyu.edu/jsimonof/AnalCatData/Data/Comma_separated/floridashark.csv

Examples

data(shark)
head(shark)

Shorten long names

Description

Removes common prefixes and suffixes from character vectors.

Usage

  Shorten(words)

Arguments

Value

The argument words is returned, after common prefixes and suffixes have been removed. If all arguments are identical then no shortening is done.

Author(s)

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

See Also

bsts.mixed.

Examples

  Shorten(c("/usr/common/foo.tex", "/usr/common/barbarian.tex"))
  # returns c("foo", "barbarian")

  Shorten(c("hello", "hellobye"))
  # returns c("", "bye")

  Shorten(c("hello", "hello"))
  # returns c("hello", "hello")

  Shorten(c("", "x", "xx"))
  # returns c("", "x", "xx")

  Shorten("abcde")
  # returns "abcde"

Simulate fake mixed frequency data

Description

Simulate a fake data set that can be used to test mixed frequency code.

Usage

SimulateFakeMixedFrequencyData(nweeks,
                               xdim,
                               number.nonzero = xdim,
                               start.date = as.Date("2009-01-03"),
                               sigma.obs = 1.0,
                               sigma.slope = .5,
                               sigma.level = .5,
                               beta.sd = 10)

Arguments

Details

The simulation begins by simulating a local linear trend model for nweeks to get the trend component.

Next a nweeks by xdim matrix of predictor variables is simulated as IID normal(0, 1) deviates, and a xdim-vector of regression coefficients is simulated as IID normal(0, beta.sd). The product of the predictor matrix and regression coefficients is added to the output of the local linear trend model to get fine.target.

Finally, fine.target is aggregated to the month level to get coarse.target.

Value

Returns a list with the following components

Author(s)

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

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.mixed, AddLocalLinearTrend,

Examples

  fake.data <- SimulateFakeMixedFrequencyData(nweeks = 100, xdim = 10)
  plot(fake.data$coarse.target)

Spike and Slab Priors for AR Processes

Description

Returns a spike and slab prior for the parameters of an AR(p) process.

Usage

SpikeSlabArPrior(
    lags,
    prior.inclusion.probabilities =
        GeometricSequence( lags, initial.value = .8, discount.factor = .8),
    prior.mean = rep(0, lags),
    prior.sd =
        GeometricSequence(lags, initial.value = .5, discount.factor = .8),
    sdy,
    prior.df = 1,
    expected.r2 = .5,
    sigma.upper.limit = Inf,
    truncate = TRUE)

Arguments

Value

A list of class SpikeSlabArPrior containing the information needed for the underlying C++ code to instantiate this prior.

Author(s)

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

Compute state dimensions

Description

Returns a vector containing the size of each state component (i.e. the state dimension) in the state vector.

Usage

  StateSizes(state.specification)

Arguments

Value

A numeric vector giving the dimension of each state component.

Author(s)

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

Examples

    y <- rnorm(1000)
    state.specification <- AddLocalLinearTrend(list(), y)
    state.specification <- AddSeasonal(state.specification, y, 7)
    StateSizes(state.specification)

Summarize a Bayesian structural time series object

Description

Print a summary of a bsts object.

Usage

  ## S3 method for class 'bsts'
summary(object, burn = SuggestBurn(.1, object), ...)

Arguments

Value

Returns a list with the following elements.

References

Harvey's goodness of fit statistic is from Harvey (1989) Forecasting, structural time series models, and the Kalman filter. Page 268.

See Also

bsts, plot.bsts, summary.lm.spike

Examples

  data(AirPassengers)
  y <- log(AirPassengers)
  ss <- AddLocalLinearTrend(list(), y)
  ss <- AddSeasonal(ss, y, nseasons = 12)
  model <- bsts(y, state.specification = ss, niter = 100)
  summary(model, burn  = 20)

Convert to POSIXt

Description

Convert an object of class Date to class POSIXct without getting bogged down in timezone calculation.

Usage

  DateToPOSIX(timestamps)
  YearMonToPOSIX(timestamps)

Arguments

Details

Calling as.POSIXct on another date/time object (e.g. Date) applies a timezone correction to the object. This can shift the time marker by a few hours, which can have the effect of shifting the day by one unit. If the day was the first or last in a month or year, then the month or year will be off by one as well.

Coercing the object to the character representation of a Date prevents this adjustment from being applied, and leaves the POSIXt return value with the intended day, month, and year.

Author(s)

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

Turkish Electricity Usage

Description

A daily time series of electricity usaage in Turkey.

Usage

data(turkish)

Format

zoo time series

Source

https://robjhyndman.com/data/turkey_elec.csv

See Also

bsts

Examples

data(turkish)
plot(turkish)

Check to see if a week contains the end of a month or quarter

Description

Returns a logical vector indicating whether the given week contains the end of a month or quarter.

Usage

  WeekEndsMonth(week.ending)
  WeekEndsQuarter(week.ending)

Arguments

Value

A logical vector indicating whether the given week contains the end of a month or a quarter.

Author(s)

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

See Also

bsts.mixed.

Examples

  week.ending <- as.Date(c("2011-10-01",
                           "2011-10-08",
                           "2011-12-03",
                           "2011-12-31"))
  WeekEndsMonth(week.ending) == c(TRUE, FALSE, TRUE, TRUE)
  WeekEndsQuarter(week.ending) == c(TRUE, FALSE, FALSE, TRUE)

Days of the Week

Description

A character vector listing the names the days of the week.

Usage

weekday.names

See Also

month.name

Convert Between Wide and Long Format

Description

Convert a multivariate time series between wide and long formats. In "wide" format there is one row per time point, with series organzied by columns. In "long" format there is one row per observation, with variables indicating the series and time point to which an observation belongs.

Usage

  WideToLong(response, na.rm = TRUE)
  LongToWide(response, series.id, timestamps)

Arguments

Value

LongToWide returns a zoo matrix with the time series in wide format. WideToLong returns a 3-column data frame with columns "time", "series", and "values".

Author(s)

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

Examples

data(gdp)
gdp.wide <- LongToWide(gdp$GDP, gdp$Country, gdp$Time)
gdp.long <- WideToLong(gdp.wide)