bssm (version 1.0.0)

predict.mcmc_output: Predictions for State Space Models

Description

Draw samples from the posterior predictive distribution given the posterior draws of hyperparameters theta and alpha_n+1.

Usage

# S3 method for mcmc_output
predict(
  object,
  future_model,
  type = "response",
  seed = sample(.Machine$integer.max, size = 1),
  nsim,
  ...
)

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

Return predictions on "mean" "response", or "state" level.

seed

Seed for RNG.

nsim

Number of samples to draw.

...

Ignored.

Value

Data frame of predicted samples.

Examples

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

mcmc_results <- run_mcmc(model, iter = 5000)
future_model <- model
future_model$y <- ts(rep(NA, 25), 
  start = tsp(model$y)[2] + 2 * deltat(model$y), 
  frequency = frequency(model$y))
pred <- predict(mcmc_results, future_model, type = "state", 
  nsim = 1000)

require("dplyr")
sumr_fit <- as.data.frame(mcmc_results, variable = "states") %>%
  group_by(time, iter) %>% 
  mutate(signal = 
      value[variable == "level"] + 
      value[variable == "seasonal_1"]) %>%
  group_by(time) %>%
  summarise(mean = mean(signal), 
    lwr = quantile(signal, 0.025), 
    upr = quantile(signal, 0.975))

sumr_pred <- pred %>% 
  group_by(time, sample) %>%
  mutate(signal = 
      value[variable == "level"] + 
      value[variable == "seasonal_1"]) %>%
  group_by(time) %>%
  summarise(mean = mean(signal),
    lwr = quantile(signal, 0.025), 
    upr = quantile(signal, 0.975)) 

require("ggplot2")
rbind(sumr_fit, sumr_pred) %>% 
  ggplot(aes(x = time, y = mean)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), 
   fill = "#92f0a8", alpha = 0.25) +
  geom_line(colour = "#92f0a8") +
  theme_bw() + 
  geom_point(data = data.frame(
    mean = log10(JohnsonJohnson), 
    time = time(JohnsonJohnson)))

# }