predict.mcmc_output {bssm}  R Documentation 
Draw samples from the posterior predictive distribution for future time points given the posterior draws of hyperparameters θ and latent state 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, nsim, type = "response", future = TRUE, seed = sample(.Machine$integer.max, size = 1), ... )
object 
Results object of class 
model 
A 
nsim 
Positive integer defining number of samples to draw. 
type 
Type of predictions. Possible choices are

future 
Default is 
seed 
Seed for RNG (positive integer). Note that this affects only the
C++ side, and 
... 
Ignored. 
A data.frame
consisting of samples from the predictive
posterior distribution.
library("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) library("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)) library("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)))