################################################################ # # 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].