Everything down to the bottom of the for loop should
be familiar, just like what we did calculating bootstrap standard errors
(sd(theta.star) would be the bootstrap standard error).

The R function quantile calculates quantiles of a data
vector, at least, what it calls quantiles. Its definition is a bit eccentric,
but is asymptotically equivalent to all other definitions of quantiles.

We only look at all these quantiles because Efron and Tibshirani did.
Usually, you only want two. For example,

quantile(theta.star, probs = c(0.025, 0.975))

for a 95% bootstrap percentile confidence interval or

quantile(theta.star, probs = c(0.05, 0.95))

for a 90% bootstrap percentile confidence interval.

An alternative method for quantiles preferred by your humble instructor
uses the following logic.

Use nboot <- 999 (or some other
value such that nboot + 1 is a round number. The reason is
that if X_{(i)} is the i-th order
statistic from a Uniform(0, 1) distribution

E{X_{(i)}} = i / (n + 1)

Another way to think of this is that the nboot data points
divide the number line into nboot + 1 intervals, which as
far as we know contain equal probability. They don't contain equal
probability because the sample is not the population, but we might as well
treat them as such for the purposes of estimation. That is,
our nboot data points should be taken as estimators of the
quantiles with denominators nboot + 1

In particular, if nboot is 999, then we take
the ordered theta.star values to be the 0.001, 0.002, . . .,
0.999 quantiles of the sampling distribution of theta.hat.
Thus

Comments

If we wanted the quantiles to print with probabilities attached we
would have to do

fred <- sort(theta.star)[(nboot + 1) * probs]
names(fred) <- probs
fred

But, of course, we don't want to do that unless we are writing
a textbook like Efron and Tibshirani. Usually we just want the endpoints
of a interval, either

sort(theta.star)[(nboot + 1) * c(0.025, 0.975)]

for a 95% bootstrap percentile confidence interval or

sort(theta.star)[(nboot + 1) * c(0.05, 0.95)]

for a 90% bootstrap percentile confidence interval.