# NOT RUN {
## 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 <- AddSemilocalLinearTrend(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 <- AddSemilocalLinearTrend(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 {
# }
# 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")
# }
# NOT RUN {
## Example 8: Non-Gaussian data
## Poisson counts of shark attacks in Florida.
data(shark)
logshark <- log1p(shark$Attacks)
ss.level <- AddLocalLevel(list(), y = logshark)
model <- bsts(shark$Attacks, ss.level, niter = 500,
ping = 250, family = "poisson", seed = 8675309)
## Poisson data can have an 'exposure' as the second column of a
## two-column matrix.
model <- bsts(cbind(shark$Attacks, shark$Population / 1000),
state.specification = ss.level, niter = 500,
family = "poisson", ping = 250, seed = 8675309)
# }
Run the code above in your browser using DataLab