Bootstrap estimation of bias

> setwd("C:/Cauldron/garage/Volatility/Learn/IntroToBootstrap")
> library(bootstrap)
> data(patch)
> mus <- apply(patch, 2, mean)
> theta <- mus[6]/mus[5]
> estimator <- function(B) {
+     n <- dim(patch)[1]
+     theta.hat.star <- numeric()
+     probs <- matrix(data = NA, nrow = B, ncol = n)
+     set.seed(1977)
+     for (i in 1:B) {
+         temp <- sample(1:n, n, T)
+         probs[i, ] <- sapply(1:n, function(x) {
+             sum(temp == x)
+         })/n
+         temp.data <- patch[temp, ]
+         mus.t <- apply(temp.data, 2, mean)
+         theta.t <- mus.t[6]/mus.t[5]
+         theta.hat.star[i] <- theta.t
+     }
+     prob.mu <- colMeans(probs)
+     bias.B.old <- mean(theta.hat.star) - theta
+     bias.B.better <- mean(theta.hat.star) - sum(prob.mu * patch[,
+         6])/sum(prob.mu * patch[, 5])
+     c(bias.B.old, bias.B.better)
+ }
> B <- c(25, 50, 100, 200, 400, 800, 1670, 3200)
> res <- sapply(B, estimator)
> res <- t(res)

> par(mfrow = c(1, 1))
> ylim <- range(res)
> plot(B, res[, 1], type = "l", log = "y", ylim = ylim, col = "red",
+     lwd = 2)
> par(new = T)
> plot(B, res[, 2], type = "l", log = "y", ylim = ylim, col = "green",
+     lwd = 2)
> abline(h = bias.inf, col = "grey", lty = "dashed", lwd = 2)
> legend("topright", legend = c("bootstrap bias estimate ", "better bootstrap bias estimate",
+     "bias limiting value"), fill = c("red", "green", "grey"))