################################################################
#
# The following is an R implementation of the moving blocks bootstrap
# following Efron and Tibshirani (sec. 8.6).
#
# We create some artificial data:
N <- 50 # length of the time series
series <- rnorm(N) # initially noise
series[-1] <- series[-1] + series[-N] # introduce auto-correlation
#
# Here is the code that collects bootstrap values of
# the auto-correlation estimate:
k <- 6 # size of moving blocks
nrep <- 1000 # number of bootstrap replications
lag.cor.bt <- rep(NA,nrep) # vessel for the boostrapped values
for(irep in 1:nrep) { # the bootstrap loop
series.bt <- rep(NA,N) # local vector for a bootstrap replication
for(i in 1:ceiling(N/k)) { # fill the vector with random blocks
endpoint <- sample(k:N, size=1) # by randomly sampling endpoints
series.bt[(i-1)*k+1:k] <- series[endpoint-(k:1)+1] # and copying blocks
}
series.bt <- series.bt[1:N] # trim overflow when k doesn't divide N
lag.cor.bt[irep] <- cor(series.bt[-1],series.bt[-N]) # the auto-cor. estimate
if(irep%%50==0) print(irep) # report every 50 bootstrap passes
}
#
# Now that we have a bootstrap distribution for the statistic of interest
# in lag.cor.bt, we analyze it.
# First visually:
hist(lag.cor.bt, col="gray", ncl=20)
#
# Then numerically:
# 1) the significance level in terms of the quantile of the value
# that represents a null assumption (zero in this case):
sum(lag.cor.bt<0)/1000
# ------------------ output begin
[1] 0.004
# ------------------ output end
# Zero is the 0.4% lower quantile, therefore, the autocorrelation
# is significant at the 1% level (2*0.4% = 0.8% < 1%).
#
# 2) a 95% quantile confidence interval:
quantile(lag.cor.bt, c(0.025,0.975))
# ------------------ output begin
2.5% 97.5%
0.096313 0.601392
# ------------------ output end
# With mild rounding the 95% confidence interval is [0.1,0.6].