SSMarima: Create a State Space Model Object of Class SSModel

Description

Function SSModel creates a state space object object of class SSModel which can be used as an input object for various functions of KFAS package.

Usage

SSMarima(
  ar = NULL,
  ma = NULL,
  d = 0,
  Q,
  stationary = TRUE,
  index,
  n = 1,
  state_names = NULL,
  ynames
)
SSMcustom(Z, T, R, Q, a1, P1, P1inf, index, n = 1, state_names = NULL)
SSMcycle(
  period,
  Q,
  type,
  index,
  a1,
  P1,
  P1inf,
  damping = 1,
  n = 1,
  state_names = NULL,
  ynames
)
SSModel(formula, data, H, u, distribution, tol = .Machine$double.eps^0.5)
SSMregression(
  rformula,
  data,
  type,
  Q,
  index,
  R,
  a1,
  P1,
  P1inf,
  n = 1,
  remove.intercept = TRUE,
  state_names = NULL,
  ynames
)
SSMseasonal(
  period,
  Q,
  sea.type = c("dummy", "trigonometric"),
  type,
  index,
  a1,
  P1,
  P1inf,
  n = 1,
  state_names = NULL,
  ynames,
  harmonics
)
SSMtrend(
  degree = 1,
  Q,
  type,
  index,
  a1,
  P1,
  P1inf,
  n = 1,
  state_names = NULL,
  ynames
)

Value

Object of class SSModel, which is a list with the following components:

y: A n x p matrix containing the observations.
Z: A p x m x 1 or p x m x n array corresponding to the system matrix of observation equation.
H: A p x p x 1 or p x p x n array corresponding to the covariance matrix of observational disturbances epsilon.
T: A m x m x 1 or m x m x n array corresponding to the first system matrix of state equation.
R: A m x k x 1 or m x k x n array corresponding to the second system matrix of state equation.
Q: A k x k x 1 or k x k x n array corresponding to the covariance matrix of state disturbances eta
a1: A m x 1 matrix containing the expected values of the initial states.
P1: A m x m matrix containing the covariance matrix of the nondiffuse part of the initial state vector.
P1inf: A m x m matrix containing the covariance matrix of the diffuse part of the initial state vector. If P1[i,i] is non-zero then P1inf[i,i] is automatically set to zero.
u: A n x p matrix of an additional parameters in case of non-Gaussian model.
distribution: A vector of length p giving the distributions of the observations.
tol: A tolerance parameter for Kalman filtering.
call: Original call to the function.

In addition, object of class

SSModel contains following attributes:

names: Names of the list components.
p, m, k, n: Integer valued scalars defining the dimensions of the model components.
state_types: Types of the states in the model.
eta_types: Types of the state disturbances in the model.
tv: Integer vector stating whether Z,H,T,R or Q is time-varying (indicated by 1 in tv and 0 otherwise). If you manually change the dimensions of the matrices you must change this attribute also.

Arguments

