bssm (version 0.1.11)

predict.mcmc_output: Predictions for State Space Models

Description

Posterior intervals of future observations or their means (success probabilities in binomial case). These are computed using the quantile method where the intervals are computed as empirical quantiles the posterior sample, or using a parametric method by Helske (2016) in a linear-Gaussian case.

Usage

# S3 method for mcmc_output
predict(object, future_model, type = "response",
  intervals = TRUE, probs = c(0.05, 0.95), nsim, return_MCSE = FALSE,
  seed = sample(.Machine$integer.max, size = 1), ...)

Arguments

object

mcmc_output object obtained from run_mcmc

future_model

Model for future observations. Should have same structure as the original model which was used in MCMC, in order to plug the posterior samples of the model parameters to the right places.

type

Compute predictions on "mean" ("confidence interval"), "response" ("prediction interval"), or "state" level. Defaults to "response".

intervals

If TRUE, intervals are returned. Otherwise samples from the posterior predictive distribution are returned.

probs

Desired quantiles. Defaults to c(0.05, 0.95). Always includes median 0.5.

nsim

Number of state samples to draw per MCMC iteration. Note that this has no effect for the time point $n+1$ (where $n$ is the length of the original series) as this is directly obtained from the MCMC output. nsim defaults to 1 except for the EKF based MCMC output of non-linear Gaussian models (see below). For linear-Gaussian models the intervals are computed based on Kalman filter so this argument has no effect if intervals is TRUE. For non-linear Gaussian models of class nlg_ssm, if nsim is 0 and intervals is TRUE, EKF based approximation is used for computing the prediction intervals.

return_MCSE

For Gaussian models, if TRUE, the Monte Carlo standard errors of the intervals are also returned.

seed

Seed for RNG.

...

Ignored.

Value

List containing the mean predictions, quantiles and Monte Carlo standard errors .

Examples

# NOT RUN {
require("graphics")
y <- log10(JohnsonJohnson)
prior <- uniform(0.01, 0, 1)
model <- bsm(window(y, end = c(1974, 4)), sd_y = prior,
sd_level = prior, sd_slope = prior, sd_seasonal = prior)

mcmc_results <- run_mcmc(model, n_iter = 5000)
future_model <- model
future_model$y <- ts(rep(NA, 25), start = end(model$y), 
  frequency = frequency(model$y))
pred_gaussian <- predict(mcmc_results, future_model, 
  probs = c(0.05, 0.1, 0.5, 0.9, 0.95))
ts.plot(log10(JohnsonJohnson), pred_gaussian$intervals, 
  col = c(1, rep(2, 5)), lty = c(1, 2, 2, 1, 2, 2))

head(pred_gaussian$intervals)
head(pred_gaussian$MCSE)

# Non-gaussian models
# }
# NOT RUN {
data("poisson_series")

model <- ng_bsm(poisson_series, sd_level = 
  halfnormal(0.1, 1), sd_slope=halfnormal(0.01, 0.1),
  distribution = "poisson")
mcmc_poisson <- run_mcmc(model, n_iter = 5000, nsim = 10)

future_model <- model
future_model$y <- ts(rep(NA, 25), start = end(model$y), 
  frequency = frequency(model$y))

pred <- predict(mcmc_poisson, future_model, 
  probs = seq(0.05,0.95, by = 0.05))

library("ggplot2")
fit <- ts(colMeans(exp(expand_sample(mcmc_poisson, 
    "alpha")$level)))
autoplot(pred, y = model$y, fit = fit)
# }