lynx.fun <- function(tsb) {
ar.fit <- ar(tsb, order.max = 25)
c(ar.fit$order, mean(tsb), tsb)
}
# the stationary bootstrap with mean block length 20
lynx.1 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "geom")
# the fixed block bootstrap with length 20
lynx.2 <- tsboot(log(lynx), lynx.fun, R = 99, l = 20, sim = "fixed")
# Now for model based resampling we need the original model
# Note that for all of the bootstraps which use the residuals as their
# data, we set orig.t to FALSE since the function applied to the residual
# time series will be meaningless.
lynx.ar <- ar(log(lynx))
lynx.model <- list(order = c(lynx.ar$order, 0, 0), ar = lynx.ar$ar)
lynx.res <- lynx.ar$resid[!is.na(lynx.ar$resid)]
lynx.res <- lynx.res - mean(lynx.res)
lynx.sim <- function(res,n.sim, ran.args) {
# random generation of replicate series using arima.sim
rg1 <- function(n, res) sample(res, n, replace = TRUE)
ts.orig <- ran.args$ts
ts.mod <- ran.args$model
mean(ts.orig)+ts(arima.sim(model = ts.mod, n = n.sim,
rand.gen = rg1, res = as.vector(res)))
}
lynx.3 <- tsboot(lynx.res, lynx.fun, R = 99, sim = "model", n.sim = 114,
orig.t = FALSE, ran.gen = lynx.sim,
ran.args = list(ts = log(lynx), model = lynx.model))
# For "post-blackening" we need to define another function
lynx.black <- function(res, n.sim, ran.args) {
ts.orig <- ran.args$ts
ts.mod <- ran.args$model
mean(ts.orig) + ts(arima.sim(model = ts.mod,n = n.sim,innov = res))
}
# Now we can run apply the two types of block resampling again but this
# time applying post-blackening.
lynx.1b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "fixed",
n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,
ran.args = list(ts = log(lynx), model = lynx.model))
lynx.2b <- tsboot(lynx.res, lynx.fun, R = 99, l = 20, sim = "geom",
n.sim = 114, orig.t = FALSE, ran.gen = lynx.black,
ran.args = list(ts = log(lynx), model = lynx.model))
# To compare the observed order of the bootstrap replicates we
# proceed as follows.
table(lynx.1$t[, 1])
table(lynx.1b$t[, 1])
table(lynx.2$t[, 1])
table(lynx.2b$t[, 1])
table(lynx.3$t[, 1])
# Notice that the post-blackened and model-based bootstraps preserve
# the true order of the model (11) in many more cases than the others.
Run the code above in your browser using DataLab