R/models.R
In helske/bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models
Defines functions
ssm_sde
ssm_nlg
ar1_lg
ar1_ng
svm
bsm_ng
bsm_lg
ssm_mng
ssm_mlg
ssm_ung
ssm_ulg
default_update_fn
default_prior_fn
Documented in
ar1_lg
ar1_ng
bsm_lg
bsm_ng
ssm_mlg
ssm_mng
ssm_nlg
ssm_sde
ssm_ulg
ssm_ung
svm
#' @srrstats {G2.3, G2.3a, G2.3b} match.arg and tolower used where applicable.
#' @srrstats {G2.7, G2.8, G2.9} Only matrix/mts/arrays as tabular data are
#' supported, not data.frame or similar objects.
#' @srrstats {G2.14, G2.14a, G2.14b, G2.14c, BS3.0} Missing observations are
#' handled automatically as per SSM theory, whereas missing values are not
#' allowed elsewhere.
#' @srrstats {BS1.0, BS1.1, BS1.2, BS1.2c} Examples and definitions of priors.
NULL
## placeholder functions for fixed models
default_prior_fn <- function(theta) {
0
}
default_update_fn <- function(theta) {
}
#'
#' General univariate linear-Gaussian state space models
#'
#' Construct an object of class \code{ssm_ulg} by directly defining the
#' corresponding terms of the model.
#'
#' The general univariate linear-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{y_t = D_t + Z_t \alpha_t + H_t \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, 1)}, \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other.
#' Here k is the number of disturbance terms which can be less than m, the
#' number of states.
#'
#' The \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{H} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_ulg},
#' \code{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.
#'
#' @inheritParams ssm_ung
#' @param H A vector of standard deviations. Either a scalar or a vector of
#' length n.
#' @param 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
#' \code{Z}, \code{H}, \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D}, and
#' \code{C}, where each element matches the dimensions of the original model
#' It's best to check the internal dimensions with \code{str(model_object)} as
#' the dimensions of input arguments can differ from the final dimensions.
#' If any of these components is missing, it is assumed to be constant wrt.
#' theta.
#' @return An object of class \code{ssm_ulg}.
#' @export
#' @examples
#'
#' # Regression model with time-varying coefficients
#' set.seed(1)
#' n <- 100
#' x1 <- rnorm(n)
#' x2 <- rnorm(n)
#' b1 <- 1 + cumsum(rnorm(n, sd = 0.5))
#' b2 <- 2 + cumsum(rnorm(n, sd = 0.1))
#' y <- 1 + b1 * x1 + b2 * x2 + rnorm(n, sd = 0.1)
#'
#' Z <- rbind(1, x1, x2)
#' H <- 0.1
#' T <- diag(3)
#' R <- diag(c(0, 1, 0.1))
#' a1 <- rep(0, 3)
#' P1 <- diag(10, 3)
#'
#' # updates the model given the current values of the parameters
#' update_fn <- function(theta) {
#' R <- diag(c(0, theta[1], theta[2]))
#' dim(R) <- c(3, 3, 1)
#' list(R = R, H = theta[3])
#' }
#' # prior for standard deviations as half-normal(1)
#' prior_fn <- function(theta) {
#' if(any(theta < 0)) {
#' log_p <- -Inf
#' } else {
#' log_p <- sum(dnorm(theta, 0, 1, log = TRUE))
#' }
#' log_p
#' }
#'
#' model <- ssm_ulg(y, Z, H, T, R, a1, P1,
#' init_theta = c(1, 0.1, 0.1),
#' update_fn = update_fn, prior_fn = prior_fn,
#' state_names = c("level", "b1", "b2"),
#' # using default values, but being explicit for testing purposes
#' C = matrix(0, 3, 1), D = numeric(1))
#'
#' out <- run_mcmc(model, iter = 5000)
#' out
#' sumr <- summary(out, variable = "state", times = 1:n)
#' sumr$true <- c(b1, b2, rep(1, n))
#' library(ggplot2)
#' ggplot(sumr, aes(x = time, y = Mean)) +
#' geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), alpha = 0.5) +
#' geom_line() +
#' geom_line(aes(y = true), colour = "red") +
#' facet_wrap(~ variable, scales = "free") +
#' theme_bw()
#'
#' # Perhaps easiest way to construct a general SSM for bssm is to use the
#' # model building functionality of KFAS:
#' library("KFAS")
#'
#' model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = 5e-4)+
#' SSMseasonal(period = 12, sea.type = "trigonometric", Q = 0) +
#' log(PetrolPrice) + law, data = Seatbelts, H = 0.005)
#'
#' # use as_bssm function for conversion, kappa defines the
#' # prior variance for diffuse states
#' model_bssm <- as_bssm(model_kfas, kappa = 100)
#'
#' # define updating function for parameter estimation
#' # we can use SSModel and as_bssm functions here as well
#' # (for large model it is more efficient to do this
#' # "manually" by constructing only necessary matrices,
#' # i.e., in this case a list with H and Q)
#'
#' prior_fn <- function(theta) {
#' if(any(theta < 0)) -Inf else sum(dnorm(theta, 0, 0.1, log = TRUE))
#' }
#'
#' update_fn <- function(theta) {
#'
#' model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
#' SSMseasonal(period = 12,
#' sea.type = "trigonometric", Q = theta[2]^2) +
#' log(PetrolPrice) + law, data = Seatbelts, H = theta[3]^2)
#'
#' # the bssm_model object is essentially list so this is fine
#' as_bssm(model_kfas, kappa = 100, init_theta = init_theta,
#' update_fn = update_fn, prior_fn = prior_fn)
#' }
#'
#' init_theta <- rep(1e-2, 3)
#' names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y")
#'
#' model_bssm <- update_fn(init_theta)
#'
#' \donttest{
#' out <- run_mcmc(model_bssm, iter = 10000, burnin = 5000)
#' out
#' }
#' # Above the regression coefficients are modelled as
#' # time-invariant latent states.
#' # Here is an alternative way where we use variable D so that the
#' # coefficients are part of parameter vector theta. Note however that the
#' # first option often preferable in order to keep the dimension of theta low.
#'
#' updatefn2 <- function(theta) {
#' # note no PetrolPrice or law variables here
#' model_kfas2 <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
#' SSMseasonal(period = 12, sea.type = "trigonometric", Q = theta[2]^2),
#' data = Seatbelts, H = theta[3]^2)
#'
#' X <- model.matrix(~ -1 + law + log(PetrolPrice), data = Seatbelts)
#' D <- t(X %*% theta[4:5])
#' as_bssm(model_kfas2, D = D, kappa = 100)
#' }
#' prior2 <- function(theta) {
#' if(any(theta[1:3] < 0)) {
#' -Inf
#' } else {
#' sum(dnorm(theta[1:3], 0, 0.1, log = TRUE)) +
#' sum(dnorm(theta[4:5], 0, 10, log = TRUE))
#' }
#' }
#' init_theta <- c(rep(1e-2, 3), 0, 0)
#' names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y", "law", "Petrol")
#' model_bssm2 <- updatefn2(init_theta)
#' model_bssm2$theta <- init_theta
#' model_bssm2$prior_fn <- prior2
#' model_bssm2$update_fn <- updatefn2
#' \donttest{
#' out2 <- run_mcmc(model_bssm2, iter = 10000, burnin = 5000)
#' out2
#' }
ssm_ulg <- function(y, Z, H, T, R, a1 = NULL, P1 = NULL,
init_theta = numeric(0),
D = NULL, C = NULL, state_names, update_fn = default_update_fn,
prior_fn = default_prior_fn) {
y <- check_y(y)
n <- length(y)
# create Z
Z <- check_Z(Z, 1L, n)
m <- dim(Z)[1]
# create T
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
H <- check_H(H, 1L, n)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
rownames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if (is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
# xreg and beta are need in C++ side in order to combine constructors
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R, a1 = a1, P1 = P1,
D = D, C = C, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
xreg = matrix(0, 0, 0), beta = numeric(0)),
class = c("ssm_ulg", "lineargaussian", "bssm_model"))
}
#' General univariate non-Gaussian state space model
#'
#' Construct an object of class \code{ssm_ung} by directly defining the
#' corresponding terms of the model.
#'
#' The general univariate non-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{p(y_t | D_t + Z_t \alpha_t), (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other,
#' and \eqn{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 \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{phi} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_ung},
#' \code{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.
#'
#' @inheritParams bsm_ng
#' @param y Observations as time series (or vector) of length \eqn{n}.
#' @param 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.
#' @param 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.
#' @param 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.
#' @param a1 Prior mean for the initial state as a vector of length m.
#' @param P1 Prior covariance matrix for the initial state as m x m matrix.
#' @param state_names A character vector defining the names of the states.
#' @param C Intercept terms \eqn{C_t} for the state equation, given as a
#' m times 1 or m times n matrix.
#' @param D Intercept terms \eqn{D_t} for the observations equation, given as a
#' scalar or vector of length n.
#' @param init_theta Initial values for the unknown hyperparameters theta
#' (i.e. unknown variables excluding latent state variables).
#' @param 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
#' \code{Z}, \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D}, \code{C}, and
#' \code{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
#' \code{str(model_object)} as the dimensions of input arguments can differ
#' from the final dimensions.
#' @param prior_fn A function which returns log of prior density
#' given input vector theta.
#' @return An object of class \code{ssm_ung}.
#' @export
#' @examples
#'
#' 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)
ssm_ung <- function(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) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, distribution = distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
Z <- check_Z(Z, 1L, n)
m <- dim(Z)[1]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
check_phi(phi)
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
rownames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
# xreg and beta are need in C++ side in order to combine constructors
structure(list(y = as.ts(y), Z = Z, T = T, R = R, a1 = a1, P1 = P1,
phi = phi, u = u, D = D, C = C, distribution = distribution,
initial_mode = initial_mode, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE,
xreg = matrix(0, 0, 0), beta = numeric(0)),
class = c("ssm_ung", "nongaussian", "bssm_model"))
}
#' General multivariate linear Gaussian state space models
#'
#' Construct an object of class \code{ssm_mlg} by directly defining the
#' corresponding terms of the model.
#'
#' The general multivariate linear-Gaussian model is defined using the
#' following observational and state equations:
#'
#' \deqn{y_t = D_t + Z_t \alpha_t + H_t \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, I_p)}, \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other.
#' Here p is the number of time series and k is the number of disturbance terms
#' (which can be less than m, the number of states).
#'
#' The \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{H} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_mlg},
#' \code{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
#'
#' @inheritParams ssm_ulg
#' @param y Observations as multivariate time series or matrix with
#' dimensions n x p.
#' @param Z System matrix Z of the observation equation as p x m matrix or
#' p x m x n array.
#' @param H Lower triangular matrix H of the observation. Either a scalar or
#' a vector of length n.
#' @param T System matrix T of the state equation. Either a m x m matrix or a
#' m x m x n array.
#' @param R Lower triangular matrix R the state equation. Either a m x k matrix
#' or a m x k x n array.
#' @param D Intercept terms for observation equation, given as a p x n matrix.
#' @param C Intercept terms for state equation, given as m x n matrix.
#' @return An object of class \code{ssm_mlg}.
#' @export
#' @examples
#'
#' data("GlobalTemp", package = "KFAS")
#' model_temp <- ssm_mlg(GlobalTemp, H = matrix(c(0.15,0.05,0, 0.05), 2, 2),
#' R = 0.05, Z = matrix(1, 2, 1), T = 1, P1 = 10,
#' state_names = "temperature",
#' # using default values, but being explicit for testing purposes
#' D = matrix(0, 2, 1), C = matrix(0, 1, 1))
#' ts.plot(cbind(model_temp$y, smoother(model_temp)$alphahat), col = 1:3)
#'
ssm_mlg <- function(y, Z, H, T, R, a1 = NULL, P1 = NULL,
init_theta = numeric(0), D = NULL, C = NULL, state_names,
update_fn = default_update_fn, prior_fn = default_prior_fn) {
# create y
y <- check_y(y, multivariate = TRUE)
n <- nrow(y)
p <- ncol(y)
# create Z
Z <- check_Z(Z, p, n, multivariate = TRUE)
m <- dim(Z)[2]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, p, n)
D <- as.matrix(D) # p = 1
C <- check_C(C, m, n)
H <- check_H(H, p, n, multivariate = TRUE)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
colnames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R, a1 = a1,
P1 = P1, D = D, C = C, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
state_names = state_names), class = c("ssm_mlg", "lineargaussian",
"bssm_model"))
}
#' General Non-Gaussian State Space Model
#'
#' Construct an object of class \code{ssm_mng} by directly defining the
#' corresponding terms of the model.
#'
#' The general multivariate non-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{p^i(y^i_t | D_t + Z_t \alpha_t), (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other, and
#' \eqn{p^i(y_t | .)} is either Poisson, binomial, gamma, Gaussian, or
#' negative binomial distribution for each observation series \eqn{i=1,...,p}.
#' Here k is the number of disturbance terms (which can be less than m,
#' the number of states).
#' @inheritParams ssm_ung
#' @param y Observations as multivariate time series or matrix with dimensions
#' n x p.
#' @param Z System matrix Z of the observation equation as p x m matrix or
#' p x m x n array.
#' @param T System matrix T of the state equation. Either a m x m matrix or a
#' m x m x n array.
#' @param R Lower triangular matrix R the state equation. Either a m x k
#' matrix or a
#' m x k x n array.
#' @param distribution A vector of distributions of the observed series.
#' Possible choices are
#' \code{"poisson"}, \code{"binomial"}, \code{"negative binomial"},
#' \code{"gamma"}, and \code{"gaussian"}.
#' @param 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.
#' @param 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)).
#' @param D Intercept terms for observation equation, given as p x n matrix.
#' @param C Intercept terms for state equation, given as m x n matrix.
#' @return An object of class \code{ssm_mng}.
#' @export
#' @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))
#'
ssm_mng <- function(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) {
# create y
y <- check_y(y, multivariate = TRUE)
n <- nrow(y)
p <- ncol(y)
if (missing(u)) u <- matrix(1, n, p)
u <- check_u(u, y, multivariate = TRUE)
if(length(distribution) == 1) distribution <- rep(distribution, p)
check_distribution(y, distribution)
if(length(phi) == 1) phi <- rep(phi, p)
for(i in 1:p) {
distribution[i] <- match.arg(tolower(distribution[i]),
c("poisson", "binomial", "negative binomial", "gamma", "gaussian"))
check_phi(phi[i])
}
# create Z
Z <- check_Z(Z, p, n, multivariate = TRUE)
m <- dim(Z)[2]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, p, n)
if (p == 1) D <- as.matrix(D)
C <- check_C(C, m, n)
initial_mode <- y
for(i in 1:p) {
initial_mode[, i] <- init_mode(y[, i], u[, i], distribution[i])
}
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
colnames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
structure(list(y = as.ts(y), Z = Z, T = T, R = R, a1 = a1, P1 = P1,
phi = phi, u = u, D = D, C = C, distribution = distribution,
initial_mode = initial_mode, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("ssm_mng", "nongaussian", "bssm_model"))
}
#' Basic Structural (Time Series) Model
#'
#' Constructs a basic structural model with local level or local trend
#' component and seasonal component.
#'
#' @inheritParams bsm_ng
#' @param sd_y Standard deviation of the noise of observation equation.
#' Should be an object of class \code{bssm_prior} or scalar
#' value defining a known value such as 0.
#' @param D Intercept terms for observation equation, given as a length n
#' numeric vector or a scalar in case of time-invariant intercept.
#' @param C Intercept terms for state equation, given as a m times n matrix
#' or m times 1 matrix in case of time-invariant intercept.
#' @return An object of class \code{bsm_lg}.
#' @export
#' @examples
#'
#' set.seed(1)
#' n <- 50
#' x <- rnorm(n)
#' level <- numeric(n)
#' level[1] <- rnorm(1)
#' for (i in 2:n) level[i] <- rnorm(1, -0.2 + level[i-1], sd = 0.1)
#' y <- rnorm(n, 2.1 + x + level)
#' model <- bsm_lg(y, sd_y = halfnormal(1, 5), sd_level = 0.1, a1 = level[1],
#' P1 = matrix(0, 1, 1), xreg = x, beta = normal(1, 0, 1),
#' D = 2.1, C = matrix(-0.2, 1, 1))
#'
#' ts.plot(cbind(fast_smoother(model), level), col = 1:2)
#'
#' prior <- uniform(0.1 * sd(log10(UKgas)), 0, 1)
#' # period here is redundant as frequency(UKgas) = 4
#' model_UKgas <- bsm_lg(log10(UKgas), sd_y = prior, sd_level = prior,
#' sd_slope = prior, sd_seasonal = prior, period = 4)
#'
#' # Note small number of iterations for CRAN checks
#' mcmc_out <- run_mcmc(model_UKgas, iter = 5000)
#' summary(mcmc_out, return_se = TRUE)
#' # Use the summary method from coda:
#' summary(expand_sample(mcmc_out, "theta"))$stat
#' mcmc_out$theta[which.max(mcmc_out$posterior), ]
#' sqrt((fit <- StructTS(log10(UKgas), type = "BSM"))$coef)[c(4, 1:3)]
#'
bsm_lg <- function(y, sd_y, sd_level, sd_slope, sd_seasonal,
beta, xreg = NULL, period, a1 = NULL, P1 = NULL, D = NULL,
C = NULL) {
y <- check_y(y)
n <- length(y)
regression_part <- create_regression(beta, xreg, n)
notfixed <- c("y" = 1, "level" = 1, "slope" = 1, "seasonal" = 1)
if (missing(sd_y) || is.null(sd_y)) {
stop("Provide either prior or fixed value for sd_y.")
} else {
if (is_prior(sd_y)) {
check_sd(sd_y$init, "y")
H <- matrix(sd_y$init)
} else {
notfixed[1] <- 0
check_sd(sd_y, "y")
H <- matrix(sd_y)
}
}
if (missing(sd_level) || is.null(sd_level)) {
stop("Provide either prior or fixed value for sd_level.")
} else {
if (is_prior(sd_level)) {
check_sd(sd_level$init, "level")
} else {
notfixed["level"] <- 0
check_sd(sd_level, "level")
}
}
if (missing(sd_slope) || is.null(sd_slope)) {
notfixed["slope"] <- 0
slope <- FALSE
sd_slope <- NULL
} else {
if (is_prior(sd_slope)) {
check_sd(sd_slope$init, "slope")
} else {
notfixed["slope"] <- 0
check_sd(sd_slope, "slope")
}
slope <- TRUE
}
if (missing(sd_seasonal) || is.null(sd_seasonal)) {
notfixed["seasonal"] <- 0
seasonal_names <- NULL
seasonal <- FALSE
sd_seasonal <- NULL
period <- 1L
} else {
if (missing(period)) period <- frequency(y)
period <- check_period(period, n)
if (is_prior(sd_seasonal)) {
check_sd(sd_seasonal$init, "seasonal")
} else {
notfixed["seasonal"] <- 0
check_sd(sd_seasonal, "seasonal")
}
seasonal <- TRUE
seasonal_names <- paste0("seasonal_", 1:(period - 1))
}
npar_R <- 1L + as.integer(slope) + as.integer(seasonal)
m <- as.integer(1L + as.integer(slope) + as.integer(seasonal) * (period - 1))
a1 <- check_a1(a1, m)
if (missing(P1)) {
P1 <- diag(100, m)
} else {
P1 <- check_P1(P1, m)
}
if (slope) {
state_names <- c("level", "slope", seasonal_names)
} else {
state_names <- c("level", seasonal_names)
}
Z <- matrix(0, m, 1)
Z[1, 1] <- 1
if (seasonal) {
Z[2 + slope, 1] <- 1
}
T <- matrix(0, m, m)
T[1, 1] <- 1
if (slope) {
T[1:2, 2] <- 1
}
if (seasonal) {
T[(2 + slope), (2 + slope):m] <- -1
T[cbind(1 + slope + 2:(period - 1), 1 + slope + 1:(period - 2))] <- 1
}
R <- matrix(0, m, max(1, npar_R))
if (notfixed["level"]) {
R[1, 1] <- sd_level$init
} else {
R[1, 1] <- sd_level
}
if (slope) {
if (notfixed["slope"]) {
R[2, 2] <- sd_slope$init
} else {
R[2, 2] <- sd_slope
}
}
if (seasonal) {
if (notfixed["seasonal"]) {
R[2 + slope, 2 + slope] <- sd_seasonal$init
} else {
R[2 + slope, 2 + slope] <- sd_seasonal
}
}
dim(H) <- 1
dim(T) <- c(m, m, 1)
dim(R) <- c(m, ncol(R), 1)
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- state_names
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(sd_y, sd_level, sd_slope, sd_seasonal),
regression_part$beta)
} else {
priors <- list(sd_y, sd_level, sd_slope, sd_seasonal,
regression_part$beta)
}
names(priors) <- c("sd_y", "sd_level", "sd_slope", "sd_seasonal",
names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R,
a1 = a1, P1 = P1, xreg = regression_part$xreg,
beta = regression_part$coefs, D = D, C = C,
slope = slope, seasonal = seasonal, period = period,
fixed = as.integer(!notfixed),
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta), class = c("bsm_lg", "ssm_ulg", "lineargaussian",
"bssm_model"))
}
#' Non-Gaussian Basic Structural (Time Series) Model
#'
#' 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).
#'
#' @param y A vector or a \code{ts} object of observations.
#' @param sd_level Standard deviation of the noise of level equation.
#' Should be an object of class \code{bssm_prior} or scalar
#' value defining a known value such as 0.
#' @param sd_slope Standard deviation of the noise of slope equation.
#' Should be an object of class \code{bssm_prior}, scalar
#' value defining a known value such as 0, or missing, in which case the slope
#' term is omitted from the model.
#' @param sd_seasonal Standard deviation of the noise of seasonal equation.
#' Should be an object of class \code{bssm_prior}, scalar
#' value defining a known value such as 0, or missing, in which case the
#' seasonal term is omitted from the model.
#' @param sd_noise A prior for the standard deviation of the additional noise
#' term to be added to linear predictor, defined as an object of class
#' \code{bssm_prior}. If missing, no additional noise term is used.
#' @param distribution Distribution of the observed time series. Possible
#' choices are \code{"poisson"}, \code{"binomial"}, \code{"gamma"}, and
#' \code{"negative binomial"}.
#' @param 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 \code{bssm_prior} or
#' a positive scalar.
#' @param 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.
#' @param beta A prior for the regression coefficients.
#' Should be an object of class \code{bssm_prior} or \code{bssm_prior_list}
#' (in case of multiple coefficients) or missing in case of no covariates.
#' @param xreg A matrix containing covariates with number of rows matching the
#' length of \code{y}. Can also be \code{ts}, \code{mts} or similar object
#' convertible to matrix.
#' @param period Length of the seasonal pattern.
#' Must be a positive value greater than 2 and less than the length of the
#' input time series. Default is \code{frequency(y)},
#' which can also return non-integer value (in which case error is given).
#' @param a1 Prior means for the initial states (level, slope, seasonals).
#' Defaults to vector of zeros.
#' @param P1 Prior covariance matrix for the initial states (level, slope,
#' seasonals).Default is diagonal matrix with 100 on the diagonal.
#' @param C Intercept terms for state equation, given as a m x n or m x 1
#' matrix.
#' @return An object of class \code{bsm_ng}.
#' @export
#' @examples
#' # Same data as in Vihola, Helske, Franks (2020)
#' data(poisson_series)
#' s <- sd(log(pmax(0.1, poisson_series)))
#' model <- bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
#' sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2),
#' distribution = "poisson")
#'
#' \donttest{
#' out <- run_mcmc(model, iter = 1e5, particles = 10)
#' summary(out, variable = "theta", return_se = TRUE)
#' # should be about 0.093 and 0.016
#' summary(out, variable = "states", return_se = TRUE,
#' states = 1, times = c(1, 100))
#' # should be about -0.075, 2.618
#' }
#'
#' 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 = 12L,
#' a1 = rep(0, 1 + 11), # level + period - 1 seasonal states
#' P1 = diag(1, 12),
#' C = matrix(0, 12, 1),
#' u = rep(1, nrow(Seatbelts)))
#'
#' \donttest{
#' 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()
#'
#' }
#' # 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))
#'
bsm_ng <- function(y, sd_level, sd_slope, sd_seasonal, sd_noise,
distribution, phi, u, beta, xreg = NULL, period,
a1 = NULL, P1 = NULL, C = NULL) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, multivariate = FALSE, distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
regression_part <- create_regression(beta, xreg, n)
notfixed <- c("level" = 1, "slope" = 1, "seasonal" = 1)
if (missing(sd_level) || missing(sd_level)) {
stop("Provide either prior or fixed value for sd_level.")
} else {
if (is_prior(sd_level)) {
check_sd(sd_level$init, "level")
} else {
notfixed["level"] <- 0
check_sd(sd_level, "level")
}
}
if (missing(sd_slope) || is.null(sd_slope)) {
notfixed["slope"] <- 0
slope <- FALSE
sd_slope <- NULL
} else {
if (is_prior(sd_slope)) {
check_sd(sd_slope$init, "slope")
} else {
notfixed["slope"] <- 0
check_sd(sd_slope, "slope")
}
slope <- TRUE
}
if (missing(sd_seasonal) || is.null(sd_seasonal)) {
notfixed["seasonal"] <- 0
seasonal_names <- NULL
seasonal <- FALSE
sd_seasonal <- NULL
period <- 1L
} else {
if (missing(period)) period <- frequency(y)
period <- check_period(period, n)
if (is_prior(sd_seasonal)) {
check_sd(sd_seasonal$init, "seasonal")
} else {
notfixed["seasonal"] <- 0
check_sd(sd_seasonal, "seasonal")
}
seasonal <- TRUE
seasonal_names <- paste0("seasonal_", 1:(period - 1))
}
if (missing(sd_noise) || is.null(sd_noise)) {
noise <- FALSE
sd_noise <- NULL
} else {
check_sd(sd_noise$init, "noise")
noise <- TRUE
}
npar_R <- 1L + as.integer(slope) + as.integer(seasonal) + as.integer(noise)
m <- as.integer(1L + as.integer(slope) +
as.integer(seasonal) * (period - 1) + as.integer(noise))
a1 <- check_a1(a1, m)
if (missing(P1)) {
P1 <- diag(100, m)
} else {
P1 <- check_P1(P1, m)
}
if (slope) {
state_names <- c("level", "slope", seasonal_names)
} else {
state_names <- c("level", seasonal_names)
}
Z <- matrix(0, m, 1)
Z[1, 1] <- 1
if (seasonal) {
Z[2 + slope, 1] <- 1
}
T <- matrix(0, m, m)
T[1, 1] <- 1
if (slope) {
T[1:2, 2] <- 1
}
if (seasonal) {
T[(2 + slope), (2 + slope):m] <- -1
T[cbind(1 + slope + 2:(period - 1), 1 + slope + 1:(period - 2))] <- 1
}
R <- matrix(0, m, max(1, npar_R))
if (notfixed["level"]) {
R[1, 1] <- sd_level$init
} else {
R[1, 1] <- sd_level
}
if (slope) {
if (notfixed["slope"]) {
R[2, 2] <- sd_slope$init
} else {
R[2, 2] <- sd_slope
}
}
if (seasonal) {
if (notfixed["seasonal"]) {
R[2 + slope, 2 + slope] <- sd_seasonal$init
} else {
R[2 + slope, 2 + slope] <- sd_seasonal
}
}
#additional noise term
if (noise) {
P1[m, m] <- sd_noise$init^2
Z[m] <- 1
state_names <- c(state_names, "noise")
R[m, max(1, ncol(R) - 1)] <- sd_noise$init
}
phi_est <- FALSE
use_phi <- distribution %in% c("negative binomial", "gamma")
if (use_phi) {
if (is_prior(phi)) {
check_phi(phi$init)
phi_est <- TRUE
} else {
check_phi(phi)
}
} else {
phi <- 1
}
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
dim(T) <- c(m, m, 1)
dim(R) <- c(m, ncol(R), 1)
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- state_names
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(sd_level, sd_slope, sd_seasonal, sd_noise, phi),
regression_part$beta)
} else {
priors <- list(sd_level, sd_slope, sd_seasonal, sd_noise, phi,
regression_part$beta)
}
names(priors) <- c("sd_level", "sd_slope", "sd_seasonal", "sd_noise", "phi",
names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
if (phi_est) {
phi <- phi$init
}
D <- numeric(1)
C <- check_C(C, m, n)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = phi, u = u, xreg = regression_part$xreg,
beta = regression_part$coefs, D = D, C = C,
slope = slope, seasonal = seasonal, noise = noise,
period = period, fixed = as.integer(!notfixed),
distribution = distribution, initial_mode = initial_mode,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta, phi_est = phi_est,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("bsm_ng", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Stochastic Volatility Model
#'
#' Constructs a simple stochastic volatility model with Gaussian errors and
#' first order autoregressive signal. See the main vignette for details.
#'
#' @param y A numeric vector or a \code{\link{ts}} object of observations.
#' @param mu A prior for mu parameter of transition equation.
#' Should be an object of class \code{bssm_prior}.
#' @param rho A prior for autoregressive coefficient.
#' Should be an object of class \code{bssm_prior}.
#' @param sd_ar A prior for the standard deviation of noise of the AR-process.
#' Should be an object of class \code{bssm_prior}.
#' @param sigma A prior for sigma parameter of observation equation, internally
#' denoted as phi. Should be an object of class \code{bssm_prior}.
#' Ignored if \code{mu} is provided. Note that typically
#' parametrization using mu is preferred due to better numerical properties and
#' availability of better Gaussian approximation.
#' Most notably the global approximation approach does not work with sigma
#' parameterization as sigma is not a parameter of the resulting approximate
#' model.
#' @return An object of class \code{svm}.
#' @export
#' @rdname svm
#' @examples
#'
#' data("exchange")
#' y <- exchange[1:100] # for faster CRAN check
#' model <- svm(y, rho = uniform(0.98, -0.999, 0.999),
#' sd_ar = halfnormal(0.15, 5), sigma = halfnormal(0.6, 2))
#'
#' obj <- function(pars) {
#' -logLik(svm(y,
#' rho = uniform(pars[1], -0.999, 0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' sigma = halfnormal(pars[3], 2)), particles = 0)
#' }
#' opt <- optim(c(0.98, 0.15, 0.6), obj,
#' lower = c(-0.999, 1e-4, 1e-4),
#' upper = c(0.999, 10, 10), method = "L-BFGS-B")
#' pars <- opt$par
#' model <- svm(y,
#' rho = uniform(pars[1],-0.999,0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' sigma = halfnormal(pars[3], 2))
#'
#' # alternative parameterization
#' model2 <- svm(y, rho = uniform(0.98,-0.999, 0.999),
#' sd_ar = halfnormal(0.15, 5), mu = normal(0, 0, 1))
#'
#' obj2 <- function(pars) {
#' -logLik(svm(y,
#' rho = uniform(pars[1], -0.999, 0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' mu = normal(pars[3], 0, 1)), particles = 0)
#' }
#' opt2 <- optim(c(0.98, 0.15, 0), obj2, lower = c(-0.999, 1e-4, -Inf),
#' upper = c(0.999, 10, Inf), method = "L-BFGS-B")
#' pars2 <- opt2$par
#' model2 <- svm(y,
#' rho = uniform(pars2[1],-0.999,0.999),
#' sd_ar = halfnormal(pars2[2], 5),
#' mu = normal(pars2[3], 0, 1))
#'
#' # sigma is internally stored in phi
#' ts.plot(cbind(model$phi * exp(0.5 * fast_smoother(model)),
#' exp(0.5 * fast_smoother(model2))), col = 1:2)
#'
svm <- function(y, mu, rho, sd_ar, sigma) {
if(!missing(sigma) && !missing(mu)) {
stop("Define either sigma or mu, but not both.")
}
y <- check_y(y)
check_rho(rho$init)
check_sd(sd_ar$init, "rho")
if(!missing(mu)) {
svm_type <- 1L
check_mu(mu$init)
initial_mode <- matrix(log(pmax(1e-4, y^2)), ncol = 1)
} else {
svm_type <- 0L
check_sd(sigma$init, "sigma", FALSE)
initial_mode <-
matrix(log(pmax(1e-4, y^2)) - 2 * log(sigma$init), ncol = 1)
}
a1 <- if(svm_type) mu$init else 0
P1 <- matrix(sd_ar$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sd_ar$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
priors <- list(rho, sd_ar, if(svm_type==0) sigma else mu)
priors <- priors[!vapply(priors, is.null, TRUE)]
names(priors) <-
c("rho", "sd_ar", if(svm_type==0) "sigma" else "mu")
C <- if (svm_type) matrix(mu$init * (1 - T[1])) else matrix(0)
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = if (svm_type == 0) sigma$init else 1,
xreg = matrix(0, 0, 0),
beta = numeric(0), D = D, C = C,
initial_mode = initial_mode,
svm_type = svm_type, distribution = "svm", u = 1,
phi_est = !as.logical(svm_type),
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("svm", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Non-Gaussian model with AR(1) latent process
#'
#' Constructs a simple non-Gaussian model where the state dynamics follow an
#' AR(1) process.
#'
#' @inheritParams bsm_ng
#' @param rho A prior for autoregressive coefficient.
#' Should be an object of class \code{bssm_prior}.
#' @param mu A fixed value or a prior for the stationary mean of the latent
#' AR(1) process. Should be an object of class \code{bssm_prior} or scalar
#' value defining a fixed mean such as 0.
#' @param sigma A prior for the standard deviation of noise of the AR-process.
#' Should be an object of class \code{bssm_prior}
#' @return An object of class \code{ar1_ng}.
#' @export
#' @rdname ar1_ng
#' @examples
#' model <- ar1_ng(discoveries, rho = uniform(0.5,-1,1),
#' sigma = halfnormal(0.1, 1), mu = normal(0, 0, 1),
#' distribution = "poisson")
#' out <- run_mcmc(model, iter = 1e4, mcmc_type = "approx",
#' output_type = "summary")
#'
#' ts.plot(cbind(discoveries, exp(out$alphahat)), col = 1:2)
#'
#' set.seed(1)
#' n <- 30
#' phi <- 2
#' rho <- 0.9
#' sigma <- 0.1
#' beta <- 0.5
#' u <- rexp(n, 0.1)
#' x <- rnorm(n)
#' z <- y <- numeric(n)
#' z[1] <- rnorm(1, 0, sigma / sqrt(1 - rho^2))
#' y[1] <- rnbinom(1, mu = u * exp(beta * x[1] + z[1]), size = phi)
#' for(i in 2:n) {
#' z[i] <- rnorm(1, rho * z[i - 1], sigma)
#' y[i] <- rnbinom(1, mu = u * exp(beta * x[i] + z[i]), size = phi)
#' }
#'
#' model <- ar1_ng(y, rho = uniform_prior(0.9, 0, 1),
#' sigma = gamma_prior(0.1, 2, 10), mu = 0.,
#' phi = gamma_prior(2, 2, 1), distribution = "negative binomial",
#' xreg = x, beta = normal_prior(0.5, 0, 1), u = u)
#'
ar1_ng <- function(y, rho, sigma, mu, distribution, phi, u, beta,
xreg = NULL) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, multivariate = FALSE, distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
regression_part <- create_regression(beta, xreg, n)
check_rho(rho$init)
check_sd(sigma$init, "rho")
if (is_prior(mu)) {
check_mu(mu$init)
mu_est <- TRUE
a1 <- mu$init
C <- matrix(mu$init * (1 - rho$init))
} else {
mu_est <- FALSE
check_mu(mu)
a1 <- mu
C <- matrix(mu * (1 - rho$init))
}
distribution <- match.arg(tolower(distribution), c("poisson", "binomial",
"negative binomial", "gamma"))
use_phi <- distribution %in% c("negative binomial", "gamma")
phi_est <- FALSE
if (use_phi) {
if (is_prior(phi)) {
check_phi(phi$init)
phi_est <- TRUE
} else {
check_phi(phi)
}
} else {
phi <- 1
}
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
P1 <- matrix(sigma$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sigma$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(rho, sigma, mu, phi), regression_part$beta)
} else {
priors <- list(rho, sigma, mu, phi, regression_part$beta)
}
names(priors) <-
c("rho", "sigma", "mu", "phi", names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
if (phi_est) {
phi <- phi$init
}
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = phi, u = u,
xreg = regression_part$xreg, beta = regression_part$coefs,
D = D, C = C,
initial_mode = initial_mode,
distribution = distribution, mu_est = mu_est, phi_est = phi_est,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters, theta = theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("ar1_ng", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Univariate Gaussian model with AR(1) latent process
#'
#' Constructs a simple Gaussian model where the state dynamics
#' follow an AR(1) process.
#'
#' @inheritParams ar1_ng
#' @param sd_y A prior for the standard deviation of observation equation.
#' @return An object of class \code{ar1_lg}.
#' @export
#' @rdname ar1_lg
#' @examples
#' set.seed(1)
#' mu <- 2
#' rho <- 0.7
#' sd_y <- 0.1
#' sigma <- 0.5
#' beta <- -1
#' x <- rnorm(30)
#' z <- y <- numeric(30)
#' z[1] <- rnorm(1, mu, sigma / sqrt(1 - rho^2))
#' y[1] <- rnorm(1, beta * x[1] + z[1], sd_y)
#' for(i in 2:30) {
#' z[i] <- rnorm(1, mu * (1 - rho) + rho * z[i - 1], sigma)
#' y[i] <- rnorm(1, beta * x[i] + z[i], sd_y)
#' }
#' model <- ar1_lg(y, rho = uniform(0.5, -1, 1),
#' sigma = halfnormal(1, 10), mu = normal(0, 0, 1),
#' sd_y = halfnormal(1, 10),
#' xreg = x, beta = normal(0, 0, 1))
#' out <- run_mcmc(model, iter = 2e4)
#' summary(out, return_se = TRUE)
#'
ar1_lg <- function(y, rho, sigma, mu, sd_y, beta, xreg = NULL) {
y <- check_y(y)
n <- length(y)
regression_part <- create_regression(beta, xreg, n)
check_rho(rho$init)
check_sd(sigma$init, "rho")
if (is_prior(mu)) {
check_mu(mu$init)
mu_est <- TRUE
a1 <- mu$init
C <- matrix(mu$init * (1 - rho$init))
} else {
mu_est <- FALSE
check_mu(mu)
a1 <- mu
C <- matrix(mu * (1 - rho$init))
}
if (is_prior(sd_y)) {
check_sd(sd_y$init, "y")
sd_y_est <- TRUE
H <- matrix(sd_y$init)
} else {
sd_y_est <- FALSE
check_sd(sd_y, "y")
H <- matrix(sd_y)
}
P1 <- matrix(sigma$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sigma$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(rho, sigma, mu, sd_y), regression_part$beta)
} else {
priors <- list(rho, sigma, mu, sd_y, regression_part$beta)
}
names(priors) <-
c("rho", "sigma", "mu", "sd_y", names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R,
a1 = a1, P1 = P1,
xreg = regression_part$xreg, beta = regression_part$coefs,
D = D, C = C,
mu_est = mu_est, sd_y_est = sd_y_est,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters, theta = theta,
max_iter = 100, conv_tol = 1e-8),
class = c("ar1_lg", "ssm_ulg", "lineargaussian", "bssm_model"))
}
#'
#' General multivariate nonlinear Gaussian state space models
#'
#' Constructs an object of class \code{ssm_nlg} by defining the corresponding
#' terms of the observation and state equation.
#'
#' The nonlinear Gaussian model is defined as
#'
#' \deqn{y_t = Z(t, \alpha_t, \theta) + H(t, \theta) \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = T(t, \alpha_t, \theta) + R(t, \theta)\eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, I_p)}, \eqn{\eta_t \sim N(0, I_m)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other, and functions
#' \eqn{Z, H, T, R} can depend on \eqn{\alpha_t} and parameter vector
#' \eqn{\theta}.
#'
#' Compared to other models, these general models need a bit more effort from
#' the user, as you must provide the several small C++ snippets which define the
#' model structure. See examples in the vignette and \code{cpp_example_model}.
#'
#' @param y Observations as multivariate time series (or matrix) of length
#' \eqn{n}.
#' @param Z,H,T,R An external pointers (object of class \code{externalptr})
#' for the C++ functions which define the corresponding model functions.
#' @param Z_gn,T_gn An external pointers (object of class \code{externalptr})
#' for the C++ functions which define the gradients of the corresponding model
#' functions.
#' @param a1 Prior mean for the initial state as object of class
#' \code{externalptr}
#' @param P1 Prior covariance matrix for the initial state as object of class
#' \code{externalptr}
#' @param theta Parameter vector passed to all model functions.
#' @param known_params A vector of known parameters passed to all model
#' functions.
#' @param known_tv_params A matrix of known parameters passed to all model
#' functions.
#' @param n_states Number of states in the model (positive integer).
#' @param n_etas Dimension of the noise term of the transition equation
#' (positive integer).
#' @param log_prior_pdf An external pointer (object of class
#' \code{externalptr}) for the C++ function which
#' computes the log-prior density given theta.
#' @param time_varying Optional logical vector of length 4, denoting whether
#' the values of
#' Z, H, T, and R vary with respect to time variable (given identical states).
#' If used, this can speed up some computations.
#' @param state_names A character vector containing names for the states.
#' @return An object of class \code{ssm_nlg}.
#' @export
#' @examples
#' \donttest{ # Takes a while on CRAN
#' set.seed(1)
#' n <- 50
#' x <- y <- numeric(n)
#' y[1] <- rnorm(1, exp(x[1]), 0.1)
#' for(i in 1:(n-1)) {
#' x[i+1] <- rnorm(1, sin(x[i]), 0.1)
#' y[i+1] <- rnorm(1, exp(x[i+1]), 0.1)
#' }
#'
#' pntrs <- cpp_example_model("nlg_sin_exp")
#'
#' model_nlg <- ssm_nlg(y = y, a1 = pntrs$a1, P1 = pntrs$P1,
#' Z = pntrs$Z_fn, H = pntrs$H_fn, T = pntrs$T_fn, R = pntrs$R_fn,
#' Z_gn = pntrs$Z_gn, T_gn = pntrs$T_gn,
#' theta = c(log_H = log(0.1), log_R = log(0.1)),
#' log_prior_pdf = pntrs$log_prior_pdf,
#' n_states = 1, n_etas = 1, state_names = "state")
#'
#' out <- ekf(model_nlg, iekf_iter = 100)
#' ts.plot(cbind(x, out$at[1:n], out$att[1:n]), col = 1:3)
#' }
ssm_nlg <- function(y, Z, H, T, R, Z_gn, T_gn, a1, P1, theta,
known_params = NA, known_tv_params = matrix(NA), n_states, n_etas,
log_prior_pdf, time_varying = rep(TRUE, 4),
state_names = paste0("state", 1:n_states)) {
if (is.null(dim(y))) {
dim(y) <- c(length(y), 1)
}
if(missing(n_etas)) {
n_etas <- n_states
}
n_states <- as.integer(n_states)
n_etas <- as.integer(n_etas)
theta <- check_theta(theta)
structure(list(y = as.ts(y), Z = Z, H = H, T = T,
R = R, Z_gn = Z_gn, T_gn = T_gn, a1 = a1, P1 = P1, theta = theta,
log_prior_pdf = log_prior_pdf, known_params = known_params,
known_tv_params = known_tv_params,
n_states = n_states, n_etas = n_etas,
time_varying = time_varying,
state_names = state_names,
max_iter = 100, conv_tol = 1e-8),
class = c("ssm_nlg", "bssm_model"))
}
#'
#' Univariate state space model with continuous SDE dynamics
#'
#' Constructs an object of class \code{ssm_sde} by defining the functions for
#' the drift, diffusion and derivative of diffusion terms of univariate SDE,
#' as well as the log-density of observation equation. We assume that the
#' observations are measured at integer times (missing values are allowed).
#'
#' As in case of \code{ssm_nlg} models, these general models need a bit more
#' effort from the user, as you must provide the several small C++ snippets
#' which define the model structure. See vignettes for an example and
#' \code{cpp_example_model}.
#'
#' @param y Observations as univariate time series (or vector) of length
#' \eqn{n}.
#' @param drift,diffusion,ddiffusion An external pointers for the C++ functions
#' which
#' define the drift, diffusion and derivative of diffusion functions of SDE.
#' @param obs_pdf An external pointer for the C++ function which
#' computes the observational log-density given the the states and parameter
#' vector theta.
#' @param prior_pdf An external pointer for the C++ function which
#' computes the prior log-density given the parameter vector theta.
#' @param theta Parameter vector passed to all model functions.
#' @param x0 Fixed initial value for SDE at time 0.
#' @param positive If \code{TRUE}, positivity constraint is
#' forced by \code{abs} in Milstein scheme.
#' @return An object of class \code{ssm_sde}.
#' @export
#' @examples
#'
#' \donttest{ # Takes a while on CRAN
#' library("sde")
#' set.seed(1)
#' # theta_0 = rho = 0.5
#' # theta_1 = nu = 2
#' # theta_2 = sigma = 0.3
#' x <- sde.sim(t0 = 0, T = 50, X0 = 1, N = 50,
#' drift = expression(0.5 * (2 - x)),
#' sigma = expression(0.3),
#' sigma.x = expression(0))
#' y <- rpois(50, exp(x[-1]))
#'
#' # source c++ snippets
#' pntrs <- cpp_example_model("sde_poisson_OU")
#'
#' sde_model <- ssm_sde(y, pntrs$drift, pntrs$diffusion,
#' pntrs$ddiffusion, pntrs$obs_density, pntrs$prior,
#' c(rho = 0.5, nu = 2, sigma = 0.3), 1, positive = FALSE)
#'
#' est <- particle_smoother(sde_model, L = 12, particles = 500)
#'
#' ts.plot(cbind(x, est$alphahat,
#' est$alphahat - 2*sqrt(c(est$Vt)),
#' est$alphahat + 2*sqrt(c(est$Vt))),
#' col = c(2, 1, 1, 1), lty = c(1, 1, 2, 2))
#'
#'
#' # Takes time with finer mesh, parallelization with IS-MCMC helps a lot
#' out <- run_mcmc(sde_model, L_c = 4, L_f = 8,
#' particles = 50, iter = 2e4,
#' threads = 4L)
#' }
#'
ssm_sde <- function(y, drift, diffusion, ddiffusion, obs_pdf,
prior_pdf, theta, x0, positive) {
y <- check_y(y)
theta <- check_theta(theta)
structure(list(y = as.ts(y), drift = drift,
diffusion = diffusion,
ddiffusion = ddiffusion, obs_pdf = obs_pdf,
prior_pdf = prior_pdf, theta = theta, x0 = x0,
positive = positive, state_names = "x"),
class = c("ssm_sde", "bssm_model"))
}
helske/bssm documentation built on Oct. 29, 2023, 6:04 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ssm_sde ssm_nlg ar1_lg ar1_ng svm bsm_ng bsm_lg ssm_mng ssm_mlg ssm_ung ssm_ulg default_update_fn default_prior_fn
Documented in ar1_lg ar1_ng bsm_lg bsm_ng ssm_mlg ssm_mng ssm_nlg ssm_sde ssm_ulg ssm_ung svm
#' @srrstats {G2.3, G2.3a, G2.3b} match.arg and tolower used where applicable.
#' @srrstats {G2.7, G2.8, G2.9} Only matrix/mts/arrays as tabular data are
#' supported, not data.frame or similar objects.
#' @srrstats {G2.14, G2.14a, G2.14b, G2.14c, BS3.0} Missing observations are
#' handled automatically as per SSM theory, whereas missing values are not
#' allowed elsewhere.
#' @srrstats {BS1.0, BS1.1, BS1.2, BS1.2c} Examples and definitions of priors.
NULL
## placeholder functions for fixed models
default_prior_fn <- function(theta) {
0
}
default_update_fn <- function(theta) {
}
#'
#' General univariate linear-Gaussian state space models
#'
#' Construct an object of class \code{ssm_ulg} by directly defining the
#' corresponding terms of the model.
#'
#' The general univariate linear-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{y_t = D_t + Z_t \alpha_t + H_t \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, 1)}, \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other.
#' Here k is the number of disturbance terms which can be less than m, the
#' number of states.
#'
#' The \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{H} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_ulg},
#' \code{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.
#'
#' @inheritParams ssm_ung
#' @param H A vector of standard deviations. Either a scalar or a vector of
#' length n.
#' @param 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
#' \code{Z}, \code{H}, \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D}, and
#' \code{C}, where each element matches the dimensions of the original model
#' It's best to check the internal dimensions with \code{str(model_object)} as
#' the dimensions of input arguments can differ from the final dimensions.
#' If any of these components is missing, it is assumed to be constant wrt.
#' theta.
#' @return An object of class \code{ssm_ulg}.
#' @export
#' @examples
#'
#' # Regression model with time-varying coefficients
#' set.seed(1)
#' n <- 100
#' x1 <- rnorm(n)
#' x2 <- rnorm(n)
#' b1 <- 1 + cumsum(rnorm(n, sd = 0.5))
#' b2 <- 2 + cumsum(rnorm(n, sd = 0.1))
#' y <- 1 + b1 * x1 + b2 * x2 + rnorm(n, sd = 0.1)
#'
#' Z <- rbind(1, x1, x2)
#' H <- 0.1
#' T <- diag(3)
#' R <- diag(c(0, 1, 0.1))
#' a1 <- rep(0, 3)
#' P1 <- diag(10, 3)
#'
#' # updates the model given the current values of the parameters
#' update_fn <- function(theta) {
#' R <- diag(c(0, theta[1], theta[2]))
#' dim(R) <- c(3, 3, 1)
#' list(R = R, H = theta[3])
#' }
#' # prior for standard deviations as half-normal(1)
#' prior_fn <- function(theta) {
#' if(any(theta < 0)) {
#' log_p <- -Inf
#' } else {
#' log_p <- sum(dnorm(theta, 0, 1, log = TRUE))
#' }
#' log_p
#' }
#'
#' model <- ssm_ulg(y, Z, H, T, R, a1, P1,
#' init_theta = c(1, 0.1, 0.1),
#' update_fn = update_fn, prior_fn = prior_fn,
#' state_names = c("level", "b1", "b2"),
#' # using default values, but being explicit for testing purposes
#' C = matrix(0, 3, 1), D = numeric(1))
#'
#' out <- run_mcmc(model, iter = 5000)
#' out
#' sumr <- summary(out, variable = "state", times = 1:n)
#' sumr$true <- c(b1, b2, rep(1, n))
#' library(ggplot2)
#' ggplot(sumr, aes(x = time, y = Mean)) +
#' geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), alpha = 0.5) +
#' geom_line() +
#' geom_line(aes(y = true), colour = "red") +
#' facet_wrap(~ variable, scales = "free") +
#' theme_bw()
#'
#' # Perhaps easiest way to construct a general SSM for bssm is to use the
#' # model building functionality of KFAS:
#' library("KFAS")
#'
#' model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = 5e-4)+
#' SSMseasonal(period = 12, sea.type = "trigonometric", Q = 0) +
#' log(PetrolPrice) + law, data = Seatbelts, H = 0.005)
#'
#' # use as_bssm function for conversion, kappa defines the
#' # prior variance for diffuse states
#' model_bssm <- as_bssm(model_kfas, kappa = 100)
#'
#' # define updating function for parameter estimation
#' # we can use SSModel and as_bssm functions here as well
#' # (for large model it is more efficient to do this
#' # "manually" by constructing only necessary matrices,
#' # i.e., in this case a list with H and Q)
#'
#' prior_fn <- function(theta) {
#' if(any(theta < 0)) -Inf else sum(dnorm(theta, 0, 0.1, log = TRUE))
#' }
#'
#' update_fn <- function(theta) {
#'
#' model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
#' SSMseasonal(period = 12,
#' sea.type = "trigonometric", Q = theta[2]^2) +
#' log(PetrolPrice) + law, data = Seatbelts, H = theta[3]^2)
#'
#' # the bssm_model object is essentially list so this is fine
#' as_bssm(model_kfas, kappa = 100, init_theta = init_theta,
#' update_fn = update_fn, prior_fn = prior_fn)
#' }
#'
#' init_theta <- rep(1e-2, 3)
#' names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y")
#'
#' model_bssm <- update_fn(init_theta)
#'
#' \donttest{
#' out <- run_mcmc(model_bssm, iter = 10000, burnin = 5000)
#' out
#' }
#' # Above the regression coefficients are modelled as
#' # time-invariant latent states.
#' # Here is an alternative way where we use variable D so that the
#' # coefficients are part of parameter vector theta. Note however that the
#' # first option often preferable in order to keep the dimension of theta low.
#'
#' updatefn2 <- function(theta) {
#' # note no PetrolPrice or law variables here
#' model_kfas2 <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
#' SSMseasonal(period = 12, sea.type = "trigonometric", Q = theta[2]^2),
#' data = Seatbelts, H = theta[3]^2)
#'
#' X <- model.matrix(~ -1 + law + log(PetrolPrice), data = Seatbelts)
#' D <- t(X %*% theta[4:5])
#' as_bssm(model_kfas2, D = D, kappa = 100)
#' }
#' prior2 <- function(theta) {
#' if(any(theta[1:3] < 0)) {
#' -Inf
#' } else {
#' sum(dnorm(theta[1:3], 0, 0.1, log = TRUE)) +
#' sum(dnorm(theta[4:5], 0, 10, log = TRUE))
#' }
#' }
#' init_theta <- c(rep(1e-2, 3), 0, 0)
#' names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y", "law", "Petrol")
#' model_bssm2 <- updatefn2(init_theta)
#' model_bssm2$theta <- init_theta
#' model_bssm2$prior_fn <- prior2
#' model_bssm2$update_fn <- updatefn2
#' \donttest{
#' out2 <- run_mcmc(model_bssm2, iter = 10000, burnin = 5000)
#' out2
#' }
ssm_ulg <- function(y, Z, H, T, R, a1 = NULL, P1 = NULL,
init_theta = numeric(0),
D = NULL, C = NULL, state_names, update_fn = default_update_fn,
prior_fn = default_prior_fn) {
y <- check_y(y)
n <- length(y)
# create Z
Z <- check_Z(Z, 1L, n)
m <- dim(Z)[1]
# create T
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
H <- check_H(H, 1L, n)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
rownames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if (is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
# xreg and beta are need in C++ side in order to combine constructors
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R, a1 = a1, P1 = P1,
D = D, C = C, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
xreg = matrix(0, 0, 0), beta = numeric(0)),
class = c("ssm_ulg", "lineargaussian", "bssm_model"))
}
#' General univariate non-Gaussian state space model
#'
#' Construct an object of class \code{ssm_ung} by directly defining the
#' corresponding terms of the model.
#'
#' The general univariate non-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{p(y_t | D_t + Z_t \alpha_t), (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other,
#' and \eqn{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 \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{phi} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_ung},
#' \code{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.
#'
#' @inheritParams bsm_ng
#' @param y Observations as time series (or vector) of length \eqn{n}.
#' @param 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.
#' @param 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.
#' @param 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.
#' @param a1 Prior mean for the initial state as a vector of length m.
#' @param P1 Prior covariance matrix for the initial state as m x m matrix.
#' @param state_names A character vector defining the names of the states.
#' @param C Intercept terms \eqn{C_t} for the state equation, given as a
#' m times 1 or m times n matrix.
#' @param D Intercept terms \eqn{D_t} for the observations equation, given as a
#' scalar or vector of length n.
#' @param init_theta Initial values for the unknown hyperparameters theta
#' (i.e. unknown variables excluding latent state variables).
#' @param 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
#' \code{Z}, \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D}, \code{C}, and
#' \code{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
#' \code{str(model_object)} as the dimensions of input arguments can differ
#' from the final dimensions.
#' @param prior_fn A function which returns log of prior density
#' given input vector theta.
#' @return An object of class \code{ssm_ung}.
#' @export
#' @examples
#'
#' 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)
ssm_ung <- function(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) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, distribution = distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
Z <- check_Z(Z, 1L, n)
m <- dim(Z)[1]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
check_phi(phi)
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
rownames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
# xreg and beta are need in C++ side in order to combine constructors
structure(list(y = as.ts(y), Z = Z, T = T, R = R, a1 = a1, P1 = P1,
phi = phi, u = u, D = D, C = C, distribution = distribution,
initial_mode = initial_mode, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE,
xreg = matrix(0, 0, 0), beta = numeric(0)),
class = c("ssm_ung", "nongaussian", "bssm_model"))
}
#' General multivariate linear Gaussian state space models
#'
#' Construct an object of class \code{ssm_mlg} by directly defining the
#' corresponding terms of the model.
#'
#' The general multivariate linear-Gaussian model is defined using the
#' following observational and state equations:
#'
#' \deqn{y_t = D_t + Z_t \alpha_t + H_t \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, I_p)}, \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other.
#' Here p is the number of time series and k is the number of disturbance terms
#' (which can be less than m, the number of states).
#'
#' The \code{update_fn} function should take only one
#' vector argument which is used to create list with elements named as
#' \code{Z}, \code{H} \code{T}, \code{R}, \code{a1}, \code{P1}, \code{D},
#' and \code{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 \code{ssm_mlg},
#' \code{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
#'
#' @inheritParams ssm_ulg
#' @param y Observations as multivariate time series or matrix with
#' dimensions n x p.
#' @param Z System matrix Z of the observation equation as p x m matrix or
#' p x m x n array.
#' @param H Lower triangular matrix H of the observation. Either a scalar or
#' a vector of length n.
#' @param T System matrix T of the state equation. Either a m x m matrix or a
#' m x m x n array.
#' @param R Lower triangular matrix R the state equation. Either a m x k matrix
#' or a m x k x n array.
#' @param D Intercept terms for observation equation, given as a p x n matrix.
#' @param C Intercept terms for state equation, given as m x n matrix.
#' @return An object of class \code{ssm_mlg}.
#' @export
#' @examples
#'
#' data("GlobalTemp", package = "KFAS")
#' model_temp <- ssm_mlg(GlobalTemp, H = matrix(c(0.15,0.05,0, 0.05), 2, 2),
#' R = 0.05, Z = matrix(1, 2, 1), T = 1, P1 = 10,
#' state_names = "temperature",
#' # using default values, but being explicit for testing purposes
#' D = matrix(0, 2, 1), C = matrix(0, 1, 1))
#' ts.plot(cbind(model_temp$y, smoother(model_temp)$alphahat), col = 1:3)
#'
ssm_mlg <- function(y, Z, H, T, R, a1 = NULL, P1 = NULL,
init_theta = numeric(0), D = NULL, C = NULL, state_names,
update_fn = default_update_fn, prior_fn = default_prior_fn) {
# create y
y <- check_y(y, multivariate = TRUE)
n <- nrow(y)
p <- ncol(y)
# create Z
Z <- check_Z(Z, p, n, multivariate = TRUE)
m <- dim(Z)[2]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, p, n)
D <- as.matrix(D) # p = 1
C <- check_C(C, m, n)
H <- check_H(H, p, n, multivariate = TRUE)
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
colnames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R, a1 = a1,
P1 = P1, D = D, C = C, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
state_names = state_names), class = c("ssm_mlg", "lineargaussian",
"bssm_model"))
}
#' General Non-Gaussian State Space Model
#'
#' Construct an object of class \code{ssm_mng} by directly defining the
#' corresponding terms of the model.
#'
#' The general multivariate non-Gaussian model is defined using the following
#' observational and state equations:
#'
#' \deqn{p^i(y^i_t | D_t + Z_t \alpha_t), (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = C_t + T_t \alpha_t + R_t \eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\eta_t \sim N(0, I_k)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other, and
#' \eqn{p^i(y_t | .)} is either Poisson, binomial, gamma, Gaussian, or
#' negative binomial distribution for each observation series \eqn{i=1,...,p}.
#' Here k is the number of disturbance terms (which can be less than m,
#' the number of states).
#' @inheritParams ssm_ung
#' @param y Observations as multivariate time series or matrix with dimensions
#' n x p.
#' @param Z System matrix Z of the observation equation as p x m matrix or
#' p x m x n array.
#' @param T System matrix T of the state equation. Either a m x m matrix or a
#' m x m x n array.
#' @param R Lower triangular matrix R the state equation. Either a m x k
#' matrix or a
#' m x k x n array.
#' @param distribution A vector of distributions of the observed series.
#' Possible choices are
#' \code{"poisson"}, \code{"binomial"}, \code{"negative binomial"},
#' \code{"gamma"}, and \code{"gaussian"}.
#' @param 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.
#' @param 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)).
#' @param D Intercept terms for observation equation, given as p x n matrix.
#' @param C Intercept terms for state equation, given as m x n matrix.
#' @return An object of class \code{ssm_mng}.
#' @export
#' @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))
#'
ssm_mng <- function(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) {
# create y
y <- check_y(y, multivariate = TRUE)
n <- nrow(y)
p <- ncol(y)
if (missing(u)) u <- matrix(1, n, p)
u <- check_u(u, y, multivariate = TRUE)
if(length(distribution) == 1) distribution <- rep(distribution, p)
check_distribution(y, distribution)
if(length(phi) == 1) phi <- rep(phi, p)
for(i in 1:p) {
distribution[i] <- match.arg(tolower(distribution[i]),
c("poisson", "binomial", "negative binomial", "gamma", "gaussian"))
check_phi(phi[i])
}
# create Z
Z <- check_Z(Z, p, n, multivariate = TRUE)
m <- dim(Z)[2]
T <- check_T(T, m, n)
# create R
R <- check_R(R, m, n)
a1 <- check_a1(a1, m)
P1 <- check_P1(P1, m)
D <- check_D(D, p, n)
if (p == 1) D <- as.matrix(D)
C <- check_C(C, m, n)
initial_mode <- y
for(i in 1:p) {
initial_mode[, i] <- init_mode(y[, i], u[, i], distribution[i])
}
if (missing(state_names)) {
state_names <- paste("state", 1:m)
}
colnames(Z) <- colnames(T) <- rownames(T) <- rownames(R) <- names(a1) <-
rownames(P1) <- colnames(P1) <- state_names
if(is.null(names(init_theta)) && length(init_theta) > 0)
names(init_theta) <- paste0("theta_", seq_along(init_theta))
structure(list(y = as.ts(y), Z = Z, T = T, R = R, a1 = a1, P1 = P1,
phi = phi, u = u, D = D, C = C, distribution = distribution,
initial_mode = initial_mode, update_fn = update_fn,
prior_fn = prior_fn, theta = init_theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("ssm_mng", "nongaussian", "bssm_model"))
}
#' Basic Structural (Time Series) Model
#'
#' Constructs a basic structural model with local level or local trend
#' component and seasonal component.
#'
#' @inheritParams bsm_ng
#' @param sd_y Standard deviation of the noise of observation equation.
#' Should be an object of class \code{bssm_prior} or scalar
#' value defining a known value such as 0.
#' @param D Intercept terms for observation equation, given as a length n
#' numeric vector or a scalar in case of time-invariant intercept.
#' @param C Intercept terms for state equation, given as a m times n matrix
#' or m times 1 matrix in case of time-invariant intercept.
#' @return An object of class \code{bsm_lg}.
#' @export
#' @examples
#'
#' set.seed(1)
#' n <- 50
#' x <- rnorm(n)
#' level <- numeric(n)
#' level[1] <- rnorm(1)
#' for (i in 2:n) level[i] <- rnorm(1, -0.2 + level[i-1], sd = 0.1)
#' y <- rnorm(n, 2.1 + x + level)
#' model <- bsm_lg(y, sd_y = halfnormal(1, 5), sd_level = 0.1, a1 = level[1],
#' P1 = matrix(0, 1, 1), xreg = x, beta = normal(1, 0, 1),
#' D = 2.1, C = matrix(-0.2, 1, 1))
#'
#' ts.plot(cbind(fast_smoother(model), level), col = 1:2)
#'
#' prior <- uniform(0.1 * sd(log10(UKgas)), 0, 1)
#' # period here is redundant as frequency(UKgas) = 4
#' model_UKgas <- bsm_lg(log10(UKgas), sd_y = prior, sd_level = prior,
#' sd_slope = prior, sd_seasonal = prior, period = 4)
#'
#' # Note small number of iterations for CRAN checks
#' mcmc_out <- run_mcmc(model_UKgas, iter = 5000)
#' summary(mcmc_out, return_se = TRUE)
#' # Use the summary method from coda:
#' summary(expand_sample(mcmc_out, "theta"))$stat
#' mcmc_out$theta[which.max(mcmc_out$posterior), ]
#' sqrt((fit <- StructTS(log10(UKgas), type = "BSM"))$coef)[c(4, 1:3)]
#'
bsm_lg <- function(y, sd_y, sd_level, sd_slope, sd_seasonal,
beta, xreg = NULL, period, a1 = NULL, P1 = NULL, D = NULL,
C = NULL) {
y <- check_y(y)
n <- length(y)
regression_part <- create_regression(beta, xreg, n)
notfixed <- c("y" = 1, "level" = 1, "slope" = 1, "seasonal" = 1)
if (missing(sd_y) || is.null(sd_y)) {
stop("Provide either prior or fixed value for sd_y.")
} else {
if (is_prior(sd_y)) {
check_sd(sd_y$init, "y")
H <- matrix(sd_y$init)
} else {
notfixed[1] <- 0
check_sd(sd_y, "y")
H <- matrix(sd_y)
}
}
if (missing(sd_level) || is.null(sd_level)) {
stop("Provide either prior or fixed value for sd_level.")
} else {
if (is_prior(sd_level)) {
check_sd(sd_level$init, "level")
} else {
notfixed["level"] <- 0
check_sd(sd_level, "level")
}
}
if (missing(sd_slope) || is.null(sd_slope)) {
notfixed["slope"] <- 0
slope <- FALSE
sd_slope <- NULL
} else {
if (is_prior(sd_slope)) {
check_sd(sd_slope$init, "slope")
} else {
notfixed["slope"] <- 0
check_sd(sd_slope, "slope")
}
slope <- TRUE
}
if (missing(sd_seasonal) || is.null(sd_seasonal)) {
notfixed["seasonal"] <- 0
seasonal_names <- NULL
seasonal <- FALSE
sd_seasonal <- NULL
period <- 1L
} else {
if (missing(period)) period <- frequency(y)
period <- check_period(period, n)
if (is_prior(sd_seasonal)) {
check_sd(sd_seasonal$init, "seasonal")
} else {
notfixed["seasonal"] <- 0
check_sd(sd_seasonal, "seasonal")
}
seasonal <- TRUE
seasonal_names <- paste0("seasonal_", 1:(period - 1))
}
npar_R <- 1L + as.integer(slope) + as.integer(seasonal)
m <- as.integer(1L + as.integer(slope) + as.integer(seasonal) * (period - 1))
a1 <- check_a1(a1, m)
if (missing(P1)) {
P1 <- diag(100, m)
} else {
P1 <- check_P1(P1, m)
}
if (slope) {
state_names <- c("level", "slope", seasonal_names)
} else {
state_names <- c("level", seasonal_names)
}
Z <- matrix(0, m, 1)
Z[1, 1] <- 1
if (seasonal) {
Z[2 + slope, 1] <- 1
}
T <- matrix(0, m, m)
T[1, 1] <- 1
if (slope) {
T[1:2, 2] <- 1
}
if (seasonal) {
T[(2 + slope), (2 + slope):m] <- -1
T[cbind(1 + slope + 2:(period - 1), 1 + slope + 1:(period - 2))] <- 1
}
R <- matrix(0, m, max(1, npar_R))
if (notfixed["level"]) {
R[1, 1] <- sd_level$init
} else {
R[1, 1] <- sd_level
}
if (slope) {
if (notfixed["slope"]) {
R[2, 2] <- sd_slope$init
} else {
R[2, 2] <- sd_slope
}
}
if (seasonal) {
if (notfixed["seasonal"]) {
R[2 + slope, 2 + slope] <- sd_seasonal$init
} else {
R[2 + slope, 2 + slope] <- sd_seasonal
}
}
dim(H) <- 1
dim(T) <- c(m, m, 1)
dim(R) <- c(m, ncol(R), 1)
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- state_names
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(sd_y, sd_level, sd_slope, sd_seasonal),
regression_part$beta)
} else {
priors <- list(sd_y, sd_level, sd_slope, sd_seasonal,
regression_part$beta)
}
names(priors) <- c("sd_y", "sd_level", "sd_slope", "sd_seasonal",
names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
D <- check_D(D, 1L, n)
C <- check_C(C, m, n)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R,
a1 = a1, P1 = P1, xreg = regression_part$xreg,
beta = regression_part$coefs, D = D, C = C,
slope = slope, seasonal = seasonal, period = period,
fixed = as.integer(!notfixed),
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta), class = c("bsm_lg", "ssm_ulg", "lineargaussian",
"bssm_model"))
}
#' Non-Gaussian Basic Structural (Time Series) Model
#'
#' 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).
#'
#' @param y A vector or a \code{ts} object of observations.
#' @param sd_level Standard deviation of the noise of level equation.
#' Should be an object of class \code{bssm_prior} or scalar
#' value defining a known value such as 0.
#' @param sd_slope Standard deviation of the noise of slope equation.
#' Should be an object of class \code{bssm_prior}, scalar
#' value defining a known value such as 0, or missing, in which case the slope
#' term is omitted from the model.
#' @param sd_seasonal Standard deviation of the noise of seasonal equation.
#' Should be an object of class \code{bssm_prior}, scalar
#' value defining a known value such as 0, or missing, in which case the
#' seasonal term is omitted from the model.
#' @param sd_noise A prior for the standard deviation of the additional noise
#' term to be added to linear predictor, defined as an object of class
#' \code{bssm_prior}. If missing, no additional noise term is used.
#' @param distribution Distribution of the observed time series. Possible
#' choices are \code{"poisson"}, \code{"binomial"}, \code{"gamma"}, and
#' \code{"negative binomial"}.
#' @param 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 \code{bssm_prior} or
#' a positive scalar.
#' @param 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.
#' @param beta A prior for the regression coefficients.
#' Should be an object of class \code{bssm_prior} or \code{bssm_prior_list}
#' (in case of multiple coefficients) or missing in case of no covariates.
#' @param xreg A matrix containing covariates with number of rows matching the
#' length of \code{y}. Can also be \code{ts}, \code{mts} or similar object
#' convertible to matrix.
#' @param period Length of the seasonal pattern.
#' Must be a positive value greater than 2 and less than the length of the
#' input time series. Default is \code{frequency(y)},
#' which can also return non-integer value (in which case error is given).
#' @param a1 Prior means for the initial states (level, slope, seasonals).
#' Defaults to vector of zeros.
#' @param P1 Prior covariance matrix for the initial states (level, slope,
#' seasonals).Default is diagonal matrix with 100 on the diagonal.
#' @param C Intercept terms for state equation, given as a m x n or m x 1
#' matrix.
#' @return An object of class \code{bsm_ng}.
#' @export
#' @examples
#' # Same data as in Vihola, Helske, Franks (2020)
#' data(poisson_series)
#' s <- sd(log(pmax(0.1, poisson_series)))
#' model <- bsm_ng(poisson_series, sd_level = uniform(0.115, 0, 2 * s),
#' sd_slope = uniform(0.004, 0, 2 * s), P1 = diag(0.1, 2),
#' distribution = "poisson")
#'
#' \donttest{
#' out <- run_mcmc(model, iter = 1e5, particles = 10)
#' summary(out, variable = "theta", return_se = TRUE)
#' # should be about 0.093 and 0.016
#' summary(out, variable = "states", return_se = TRUE,
#' states = 1, times = c(1, 100))
#' # should be about -0.075, 2.618
#' }
#'
#' 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 = 12L,
#' a1 = rep(0, 1 + 11), # level + period - 1 seasonal states
#' P1 = diag(1, 12),
#' C = matrix(0, 12, 1),
#' u = rep(1, nrow(Seatbelts)))
#'
#' \donttest{
#' 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()
#'
#' }
#' # 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))
#'
bsm_ng <- function(y, sd_level, sd_slope, sd_seasonal, sd_noise,
distribution, phi, u, beta, xreg = NULL, period,
a1 = NULL, P1 = NULL, C = NULL) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, multivariate = FALSE, distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
regression_part <- create_regression(beta, xreg, n)
notfixed <- c("level" = 1, "slope" = 1, "seasonal" = 1)
if (missing(sd_level) || missing(sd_level)) {
stop("Provide either prior or fixed value for sd_level.")
} else {
if (is_prior(sd_level)) {
check_sd(sd_level$init, "level")
} else {
notfixed["level"] <- 0
check_sd(sd_level, "level")
}
}
if (missing(sd_slope) || is.null(sd_slope)) {
notfixed["slope"] <- 0
slope <- FALSE
sd_slope <- NULL
} else {
if (is_prior(sd_slope)) {
check_sd(sd_slope$init, "slope")
} else {
notfixed["slope"] <- 0
check_sd(sd_slope, "slope")
}
slope <- TRUE
}
if (missing(sd_seasonal) || is.null(sd_seasonal)) {
notfixed["seasonal"] <- 0
seasonal_names <- NULL
seasonal <- FALSE
sd_seasonal <- NULL
period <- 1L
} else {
if (missing(period)) period <- frequency(y)
period <- check_period(period, n)
if (is_prior(sd_seasonal)) {
check_sd(sd_seasonal$init, "seasonal")
} else {
notfixed["seasonal"] <- 0
check_sd(sd_seasonal, "seasonal")
}
seasonal <- TRUE
seasonal_names <- paste0("seasonal_", 1:(period - 1))
}
if (missing(sd_noise) || is.null(sd_noise)) {
noise <- FALSE
sd_noise <- NULL
} else {
check_sd(sd_noise$init, "noise")
noise <- TRUE
}
npar_R <- 1L + as.integer(slope) + as.integer(seasonal) + as.integer(noise)
m <- as.integer(1L + as.integer(slope) +
as.integer(seasonal) * (period - 1) + as.integer(noise))
a1 <- check_a1(a1, m)
if (missing(P1)) {
P1 <- diag(100, m)
} else {
P1 <- check_P1(P1, m)
}
if (slope) {
state_names <- c("level", "slope", seasonal_names)
} else {
state_names <- c("level", seasonal_names)
}
Z <- matrix(0, m, 1)
Z[1, 1] <- 1
if (seasonal) {
Z[2 + slope, 1] <- 1
}
T <- matrix(0, m, m)
T[1, 1] <- 1
if (slope) {
T[1:2, 2] <- 1
}
if (seasonal) {
T[(2 + slope), (2 + slope):m] <- -1
T[cbind(1 + slope + 2:(period - 1), 1 + slope + 1:(period - 2))] <- 1
}
R <- matrix(0, m, max(1, npar_R))
if (notfixed["level"]) {
R[1, 1] <- sd_level$init
} else {
R[1, 1] <- sd_level
}
if (slope) {
if (notfixed["slope"]) {
R[2, 2] <- sd_slope$init
} else {
R[2, 2] <- sd_slope
}
}
if (seasonal) {
if (notfixed["seasonal"]) {
R[2 + slope, 2 + slope] <- sd_seasonal$init
} else {
R[2 + slope, 2 + slope] <- sd_seasonal
}
}
#additional noise term
if (noise) {
P1[m, m] <- sd_noise$init^2
Z[m] <- 1
state_names <- c(state_names, "noise")
R[m, max(1, ncol(R) - 1)] <- sd_noise$init
}
phi_est <- FALSE
use_phi <- distribution %in% c("negative binomial", "gamma")
if (use_phi) {
if (is_prior(phi)) {
check_phi(phi$init)
phi_est <- TRUE
} else {
check_phi(phi)
}
} else {
phi <- 1
}
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
dim(T) <- c(m, m, 1)
dim(R) <- c(m, ncol(R), 1)
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- state_names
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(sd_level, sd_slope, sd_seasonal, sd_noise, phi),
regression_part$beta)
} else {
priors <- list(sd_level, sd_slope, sd_seasonal, sd_noise, phi,
regression_part$beta)
}
names(priors) <- c("sd_level", "sd_slope", "sd_seasonal", "sd_noise", "phi",
names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
if (phi_est) {
phi <- phi$init
}
D <- numeric(1)
C <- check_C(C, m, n)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = phi, u = u, xreg = regression_part$xreg,
beta = regression_part$coefs, D = D, C = C,
slope = slope, seasonal = seasonal, noise = noise,
period = period, fixed = as.integer(!notfixed),
distribution = distribution, initial_mode = initial_mode,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta, phi_est = phi_est,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("bsm_ng", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Stochastic Volatility Model
#'
#' Constructs a simple stochastic volatility model with Gaussian errors and
#' first order autoregressive signal. See the main vignette for details.
#'
#' @param y A numeric vector or a \code{\link{ts}} object of observations.
#' @param mu A prior for mu parameter of transition equation.
#' Should be an object of class \code{bssm_prior}.
#' @param rho A prior for autoregressive coefficient.
#' Should be an object of class \code{bssm_prior}.
#' @param sd_ar A prior for the standard deviation of noise of the AR-process.
#' Should be an object of class \code{bssm_prior}.
#' @param sigma A prior for sigma parameter of observation equation, internally
#' denoted as phi. Should be an object of class \code{bssm_prior}.
#' Ignored if \code{mu} is provided. Note that typically
#' parametrization using mu is preferred due to better numerical properties and
#' availability of better Gaussian approximation.
#' Most notably the global approximation approach does not work with sigma
#' parameterization as sigma is not a parameter of the resulting approximate
#' model.
#' @return An object of class \code{svm}.
#' @export
#' @rdname svm
#' @examples
#'
#' data("exchange")
#' y <- exchange[1:100] # for faster CRAN check
#' model <- svm(y, rho = uniform(0.98, -0.999, 0.999),
#' sd_ar = halfnormal(0.15, 5), sigma = halfnormal(0.6, 2))
#'
#' obj <- function(pars) {
#' -logLik(svm(y,
#' rho = uniform(pars[1], -0.999, 0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' sigma = halfnormal(pars[3], 2)), particles = 0)
#' }
#' opt <- optim(c(0.98, 0.15, 0.6), obj,
#' lower = c(-0.999, 1e-4, 1e-4),
#' upper = c(0.999, 10, 10), method = "L-BFGS-B")
#' pars <- opt$par
#' model <- svm(y,
#' rho = uniform(pars[1],-0.999,0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' sigma = halfnormal(pars[3], 2))
#'
#' # alternative parameterization
#' model2 <- svm(y, rho = uniform(0.98,-0.999, 0.999),
#' sd_ar = halfnormal(0.15, 5), mu = normal(0, 0, 1))
#'
#' obj2 <- function(pars) {
#' -logLik(svm(y,
#' rho = uniform(pars[1], -0.999, 0.999),
#' sd_ar = halfnormal(pars[2], 5),
#' mu = normal(pars[3], 0, 1)), particles = 0)
#' }
#' opt2 <- optim(c(0.98, 0.15, 0), obj2, lower = c(-0.999, 1e-4, -Inf),
#' upper = c(0.999, 10, Inf), method = "L-BFGS-B")
#' pars2 <- opt2$par
#' model2 <- svm(y,
#' rho = uniform(pars2[1],-0.999,0.999),
#' sd_ar = halfnormal(pars2[2], 5),
#' mu = normal(pars2[3], 0, 1))
#'
#' # sigma is internally stored in phi
#' ts.plot(cbind(model$phi * exp(0.5 * fast_smoother(model)),
#' exp(0.5 * fast_smoother(model2))), col = 1:2)
#'
svm <- function(y, mu, rho, sd_ar, sigma) {
if(!missing(sigma) && !missing(mu)) {
stop("Define either sigma or mu, but not both.")
}
y <- check_y(y)
check_rho(rho$init)
check_sd(sd_ar$init, "rho")
if(!missing(mu)) {
svm_type <- 1L
check_mu(mu$init)
initial_mode <- matrix(log(pmax(1e-4, y^2)), ncol = 1)
} else {
svm_type <- 0L
check_sd(sigma$init, "sigma", FALSE)
initial_mode <-
matrix(log(pmax(1e-4, y^2)) - 2 * log(sigma$init), ncol = 1)
}
a1 <- if(svm_type) mu$init else 0
P1 <- matrix(sd_ar$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sd_ar$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
priors <- list(rho, sd_ar, if(svm_type==0) sigma else mu)
priors <- priors[!vapply(priors, is.null, TRUE)]
names(priors) <-
c("rho", "sd_ar", if(svm_type==0) "sigma" else "mu")
C <- if (svm_type) matrix(mu$init * (1 - T[1])) else matrix(0)
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = if (svm_type == 0) sigma$init else 1,
xreg = matrix(0, 0, 0),
beta = numeric(0), D = D, C = C,
initial_mode = initial_mode,
svm_type = svm_type, distribution = "svm", u = 1,
phi_est = !as.logical(svm_type),
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters,
theta = theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("svm", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Non-Gaussian model with AR(1) latent process
#'
#' Constructs a simple non-Gaussian model where the state dynamics follow an
#' AR(1) process.
#'
#' @inheritParams bsm_ng
#' @param rho A prior for autoregressive coefficient.
#' Should be an object of class \code{bssm_prior}.
#' @param mu A fixed value or a prior for the stationary mean of the latent
#' AR(1) process. Should be an object of class \code{bssm_prior} or scalar
#' value defining a fixed mean such as 0.
#' @param sigma A prior for the standard deviation of noise of the AR-process.
#' Should be an object of class \code{bssm_prior}
#' @return An object of class \code{ar1_ng}.
#' @export
#' @rdname ar1_ng
#' @examples
#' model <- ar1_ng(discoveries, rho = uniform(0.5,-1,1),
#' sigma = halfnormal(0.1, 1), mu = normal(0, 0, 1),
#' distribution = "poisson")
#' out <- run_mcmc(model, iter = 1e4, mcmc_type = "approx",
#' output_type = "summary")
#'
#' ts.plot(cbind(discoveries, exp(out$alphahat)), col = 1:2)
#'
#' set.seed(1)
#' n <- 30
#' phi <- 2
#' rho <- 0.9
#' sigma <- 0.1
#' beta <- 0.5
#' u <- rexp(n, 0.1)
#' x <- rnorm(n)
#' z <- y <- numeric(n)
#' z[1] <- rnorm(1, 0, sigma / sqrt(1 - rho^2))
#' y[1] <- rnbinom(1, mu = u * exp(beta * x[1] + z[1]), size = phi)
#' for(i in 2:n) {
#' z[i] <- rnorm(1, rho * z[i - 1], sigma)
#' y[i] <- rnbinom(1, mu = u * exp(beta * x[i] + z[i]), size = phi)
#' }
#'
#' model <- ar1_ng(y, rho = uniform_prior(0.9, 0, 1),
#' sigma = gamma_prior(0.1, 2, 10), mu = 0.,
#' phi = gamma_prior(2, 2, 1), distribution = "negative binomial",
#' xreg = x, beta = normal_prior(0.5, 0, 1), u = u)
#'
ar1_ng <- function(y, rho, sigma, mu, distribution, phi, u, beta,
xreg = NULL) {
distribution <- match.arg(tolower(distribution),
c("poisson", "binomial", "negative binomial", "gamma"))
y <- check_y(y, multivariate = FALSE, distribution)
n <- length(y)
if (missing(u)) u <- rep(1, n)
u <- check_u(u, y)
regression_part <- create_regression(beta, xreg, n)
check_rho(rho$init)
check_sd(sigma$init, "rho")
if (is_prior(mu)) {
check_mu(mu$init)
mu_est <- TRUE
a1 <- mu$init
C <- matrix(mu$init * (1 - rho$init))
} else {
mu_est <- FALSE
check_mu(mu)
a1 <- mu
C <- matrix(mu * (1 - rho$init))
}
distribution <- match.arg(tolower(distribution), c("poisson", "binomial",
"negative binomial", "gamma"))
use_phi <- distribution %in% c("negative binomial", "gamma")
phi_est <- FALSE
if (use_phi) {
if (is_prior(phi)) {
check_phi(phi$init)
phi_est <- TRUE
} else {
check_phi(phi)
}
} else {
phi <- 1
}
initial_mode <- matrix(init_mode(y, u, distribution), ncol = 1)
P1 <- matrix(sigma$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sigma$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(rho, sigma, mu, phi), regression_part$beta)
} else {
priors <- list(rho, sigma, mu, phi, regression_part$beta)
}
names(priors) <-
c("rho", "sigma", "mu", "phi", names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
if (phi_est) {
phi <- phi$init
}
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, T = T, R = R,
a1 = a1, P1 = P1, phi = phi, u = u,
xreg = regression_part$xreg, beta = regression_part$coefs,
D = D, C = C,
initial_mode = initial_mode,
distribution = distribution, mu_est = mu_est, phi_est = phi_est,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters, theta = theta,
max_iter = 100, conv_tol = 1e-8, local_approx = TRUE),
class = c("ar1_ng", "ssm_ung", "nongaussian", "bssm_model"))
}
#' Univariate Gaussian model with AR(1) latent process
#'
#' Constructs a simple Gaussian model where the state dynamics
#' follow an AR(1) process.
#'
#' @inheritParams ar1_ng
#' @param sd_y A prior for the standard deviation of observation equation.
#' @return An object of class \code{ar1_lg}.
#' @export
#' @rdname ar1_lg
#' @examples
#' set.seed(1)
#' mu <- 2
#' rho <- 0.7
#' sd_y <- 0.1
#' sigma <- 0.5
#' beta <- -1
#' x <- rnorm(30)
#' z <- y <- numeric(30)
#' z[1] <- rnorm(1, mu, sigma / sqrt(1 - rho^2))
#' y[1] <- rnorm(1, beta * x[1] + z[1], sd_y)
#' for(i in 2:30) {
#' z[i] <- rnorm(1, mu * (1 - rho) + rho * z[i - 1], sigma)
#' y[i] <- rnorm(1, beta * x[i] + z[i], sd_y)
#' }
#' model <- ar1_lg(y, rho = uniform(0.5, -1, 1),
#' sigma = halfnormal(1, 10), mu = normal(0, 0, 1),
#' sd_y = halfnormal(1, 10),
#' xreg = x, beta = normal(0, 0, 1))
#' out <- run_mcmc(model, iter = 2e4)
#' summary(out, return_se = TRUE)
#'
ar1_lg <- function(y, rho, sigma, mu, sd_y, beta, xreg = NULL) {
y <- check_y(y)
n <- length(y)
regression_part <- create_regression(beta, xreg, n)
check_rho(rho$init)
check_sd(sigma$init, "rho")
if (is_prior(mu)) {
check_mu(mu$init)
mu_est <- TRUE
a1 <- mu$init
C <- matrix(mu$init * (1 - rho$init))
} else {
mu_est <- FALSE
check_mu(mu)
a1 <- mu
C <- matrix(mu * (1 - rho$init))
}
if (is_prior(sd_y)) {
check_sd(sd_y$init, "y")
sd_y_est <- TRUE
H <- matrix(sd_y$init)
} else {
sd_y_est <- FALSE
check_sd(sd_y, "y")
H <- matrix(sd_y)
}
P1 <- matrix(sigma$init^2 / (1 - rho$init^2))
Z <- matrix(1)
T <- array(rho$init, c(1, 1, 1))
R <- array(sigma$init, c(1, 1, 1))
names(a1) <- rownames(P1) <- colnames(P1) <- rownames(Z) <-
rownames(T) <- colnames(T) <- rownames(R) <- "signal"
if(ncol(regression_part$xreg) > 1) {
priors <- c(list(rho, sigma, mu, sd_y), regression_part$beta)
} else {
priors <- list(rho, sigma, mu, sd_y, regression_part$beta)
}
names(priors) <-
c("rho", "sigma", "mu", "sd_y", names(regression_part$coefs))
priors <- priors[vapply(priors, is_prior, TRUE)]
D <- matrix(0)
theta <- if (length(priors) > 0) {
vapply(priors, "[[", "init", FUN.VALUE = 1)
} else numeric(0)
priors <- combine_priors(priors)
structure(list(y = as.ts(y), Z = Z, H = H, T = T, R = R,
a1 = a1, P1 = P1,
xreg = regression_part$xreg, beta = regression_part$coefs,
D = D, C = C,
mu_est = mu_est, sd_y_est = sd_y_est,
prior_distributions = priors$prior_distribution,
prior_parameters = priors$parameters, theta = theta,
max_iter = 100, conv_tol = 1e-8),
class = c("ar1_lg", "ssm_ulg", "lineargaussian", "bssm_model"))
}
#'
#' General multivariate nonlinear Gaussian state space models
#'
#' Constructs an object of class \code{ssm_nlg} by defining the corresponding
#' terms of the observation and state equation.
#'
#' The nonlinear Gaussian model is defined as
#'
#' \deqn{y_t = Z(t, \alpha_t, \theta) + H(t, \theta) \epsilon_t,
#' (\textrm{observation equation})}
#' \deqn{\alpha_{t+1} = T(t, \alpha_t, \theta) + R(t, \theta)\eta_t,
#' (\textrm{transition equation})}
#'
#' where \eqn{\epsilon_t \sim N(0, I_p)}, \eqn{\eta_t \sim N(0, I_m)} and
#' \eqn{\alpha_1 \sim N(a_1, P_1)} independently of each other, and functions
#' \eqn{Z, H, T, R} can depend on \eqn{\alpha_t} and parameter vector
#' \eqn{\theta}.
#'
#' Compared to other models, these general models need a bit more effort from
#' the user, as you must provide the several small C++ snippets which define the
#' model structure. See examples in the vignette and \code{cpp_example_model}.
#'
#' @param y Observations as multivariate time series (or matrix) of length
#' \eqn{n}.
#' @param Z,H,T,R An external pointers (object of class \code{externalptr})
#' for the C++ functions which define the corresponding model functions.
#' @param Z_gn,T_gn An external pointers (object of class \code{externalptr})
#' for the C++ functions which define the gradients of the corresponding model
#' functions.
#' @param a1 Prior mean for the initial state as object of class
#' \code{externalptr}
#' @param P1 Prior covariance matrix for the initial state as object of class
#' \code{externalptr}
#' @param theta Parameter vector passed to all model functions.
#' @param known_params A vector of known parameters passed to all model
#' functions.
#' @param known_tv_params A matrix of known parameters passed to all model
#' functions.
#' @param n_states Number of states in the model (positive integer).
#' @param n_etas Dimension of the noise term of the transition equation
#' (positive integer).
#' @param log_prior_pdf An external pointer (object of class
#' \code{externalptr}) for the C++ function which
#' computes the log-prior density given theta.
#' @param time_varying Optional logical vector of length 4, denoting whether
#' the values of
#' Z, H, T, and R vary with respect to time variable (given identical states).
#' If used, this can speed up some computations.
#' @param state_names A character vector containing names for the states.
#' @return An object of class \code{ssm_nlg}.
#' @export
#' @examples
#' \donttest{ # Takes a while on CRAN
#' set.seed(1)
#' n <- 50
#' x <- y <- numeric(n)
#' y[1] <- rnorm(1, exp(x[1]), 0.1)
#' for(i in 1:(n-1)) {
#' x[i+1] <- rnorm(1, sin(x[i]), 0.1)
#' y[i+1] <- rnorm(1, exp(x[i+1]), 0.1)
#' }
#'
#' pntrs <- cpp_example_model("nlg_sin_exp")
#'
#' model_nlg <- ssm_nlg(y = y, a1 = pntrs$a1, P1 = pntrs$P1,
#' Z = pntrs$Z_fn, H = pntrs$H_fn, T = pntrs$T_fn, R = pntrs$R_fn,
#' Z_gn = pntrs$Z_gn, T_gn = pntrs$T_gn,
#' theta = c(log_H = log(0.1), log_R = log(0.1)),
#' log_prior_pdf = pntrs$log_prior_pdf,
#' n_states = 1, n_etas = 1, state_names = "state")
#'
#' out <- ekf(model_nlg, iekf_iter = 100)
#' ts.plot(cbind(x, out$at[1:n], out$att[1:n]), col = 1:3)
#' }
ssm_nlg <- function(y, Z, H, T, R, Z_gn, T_gn, a1, P1, theta,
known_params = NA, known_tv_params = matrix(NA), n_states, n_etas,
log_prior_pdf, time_varying = rep(TRUE, 4),
state_names = paste0("state", 1:n_states)) {
if (is.null(dim(y))) {
dim(y) <- c(length(y), 1)
}
if(missing(n_etas)) {
n_etas <- n_states
}
n_states <- as.integer(n_states)
n_etas <- as.integer(n_etas)
theta <- check_theta(theta)
structure(list(y = as.ts(y), Z = Z, H = H, T = T,
R = R, Z_gn = Z_gn, T_gn = T_gn, a1 = a1, P1 = P1, theta = theta,
log_prior_pdf = log_prior_pdf, known_params = known_params,
known_tv_params = known_tv_params,
n_states = n_states, n_etas = n_etas,
time_varying = time_varying,
state_names = state_names,
max_iter = 100, conv_tol = 1e-8),
class = c("ssm_nlg", "bssm_model"))
}
#'
#' Univariate state space model with continuous SDE dynamics
#'
#' Constructs an object of class \code{ssm_sde} by defining the functions for
#' the drift, diffusion and derivative of diffusion terms of univariate SDE,
#' as well as the log-density of observation equation. We assume that the
#' observations are measured at integer times (missing values are allowed).
#'
#' As in case of \code{ssm_nlg} models, these general models need a bit more
#' effort from the user, as you must provide the several small C++ snippets
#' which define the model structure. See vignettes for an example and
#' \code{cpp_example_model}.
#'
#' @param y Observations as univariate time series (or vector) of length
#' \eqn{n}.
#' @param drift,diffusion,ddiffusion An external pointers for the C++ functions
#' which
#' define the drift, diffusion and derivative of diffusion functions of SDE.
#' @param obs_pdf An external pointer for the C++ function which
#' computes the observational log-density given the the states and parameter
#' vector theta.
#' @param prior_pdf An external pointer for the C++ function which
#' computes the prior log-density given the parameter vector theta.
#' @param theta Parameter vector passed to all model functions.
#' @param x0 Fixed initial value for SDE at time 0.
#' @param positive If \code{TRUE}, positivity constraint is
#' forced by \code{abs} in Milstein scheme.
#' @return An object of class \code{ssm_sde}.
#' @export
#' @examples
#'
#' \donttest{ # Takes a while on CRAN
#' library("sde")
#' set.seed(1)
#' # theta_0 = rho = 0.5
#' # theta_1 = nu = 2
#' # theta_2 = sigma = 0.3
#' x <- sde.sim(t0 = 0, T = 50, X0 = 1, N = 50,
#' drift = expression(0.5 * (2 - x)),
#' sigma = expression(0.3),
#' sigma.x = expression(0))
#' y <- rpois(50, exp(x[-1]))
#'
#' # source c++ snippets
#' pntrs <- cpp_example_model("sde_poisson_OU")
#'
#' sde_model <- ssm_sde(y, pntrs$drift, pntrs$diffusion,
#' pntrs$ddiffusion, pntrs$obs_density, pntrs$prior,
#' c(rho = 0.5, nu = 2, sigma = 0.3), 1, positive = FALSE)
#'
#' est <- particle_smoother(sde_model, L = 12, particles = 500)
#'
#' ts.plot(cbind(x, est$alphahat,
#' est$alphahat - 2*sqrt(c(est$Vt)),
#' est$alphahat + 2*sqrt(c(est$Vt))),
#' col = c(2, 1, 1, 1), lty = c(1, 1, 2, 2))
#'
#'
#' # Takes time with finer mesh, parallelization with IS-MCMC helps a lot
#' out <- run_mcmc(sde_model, L_c = 4, L_f = 8,
#' particles = 50, iter = 2e4,
#' threads = 4L)
#' }
#'
ssm_sde <- function(y, drift, diffusion, ddiffusion, obs_pdf,
prior_pdf, theta, x0, positive) {
y <- check_y(y)
theta <- check_theta(theta)
structure(list(y = as.ts(y), drift = drift,
diffusion = diffusion,
ddiffusion = ddiffusion, obs_pdf = obs_pdf,
prior_pdf = prior_pdf, theta = theta, x0 = x0,
positive = positive, state_names = "x"),
class = c("ssm_sde", "bssm_model"))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.