ar: For arima component, a numeric vector containing the autoregressive coeffients.
ma: For arima component, a numericvector containing the moving average coeffients.
d: For arima component, a degree of differencing.
Q: For arima, cycle and seasonal component, a $p \times p$ covariance matrix of the disturbances (or in the time varying case $p \times p \times n$ array), where where $p$ = length(index). For trend component, list of length degree containing the $p \times p$ or $p \times p \times n$ covariance matrices. For a custom component, arbitrary covariance matrix or array of disturbance terms $\eta_t$
stationary: For arima component, logical value indicating whether a stationarity of the arima part is assumed. Defaults to TRUE.
index: A vector indicating for which series the corresponding components are constructed.
n: Length of the series, only used internally for dimensionality check.
state_names: A character vector giving the state names.
ynames: names of the times series, used internally.
Z: For a custom component, system matrix or array of observation equation.
T: For a custom component, system matrix or array of transition equation.
R: For a custom and regression components, optional $m \times k$ system matrix or array of transition equation.
a1: Optional $m \times 1$ matrix giving the expected value of the initial state vector $\alpha_1$.
P1: Optional $m \times m$ matrix giving the covariance matrix of $\alpha_1$. In the diffuse case the non-diffuse part of $P_1$.
P1inf: Optional $m \times m$ matrix giving the diffuse part of $P_1$. Diagonal matrix with ones on diagonal elements which correspond to the diffuse initial states. If P1inf[i,i]>0, corresponding row and column of P1 should be zero.
period: For a cycle and seasonal components, the length of the cycle/seasonal pattern.
type: For cycle, seasonal, trend and regression components, character string defining if "distinct" or "common" states are used for different series.
damping: A damping factor for cycle component. Defaults to 1. Note that there are no checks for the range of the factor.
formula: An object of class formula containing the symbolic description of the model. The intercept term can be removed with -1 as in lm. In case of trend or differenced arima component the intercept is removed automatically in order to keep the model identifiable. See package vignette and examples in KFAS for special functions used in model construction.
data: An optional data frame, list or environment containing the variables in the model.
H: Covariance matrix or array of disturbance terms $\epsilon_t$ of observation equation. Defaults to an identity matrix. Omitted in case of non-Gaussian distributions (augment the state vector if you want to add additional noise).
u: Additional parameters for non-Gaussian models. See details in KFAS.
distribution: A vector of distributions of the observations. Default is rep("gaussian", p), where p is the number of series.
tol: A tolerance parameter used in checking whether Finf or F is numerically zero. Defaults to .Machine$double.eps^0.5. If F < tol * max(abs(Z[Z > 0]))^2, then F is deemed to be zero (i.e. differences are due to numerical precision). This has mostly effect only on determining when to end exact diffuse phase. Tweaking this and/or scaling model parameters/observations can sometimes help with numerical issues.
rformula: For regression component, right hand side formula or list of of such formulas defining the custom regression part.
remove.intercept: Remove intercept term from regression model. Default is TRUE. This tries to ensure that there are no extra intercept terms in the model.
sea.type: For seasonal component, character string defining whether to use "dummy" or "trigonometric" form of the seasonal component.
harmonics: For univariate trigonometric seasonal, argument harmonics can be used to specify which subharmonics are added to the model. Note that for multivariate model you can call SSMseasonal multiple times with different values of index.
degree: For trend component, integer defining the degree of the polynomial trend. 1 corresponds to local level, 2 for local linear trend and so forth.

Details

Formula of the model can contain the usual regression part and additional functions defining different types of components of the model, named as SSMarima, SSMcustom, SSMcycle, SSMregression, SSMseasonal and SSMtrend.

For more details, see package vignette (the mathematical notation is somewhat non-readable in ASCII).

Examples

Run this code

# add intercept to state equation by augmenting the state vector:
# diffuse initialization for the intercept, gets estimated like other states:
# for known fixed intercept, just set P1 = P1inf = 0 (default in SSMcustom).
intercept <- 0
model_int <- SSModel(Nile ~ SSMtrend(1, Q = 1469) +
SSMcustom(Z = 0, T = 1, Q = 0, a1 = intercept, P1inf = 1), H = 15099)

model_int$T
model_int$T[1, 2, 1] <- 1 # add the intercept value to level
out <- KFS(model_int)

# An example of a time-varying variance

model_drivers <- SSModel(log(cbind(front, rear)) ~ SSMtrend(1, Q = list(diag(2))),
data = Seatbelts, H = array(NA, c(2, 2, 192)))

ownupdatefn <- function(pars, model){
  diag(model$Q[, , 1]) <- exp(pars[1:2])
  model$H[,,1:169] <- diag(exp(pars[3:4])) # break in variance
  model$H[,,170:192] <- diag(exp(pars[5:6]))
  model
}

fit_drivers <- fitSSM(model_drivers, inits = rep(-1, 6),
  updatefn = ownupdatefn, method = "BFGS")
fit_drivers$model$H[,,1]
fit_drivers$model$H[,,192]

# An example of shift in the level component

