Description Usage Arguments Value Examples
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).
1 2 3 4 5 6 7 8 9 10 
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  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)))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.