bssm (version 1.1.7-1)

bsm_ng: 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 = frequency(y),
  a1 = NULL,
  P1 = NULL,
  C = NULL
)

Arguments

y

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

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

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

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

Matrix containing covariates with number of rows matching the length of y.

period

Length of the seasonal pattern. Default is frequency(y). Must be a positive integer greater than 2 and less than the length of the input time series.

a1

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

P1

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

C

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

Value

Object of class bsm_ng.

Examples

# NOT RUN {
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 = 12,
  a1 = rep(0, 1 + 11), # level + period - 1 seasonal states
  P1 = diag(1, 12),
  C = matrix(0, 12, 1),
  u = rep(1, nrow(Seatbelts)))

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

# }