Predictions for State Space Models
Draw samples from the posterior predictive distribution for future time points given the posterior draws of hyperparameters θ and alpha_{n+1}. Function can also be used to draw samples from the posterior predictive distribution p(\tilde y_1, …, \tilde y_n | y_1,…, y_n).
## S3 method for class 'mcmc_output' predict( object, model, type = "response", nsim, future = TRUE, seed = sample(.Machine$integer.max, size = 1), ... )
object |
mcmc_output object obtained from
|
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.
It is also possible to input the original model, which can be useful for example for
posterior predictive checks. In this case, set argument |
type |
Return predictions on |
nsim |
Number of samples to draw. |
future |
Default is |
seed |
Seed for RNG. |
... |
Ignored. |
Data frame of predicted samples.
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))
# use "state" for illustrative purposes, we could use type = "mean" directly
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))
# If we used type = "mean", we could do
# sumr_pred <- pred %>%
# group_by(time) %>%
# summarise(mean = mean(value),
# lwr = quantile(value, 0.025),
# upr = quantile(value, 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)))
# Posterior predictions for past observations:
yrep <- predict(mcmc_results, model, type = "response",
future = FALSE, nsim = 1000)
meanrep <- predict(mcmc_results, model, type = "mean",
future = FALSE, nsim = 1000)
sumr_yrep <- yrep %>%
group_by(time) %>%
summarise(earnings = mean(value),
lwr = quantile(value, 0.025),
upr = quantile(value, 0.975)) %>%
mutate(interval = "Observations")
sumr_meanrep <- meanrep %>%
group_by(time) %>%
summarise(earnings = mean(value),
lwr = quantile(value, 0.025),
upr = quantile(value, 0.975)) %>%
mutate(interval = "Mean")
rbind(sumr_meanrep, sumr_yrep) %>%
mutate(interval = factor(interval, levels = c("Observations", "Mean"))) %>%
ggplot(aes(x = time, y = earnings)) +
geom_ribbon(aes(ymin = lwr, ymax = upr, fill = interval),
alpha = 0.75) +
theme_bw() +
geom_point(data = data.frame(
earnings = model$y,
time = time(model$y)))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.