bsm_ng {bssm}

Non-Gaussian Basic Structural (Time Series) Model

Description

Constructs a non-Gaussian basic structural model with local level or local trend component, a seasonal component, and regression component (or subset of these components).

Usage

bsm_ng(
  y,
  sd_level,
  sd_slope,
  sd_seasonal,
  sd_noise,
  distribution,
  phi,
  u,
  beta,
  xreg = NULL,
  period,
  a1 = NULL,
  P1 = NULL,
  C = NULL
)

Arguments

y

A vector or a ts object of observations.

sd_level

Standard deviation of the noise of level equation. Should be an object of class bssm_prior or scalar value defining a known value such as 0.

sd_slope

Standard deviation of the noise of slope equation. Should be an object of class bssm_prior, scalar value defining a known value such as 0, or missing, in which case the slope term is omitted from the model.

sd_seasonal

Standard deviation of the noise of seasonal equation. Should be an object of class bssm_prior, scalar value defining a known value such as 0, or missing, in which case the seasonal term is omitted from the model.

sd_noise

A prior for the standard deviation of the additional noise term to be added to linear predictor, defined as an object of class bssm_prior. If missing, no additional noise term is used.

distribution

Distribution of the observed time series. Possible choices are "poisson", "binomial", "gamma", and "negative binomial".

phi

Additional parameter relating to the non-Gaussian distribution. For negative binomial distribution this is the dispersion term, for gamma distribution this is the shape parameter, and for other distributions this is ignored. Should an object of class bssm_prior or a positive scalar.

u

A vector of positive constants for non-Gaussian models. For Poisson, gamma, and negative binomial distribution, this corresponds to the offset term. For binomial, this is the number of trials.

beta

A prior for the regression coefficients. Should be an object of class bssm_prior or bssm_prior_list (in case of multiple coefficients) or missing in case of no covariates.

xreg

A matrix containing covariates with number of rows matching the length of y. Can also be ts, mts or similar object convertible to matrix.

period

Length of the seasonal pattern. Must be a positive value greater than 2 and less than the length of the input time series. Default is frequency(y), which can also return non-integer value (in which case error is given).

a1

Prior means for the initial states (level, slope, seasonals). Defaults to vector of zeros.

P1

Prior covariance matrix for the initial states (level, slope, seasonals).Default is diagonal matrix with 100 on the diagonal.

C

Intercept terms for state equation, given as a m x n or m x 1 matrix.

Value

An object of class bsm_ng.

Examples

# Same data as in Vihola, Helske, Franks (2020)
data(poisson_series)
s <- sd(log(pmax(0.1, poisson_series)))
model <- bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
 sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2), 
 distribution = "poisson")


out <- run_mcmc(model, iter = 1e5, particles = 10)
summary(out, variable = "theta", return_se = TRUE)
# should be about 0.093 and 0.016
summary(out, variable = "states", return_se = TRUE, 
 states = 1, times = c(1, 100))
# should be about -0.075, 2.618


model <- bsm_ng(Seatbelts[, "VanKilled"], distribution = "poisson",
  sd_level = halfnormal(0.01, 1),
  sd_seasonal = halfnormal(0.01, 1),
  beta = normal(0, 0, 10),
  xreg = Seatbelts[, "law"],
  # default values, just for illustration
  period = 12L,
  a1 = rep(0, 1 + 11), # level + period - 1 seasonal states
  P1 = diag(1, 12),
  C = matrix(0, 12, 1),
  u = rep(1, nrow(Seatbelts)))


set.seed(123)
mcmc_out <- run_mcmc(model, iter = 5000, particles = 10, mcmc_type = "da")
mcmc_out$acceptance_rate
theta <- expand_sample(mcmc_out, "theta")
plot(theta)
summary(theta)

library("ggplot2")
ggplot(as.data.frame(theta[,1:2]), aes(x = sd_level, y = sd_seasonal)) +
  geom_point() + stat_density2d(aes(fill = ..level.., alpha = ..level..),
  geom = "polygon") + scale_fill_continuous(low = "green", high = "blue") +
  guides(alpha = "none")

# Traceplot using as.data.frame method for MCMC output
library("dplyr")
as.data.frame(mcmc_out) |> 
  filter(variable == "sd_level") |> 
  ggplot(aes(y = value, x = iter)) + geom_line()
  

# Model with slope term and additional noise to linear predictor to capture 
# excess variation   
model2 <- bsm_ng(Seatbelts[, "VanKilled"], distribution = "poisson",
  sd_level = halfnormal(0.01, 1),
  sd_seasonal = halfnormal(0.01, 1),
  beta = normal(0, 0, 10),
  xreg = Seatbelts[, "law"],
  sd_slope = halfnormal(0.01, 0.1),
  sd_noise = halfnormal(0.01, 1))

# instead of extra noise term, model using negative binomial distribution:
model3 <- bsm_ng(Seatbelts[, "VanKilled"], 
  distribution = "negative binomial",
  sd_level = halfnormal(0.01, 1),
  sd_seasonal = halfnormal(0.01, 1),
  beta = normal(0, 0, 10),
  xreg = Seatbelts[, "law"],
  sd_slope = halfnormal(0.01, 0.1),
  phi = gamma_prior(1, 5, 5))