Tt <- array(diag(2), c(2, 2, 100))
Tt[1,2,28] <- 1
Z <- matrix(c(1,0), 1, 2)
Q <- diag(c(NA, 0), 2)
model <- SSModel(Nile ~ -1 + SSMcustom(Z, Tt, Q = Q, P1inf = diag(2)),
  H = matrix(NA))

model <- fitSSM(model, c(10,10), method = "BFGS")$model
model$Q
model$H

conf_Nile <- predict(model, interval = "confidence", level = 0.9)
pred_Nile <- predict(model, interval = "prediction", level = 0.9)

ts.plot(cbind(Nile, pred_Nile, conf_Nile[, -1]), col = c(1:2, 3, 3, 4, 4),
        ylab = "Predicted Annual flow", main = "River Nile")

# dynamic regression model

set.seed(1)
x1 <- rnorm(100)
x2 <- rnorm(100)
b1 <- 1 + cumsum(rnorm(100, sd = 1))
b2 <- 2 + cumsum(rnorm(100, sd = 0.1))
y <- 1 + b1 * x1 + b2 * x2 + rnorm(100, sd = 0.1)

model <- SSModel(y ~ SSMregression(~ x1 + x2, Q = diag(NA,2)), H = NA)

fit <- fitSSM(model, inits = c(0,0,0), method = "BFGS")

model <- fit$model
model$Q
model$H
out <- KFS(model)

ts.plot(out$alphahat[,-1], b1, b2, col = 1:4)

# SSMregression with multivariate observations

x <- matrix(rnorm(30), 10, 3) # one variable per each series
y <- x + rnorm(30)
model <- SSModel(y ~ SSMregression(list(~ X1, ~ X2, ~ X3), data = data.frame(x)))
# more generally SSMregression(sapply(1:3, function(i) formula(paste0("~ X",i))), ...)

# three covariates per series, with same coefficients:
y <- x[,1] + x[,2] + x[,3] + matrix(rnorm(30), 10, 3)
model <- SSModel(y ~ -1 + SSMregression(~ X1 + X2 + X3, remove.intercept = FALSE,
  type = "common", data = data.frame(x)))

# the above cases can be combined in various ways, you can call SSMregression multiple times:
model <- SSModel(y ~  SSMregression(~ X1 + X2, type = "common") +
  SSMregression(~ X2), data = data.frame(x))

# examples of using data argument
y <- x <- rep(1, 3)
data1 <- data.frame(x = rep(2, 3))
data2 <- data.frame(x = rep(3, 3))

f <- formula(~ -1 + x)
# With data missing the environment of formula is checked,
# and if not found in there a calling environment via parent.frame is checked.

c(SSModel(y ~ -1 + x)["Z"]) # 1
c(SSModel(y ~ -1 + x, data = data1)["Z"]) # 2

c(SSModel(y ~ -1 + SSMregression(~ -1 + x))["Z"]) # 1
c(SSModel(y ~ -1 + SSMregression(~ -1 + x, data = data1))["Z"]) # 2
c(SSModel(y ~ -1 + SSMregression(~ -1 + x), data = data1)["Z"]) # 2
SSModel(y ~ -1 + x + SSMregression(~ -1 + x, data = data1))["Z"] # 1 and 2
SSModel(y ~ -1 + x + SSMregression(~ -1 + x), data = data1)["Z"] # both are 2
SSModel(y ~ -1 + x + SSMregression(~ -1 + x, data = data1), data = data2)["Z"] # 3 and 2

SSModel(y ~ -1 + x + SSMregression(f))["Z"] # 1 and 1
SSModel(y ~ -1 + x + SSMregression(f), data = data1)["Z"] # 2 and 1
SSModel(y ~ -1 + x + SSMregression(f,data = data1))["Z"] # 1 and 2

rm(x)
c(SSModel(y ~ -1 + SSMregression(f, data = data1))$Z) # 2
if (FALSE) {
# This fails as there is no x in the environment of f
try(c(SSModel(y ~ -1 + SSMregression(f), data = data1)$Z))
}