B <- 199
nMC <- 1000
coverage <- rep(0,nMC)
for (m in 1:nMC){
nobs <- 1000
y <- rnorm(nobs, 0, 1)
mu_bootstrapped <- rep(0,B)
for (b in 1:B){
y_b <- sample(y, replace = TRUE)
mu_bootstrapped[b] <- sd(y_b)
}
ub <- sd(y) + 1.96 * sd(mu_bootstrapped)
lb <- sd(y) - 1.96 * sd(mu_bootstrapped)
coverage[m] <- (ub > 1) * (lb < 1)
}
mean(coverage)
## [1] 0.943