bulls <- read.csv("http://www.stat.ufl.edu/~winner/sta4702/mj_bulls.csv")
attach(bulls); names(bulls)

X <- bulls[,-c(4,9)]
X <- as.matrix(X)
(N <- nrow(X))

## Obtain Sigma and Rho, based on the population of games


lambda <- eigen()$value    ## Population eigenvalues



## Set up sampling and result holders
set.seed(xxxx) 
num.sim <-               ## number of samples          
n <-                     ## sample size for each sample
k <-                     ## of eigenvalues to save per sample
lambda.hat <- 
    matrix(rep(0,num.sim*k), ncol=k)  ## Matrix of saved eigenvalues

## Begin sampling
for (i1 in 1:num.sim) {
    games <- sample(N, n, replace=F)         # Sample n integers from 1:N
    lambda.s <- eigen(cov(X[games,]))$value   # Obtain eigenvalues of sample S
       for (i2 in 1:k) {
          lambda.hat[i1,i2] <- lambda.s[i2]
       }
}

l1.hat <- lambda.hat[,1]
l2.hat <- lambda.hat[,2]

## Construct 95% CI's for lambda1 
Bon.z <- qnorm(1-.05/(2*k),0,1)
lambda1.LB <-                 ## Obtain Lower Confidence Limit for lambda1
lambda1.UB <-                 ## Obtain Upper Confidence Limit for lambda1

## Construct 95% CI's for lambda2 
lambda2.LB <-                 ## Obtain Lower Confidence Limit for lambda2
lambda2.UB <-                 ## Obtain Upper Confidence Limit for lambda2

## Obtain Coverage probabilities
cover1 <- 
   sum((lambda1.LB <= lambda[1]) & (lambda1.UB >= lambda[1])) / num.sim
cover2 <- 
   sum((lambda2.LB <= lambda[2]) & (lambda2.UB >= lambda[2])) / num.sim
cover12 <- 
   sum((lambda1.LB <= lambda[1]) & (lambda1.UB >= lambda[1]) &
       (lambda2.LB <= lambda[2]) & (lambda2.UB >= lambda[2])) / num.sim

l.CI.out <- cbind(lambda.S.val[1],lambda.S.val[2],mean(l1.hat),mean(l2.hat),
             cover1,cover2,cover12)
colnames(l.CI.out) <- c("Lambda 1", "Lambda 2", "Mean l1.hat", "Mean l2.hat",
     "Cover L1", "Cover L2", "Cover L1 & L2")
round(l.CI.out,3)