R
bootstrap replicates of a statistic applied to a time series. The
replicate time series can be generated using fixed or random block lengths or
can be model based replicates.tsboot(tseries, statistic, R, l=NULL, sim="model", endcorr=TRUE,
n.sim=length(tseries), orig.t=TRUE, ran.gen=NULL,
ran.args=NULL, norm=TRUE, ...)
rts
or cts
(S-Plus version 3.2 and later) or the earlier
ts
function. Irregular time series, with class "its"
, mtseries
returns a vector containing the
statistic(s) of interest. Each time statistic
is called it is passed
a time series of length n.sim
which is of the same class as the original
"boot"
with the following components.orig.t
is TRUE
then t0
is the result of statistic(tseries,...
)NULL
.tsboot
sim
is "fixed"
then each replicate time series is found by taking
blocks of length l
, from the original time series and putting them
end-to-end until a new series of length n.sim
is created. When sim
is
"geom"
a similar approach is taken except that now the block lengths are
generated from a geometric distribution with mean l
. Post-blackening can
be carried out on these replicate time series by including the function
ran.gen
in the call to tsboot
and having tseries
as a time series of
residuals. Model based resampling is very similar to the parametric bootstrap and all simulation must be in one of the user specified functions. This avoids the complicated problem of choosing the block length but relies on an accurate model choice being made.
Phase scrambling is described in Section 8.2.4 of Davison and Hinkley (1997).
The types of statistic for which this method produces reasonable results is
very limited and the other methods seem to do better in most situations.
Other types of resampling in the frequency domain
can be accomplished using the function boot
with the argument
sim="parametric"
.
Kunsch, H.R. (1989) The jackknife and the bootstrap for general stationary observations. Annals of Statistics, 17, 1217--1241.
Politis, D.N. and Romano, J.P. (1994) The stationary bootstrap. Journal of the American Statistical Association, 89, 1303--1313.
boot
, arima.sim
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