add.trig {bsts}R Documentation

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

state.specification

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

y

The time series to be modeled, as a numeric vector.

period

A positive scalar giving the number of time steps required for the longest cycle to repeat.

frequencies

A vector of positive real numbers giving the number of times each cyclic component repeats in a period. One sine and one cosine term will be added for each frequency.

sigma.prior

An object created by SdPrior describing the prior distribution for the standard deviation of the increments for the harmonic coefficients.

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.

method

The method of including the sinusoids. The "harmonic" method is strongly preferred, with "direct" offered mainly for teaching purposes.

Details

Harmonic Method

Each frequency lambdaj=2πj/Slambda_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: γjt\gamma_{jt}, and gammajtgamma^*_{jt}. They evolve through time as

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

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

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

The transition matrix multiplies the function by exp(iλj\exp(i \lambda_j, so that after 't' steps the harmonic's value is γexp(iλjt)\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 (γjt,γjt)(\gamma_{jt}, \gamma^*_{jt}), for j = 1, ... number of frequencies.

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

cos(λj)sin(λj)\cos(\lambda_j) \sin(\lambda_j)

sin(λj)cos(λ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.

yt=jβjtsin(fjt)+γjtcos(fjt)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].

βjt=βjt1+N(0,sigmasj2)γjt=γj1+N(0,sigmacj2)\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)


[Package bsts version 0.9.10 Index]