library(ts)
data(lynx)
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