Learn R Programming

bssm (version 2.0.2)

ssm_ung: General univariate non-Gaussian state space model

Description

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

Usage

ssm_ung(
  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
)

Value

An object of class ssm_ung.

Arguments

y

Observations as time series (or vector) of length \(n\).

Z

System matrix Z of the observation equation. Either a vector of length m, a m x n matrix, or object which can be coerced to such.

T

System matrix T of the state equation. Either a m x m matrix or a m x m x n array, or object which can be coerced to such.

R

Lower triangular matrix R the state equation. Either a m x k matrix or a m x k x n array, or object which can be coerced to such.

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

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.

init_theta

Initial values for the unknown hyperparameters theta (i.e. unknown variables excluding latent state variables).

D

Intercept terms \(D_t\) for the observations equation, given as a scalar or vector of length n.

C

Intercept terms \(C_t\) for the state equation, given as a m times 1 or m times 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 univariate non-Gaussian model is defined using the following observational and state equations:

$$p(y_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(y_t | .)\) is either Poisson, binomial, gamma, or negative binomial distribution. Here k is the number of disturbance terms which can be less than m, the number of states.

The update_fn function should take only one vector argument which is used to create list with elements named as Z, phi T, R, a1, P1, D, and C, 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. Note that while you can input say R as m x k matrix for ssm_ung, update_fn should return R as m x k x 1 in this case. It might be useful to first construct the model without updating function and then check the expected structure of the model components from the output.

Examples

Run this code

data("drownings", package = "bssm")
model <- ssm_ung(drownings[, "deaths"], Z = 1, T = 1, R = 0.2, 
  a1 = 0, P1 = 10, distribution = "poisson", u = drownings[, "population"])

# approximate results based on Gaussian approximation
out <- smoother(model)
ts.plot(cbind(model$y / model$u, exp(out$alphahat)), col = 1:2)

Run the code above in your browser using DataLab