ssm_mng: General Non-Gaussian State Space Model

In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

General Non-Gaussian State Space Model

Description

Construct an object of class ssm_mng by directly defining the corresponding terms of the model.

Usage

ssm_mng(
  y,
  Z,
  T,
  R,
  a1 = NULL,
  P1 = NULL,
  distribution,
  phi = 1,
  u,
  init_theta = numeric(0),
  D = NULL,
  C = NULL,
  state_names,
  update_fn = default_update_fn,
  prior_fn = default_prior_fn
)

Arguments

y	Observations as multivariate time series or matrix with dimensions n x p.
Z	System matrix Z of the observation equation as p x m matrix or p x m x n array.
T	System matrix T of the state equation. Either a m x m matrix or a m x m x n array.
R	Lower triangular matrix R the state equation. Either a m x k matrix or a m x k x n array.
a1	Prior mean for the initial state as a vector of length m.
P1	Prior covariance matrix for the initial state as m x m matrix.
distribution	A vector of distributions of the observed series. Possible choices are "poisson", "binomial", "negative binomial", "gamma", and "gaussian".
phi	Additional parameters relating to the non-Gaussian distributions. For negative binomial distribution this is the dispersion term, for gamma distribution this is the shape parameter, for Gaussian this is standard deviation, and for other distributions this is ignored.
u	A matrix of positive constants for non-Gaussian models (of same dimensions as y). For Poisson, gamma, and negative binomial distribution, this corresponds to the offset term. For binomial, this is the number of trials (and as such should be integer(ish)).
init_theta	Initial values for the unknown hyperparameters theta (i.e. unknown variables excluding latent state variables).
D	Intercept terms for observation equation, given as p x n matrix.
C	Intercept terms for state equation, given as m x n matrix.
state_names	A character vector defining the names of the states.
update_fn	A function which returns list of updated model components given input vector theta. This function should take only one vector argument which is used to create list with elements named as Z, T, R, a1, P1, D, C, and phi, where each element matches the dimensions of the original model. If any of these components is missing, it is assumed to be constant wrt. theta. It's best to check the internal dimensions with str(model_object) as the dimensions of input arguments can differ from the final dimensions.
prior_fn	A function which returns log of prior density given input vector theta.

Details

The general multivariate non-Gaussian model is defined using the following observational and state equations:

p^i(y^i_t | D_t + Z_t \alpha_t), (\textrm{observation equation})

\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t, (\textrm{transition equation})

where \eta_t \sim N(0, I_k) and \alpha_1 \sim N(a_1, P_1) independently of each other, and p^i(y_t | .) is either Poisson, binomial, gamma, Gaussian, or negative binomial distribution for each observation series i=1,...,p. Here k is the number of disturbance terms (which can be less than m, the number of states).

Value

An object of class ssm_mng.

Examples

 
set.seed(1)
n <- 20
x <- cumsum(rnorm(n, sd = 0.5))
phi <- 2
y <- cbind(
  rgamma(n, shape = phi, scale = exp(x) / phi),
  rbinom(n, 10, plogis(x)))

Z <- matrix(1, 2, 1)
T <- 1
R <- 0.5
a1 <- 0
P1 <- 1

update_fn <- function(theta) {
  list(R = array(theta[1], c(1, 1, 1)), phi = c(theta[2], 1))
}

prior_fn <- function(theta) {
  ifelse(all(theta > 0), sum(dnorm(theta, 0, 1, log = TRUE)), -Inf)
}

model <- ssm_mng(y, Z, T, R, a1, P1, phi = c(2, 1), 
  init_theta = c(0.5, 2), 
  distribution = c("gamma", "binomial"),
  u = cbind(1, rep(10, n)),
  update_fn = update_fn, prior_fn = prior_fn,
  state_names = "random_walk",
  # using default values, but being explicit for testing purposes
  D = matrix(0, 2, 1), C = matrix(0, 1, 1))

# smoothing based on approximating gaussian model
ts.plot(cbind(y, fast_smoother(model)), 
  col = 1:3, lty = c(1, 1, 2))

bssm documentation built on Nov. 2, 2023, 6:25 p.m.