## Example 1: Time series (ts) data
data(AirPassengers)
y <- log(AirPassengers)
ss <- AddLocalLinearTrend(list(), y)
ss <- AddSeasonal(ss, y, nseasons = 12)
model <- bsts(y, state.specification = ss, niter = 500)
pred <- predict(model, horizon = 12, burn = 100)
par(mfrow = c(1,2))
plot(model)
plot(pred)
## Not run: ------------------------------------
#
# MakePlots <- function(model, ask = TRUE) {
# ## Make all the plots callable by plot.bsts.
# opar <- par(ask = ask)
# on.exit(par(opar))
# plot.types <- c("state", "components", "residuals",
# "prediction.errors", "forecast.distribution")
# for (plot.type in plot.types) {
# plot(model, plot.type)
# }
# if (model$has.regression) {
# regression.plot.types <- c("coefficients", "predictors", "size")
# for (plot.type in regression.plot.types) {
# plot(model, plot.type)
# }
# }
# }
#
# ## Example 2: GOOG is the Google stock price, an xts series of daily
# ## data.
# data(goog)
# ss <- AddGeneralizedLocalLinearTrend(list(), goog)
# model <- bsts(goog, state.specification = ss, niter = 500)
# MakePlots(model)
#
# ## Example 3: Change GOOG to be zoo, and not xts.
# goog <- zoo(goog, index(goog))
# ss <- AddGeneralizedLocalLinearTrend(list(), goog)
# model <- bsts(goog, state.specification = ss, niter = 500)
# MakePlots(model)
#
# ## Example 4: Naked numeric data works too
# y <- rnorm(100)
# ss <- AddLocalLinearTrend(list(), y)
# model <- bsts(y, state.specification = ss, niter = 500)
# MakePlots(model)
#
# ## Example 5: zoo data with intra-day measurements
# y <- zoo(rnorm(100),
# seq(from = as.POSIXct("2012-01-01 7:00 EST"), len = 100, by = 100))
# ss <- AddLocalLinearTrend(list(), y)
# model <- bsts(y, state.specification = ss, niter = 500)
# MakePlots(model)
#
# \dontrun{
# ## Example 6: Including regressors
# data(iclaims)
# ss <- AddLocalLinearTrend(list(), initial.claims$iclaimsNSA)
# ss <- AddSeasonal(ss, initial.claims$iclaimsNSA, nseasons = 52)
# model <- bsts(iclaimsNSA ~ ., state.specification = ss, data =
# initial.claims, niter = 1000)
# plot(model)
# plot(model, "components")
# plot(model, "coefficients")
# plot(model, "predictors")
# }
## ---------------------------------------------
## Not run: ------------------------------------
# ## Example 7: Regressors with multiple time stamps.
# number.of.time.points <- 50
# sample.size.per.time.point <- 10
# total.sample.size <- number.of.time.points * sample.size.per.time.point
# sigma.level <- .1
# sigma.obs <- 1
#
# ## Simulate some fake data with a local level state component.
# trend <- cumsum(rnorm(number.of.time.points, 0, sigma.level))
# predictors <- matrix(rnorm(total.sample.size * 2), ncol = 2)
# colnames(predictors) <- c("X1", "X2")
# coefficients <- c(-10, 10)
# regression <- as.numeric(predictors <!-- %*% coefficients) -->
# y.hat <- rep(trend, sample.size.per.time.point) + regression
# y <- rnorm(length(y.hat), y.hat, sigma.obs)
#
# ## Create some time stamps, with multiple observations per time stamp.
# first <- as.POSIXct("2013-03-24")
# dates <- seq(from = first, length = number.of.time.points, by = "month")
# timestamps <- rep(dates, sample.size.per.time.point)
#
# ## Run the model with a local level trend, and an unnecessary seasonal component.
# ss <- AddLocalLevel(list(), y)
# ss <- AddSeasonal(ss, y, nseasons = 7)
# model <- bsts(y ~ predictors, ss, niter = 250, timestamps = timestamps,
# seed = 8675309)
# plot(model)
# plot(model, "components")
## ---------------------------------------------
Run the code above in your browser using DataLab