SSModel: Create a State Space Model Object of Class SSModel

In KFAS: Kalman Filter and Smoother for Exponential Family State Space Models

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
)

SSMbespoke(f, index, n)

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
)

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.
f	An user-defined function which should return the output from SSMcustom. See examples.
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.

SSMbespoke function is similar to SSMcustom but instead of defining the model component directly via system matrices, it wraps an arbitrary user-defined function which should return the output from SSMcustom. For more details, see package vignette.

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.

See Also

artransform

KFAS for more examples and link{SSMbespoke} for an alternative way of creating custom components.

Examples

# Example on bespoke function for time-varying trend
trend <- function(sigma, n) {
  Z <- array(seq_len(n), c(1, 1, n))
  T <- R <- matrix(1, 1, 1)
  Q <- matrix(sigma^2, 1, 1)
  a1 <- 0
  P1 <- 10
  SSMcustom(Z, T, R, Q, a1, P1, n = n, state_names = "timevarying trend")
}

model <- SSModel(Nile ~ SSMbespoke(trend(NA, length(Nile))), H = NA)
updatefn <- function(pars, model){
  model$Q[1, 1, 1] <- exp(0.5 * pars[1])
  model$H[1, 1, 1] <- exp(0.5 * pars[2])
  model
}

fit <- fitSSM(model, c(1, 20), updatefn, method = "BFGS")
conf_intv <- predict(fit$model, interval = "confidence", level = 0.95)

ts.plot(
  cbind(Nile, conf_intv), 
  col = c(1, 2, 2, 2),
  ylab = "Predicted Annual flow", 
  main = "River Nile"
) 
# 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
## Not run: 
# This fails as there is no x in the environment of f
try(c(SSModel(y ~ -1 + SSMregression(f), data = data1)$Z))

## End(Not run)

KFAS documentation built on June 8, 2025, 11:26 a.m.