# This R code pertains to the time series analysis of
# Sec. 8.6 of E&T.

# I'll read in the 48 level values from p. 92 of E&T.
# (They are in a file I created.) 
data <- read.table("http://mason.gmu.edu/~csutton/EandTT81.txt", header=FALSE)
level <- data[,1]

# If I use blocks of 3 consecutive observations, as is done in E&T, then I
# just need to resample 16 times in order to get a time series of length 48.
# I'll resample from the integers 1 through 46 (since a block cannot start 
# at the 47th or 48th observation and still be of length 3.
blockstart <- c(1:46)
set.seed(321)
drawnblockstart <- sample( blockstart, size=200*16, replace=T )

# Now I'll convert the 200*16 blockstarts into 200*48 level values,
# and then put them into a 200 by 48 matrix.
levelval <- numeric()
for (k in 1:3200) levelval[(3*k-2):(3*k)] <- level[drawnblockstart[k]:(drawnblockstart[k]+2)]
bss <- matrix( levelval, nrow=200, byrow=T )
# Note: It's very important to have byrow=T in order to make sure the level
# values in levelval are put into the matrix by rows.  The default is to put 
# by columns, which would leave the time series in the rows lacking the proper
# correlation structure. 

# Now I can obtain the desired bootstrap replicates.
bhrep <- numeric()
row.mean <- apply( bss, 1, mean )
for (i in 1:200) {
  ztstar <- bss[i,2:48] - row.mean[i] 
  ztminus1star <- bss[i,1:47] - row.mean[i]
  bhrep[i] <- sum( ztstar*ztminus1star )/sum( ztminus1star^2 )
 }

# Here's the bootstrap estimate of the standard error of the 
# estimator of beta:
sd(bhrep)

# Here's the mean of the bootstrap replicates and a histogram:
mean( bhrep )
hist( bhrep )

# With a small proportion of the replicates being as large as the
# estimate of beta based on the original time series, it can be
# concluded that the bootstrapped time series tend to be not as
# strongly correlated as the original series.  This makes sense,
# as the first-order correlation will be weak where the blocks end
# (every third observation).  If blocks of size 12 are used the 
# bootstrap time series ought to resemble the original series more.
# I'll modify what I have above and check this out below. 
blockstart <- c(1:37)
drawnblockstart <- sample( blockstart, size=200*4, replace=T )
for (k in 1:800) levelval[(12*k-11):(12*k)] <- level[drawnblockstart[k]:(drawnblockstart[k]+11)]
bss <- matrix( levelval, nrow=200, byrow=T )
bhrep <- numeric()
row.mean <- apply( bss, 1, mean )
for (i in 1:200) {
  ztstar <- bss[i,2:48] - row.mean[i] 
  ztminus1star <- bss[i,1:47] - row.mean[i]
  bhrep[i] <- sum( ztstar*ztminus1star )/sum( ztminus1star^2 )
 }

# Here's the bootstrap estimate of the standard error of the 
# estimator of beta:
sd(bhrep)

# Here's the mean of the bootstrap replicates and a histogram:
mean( bhrep )
hist( bhrep )