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