# Monte Carlo comparison of performances of different confidence intervals
# for the estimation setting pertaining to Table 13.3 on p. 175 of E&T.

set.seed(321)

# I'll initialize the results matrix.
results <- matrix( rep( 0, 9 ), nrow=3 )
rownames(results) <- c("standard","percentile","bootstrap t")
colnames(results) <- c("miss left","miss right","miss total") 
results

# I'll go through a loop 2500 times.  On each pass through the loop I'll
# generate a new data vector, compute 3 c.i's, and for each one determine 
# it it misses the estimand by falling to the left of it or falling to the 
# right of it.
for (i in 1:2500) {

  z <- rnorm(10)

  theta <- function(x, xdata, ...) { exp( mean( x ) ) }
  ci.result <- boott( z, theta, VS=TRUE, v.nbootg=150, v.nbootsd=200, v.nboott=999, perc=c(0.05,0.95) )
  if ( ci.result$confpoints[1] > 1 ) { 
    results[3,2] <- results[3,2] + 1 
   } 
  else {
    if ( ci.result$confpoints[2] < 1 ) { results[3,1] <- results[3,1] + 1 } 
   }

  bss <- matrix( sample( z, 10*999, replace=T ), nrow=999 )
  mean.rep <- apply( bss, 1, mean )
  theta.rep <- exp( mean.rep )
  sorted.theta.rep <- sort( theta.rep )
  if ( sorted.theta.rep[50] > 1 ) { 
    results[2,2] <- results[2,2] + 1 
   } 
  else {
    if ( sorted.theta.rep[950] < 1 ) { results[2,1] <- results[2,1] + 1 } 
   }

  theta.hat <- exp( mean( z ) )
  se.hat <- sd( theta.rep )
  if ( theta.hat - 1.6449*se.hat > 1 ) { 
    results[1,2] <- results[1,2] + 1 
   } 
  else {
    if ( theta.hat + 1.6449*se.hat < 1 ) { results[1,1] <- results[1,1] + 1 } 
   }
 
 }

results[,3] <- results[,1] + results[,2]
results/2500