bsts (version 0.9.5)

add.random.walk.holiday: 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

state.specification

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

y

The time series to be modeled, as a numeric vector convertible to xts. This state model assumes y contains daily data.

holiday

An object of class Holiday describing the influence window of the holiday being modeled.

time0

An object convertible to Date containing the date of the initial observation in the training data. If omitted and y is a zoo or xts object, then time0 will be obtained from the index of y[1].

sigma.prior

An object created by SdPrior describing the prior distribution for the standard deviation of the random walk increments.

initial.state.prior

An object created using NormalPrior, describing the prior distribution of the the initial state vector (at time 1).

sdy

The standard deviation of the series to be modeled. This will be ignored if y is provided, or if all the required prior distributions are supplied directly.

Value

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

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

# NOT RUN {
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)
# }