Bootstrap Standard Errors: should I divide the sampling standard deviation by $\sqrt{n}$? Suppose I am bootstrapping an OLS regression and want the standard error of the coefficient $\beta_1$. I estimate the following regression on 1000 resamples of the data (where $B$ indexes the bootstrap resample): $$y^B_i=\beta^B_0+\beta^B_1 x^B_i+\varepsilon^B_i$$
Is the standard error of $\beta_1$ given by  $sd(\beta^B_1)$ or $\frac{sd(\beta^B_1)}{\sqrt{n-1}}$ (where sd is the standard deviation)?
It seems correct to divide by $\sqrt{n-1}$ since that is how we usually define a standard error, but I have seen many references that instruct me to only take the standard deviation without dividing by $\sqrt{n-1}$.
 A: The standard error of the mean, with its $\sqrt{n-1}$ in the denominator, is a particular closed form of the standard error (SE), obtained exactly for the mean as a considered statistic. For any statistic, SE is defined as the standard deviation (SD) of its sampling distribution.
Thus, for the bootstrap, SE is the SD of statistic values obtained for bootstrap samples. In your case:
$$se=sd(\beta_1^B)$$ not the $$\frac{sd(\beta_1^B)}{\sqrt{n-1}}$$
Regarding the SE for mean, recall that $Var(aX) = a^2Var(X)$, where $a$ is some constant. Also, there is $Var(\sum_i X_i)$ = $\sum_i Var(X_i)$ for independent $X_i$. Therefore, we got:
$$Var(\bar X ) = Var(\frac{1}{n}\sum_i^n X_i) = \frac{1}{n^2}\sum_i^nVar(X_i) =\frac{1}{n^2}n\sigma^2 = \frac{\sigma^2}{n}$$
where $\sigma$ is the variance of population from each $X_i$ is drawn.
How is this related to sampling distribution?
We can safely stop here and state that whenever there is a sampling of population by taking $n$ elements of it and statistic of interest  is the mean of those elements, the variance of that mean is given by $\frac{\sigma^2}{n}$. So, the well-known form for the SE of mean is:
$$se_{mean}=\frac{\sigma}{\sqrt n}$$ and with the Bessel's correction: $$se_{mean}=\frac{\sigma}{\sqrt{n-1}}$$ .
But you also can find the following explanation. Let's denote vector of even sub-samples of the investigated population by $X_s = (X_{s_1}, X_{s_2}, ..., X_{s_k})$. If the population mean is estimated as total of all variates in $X_s$ divided by a total number of variates $n$, we got:
$$\bar X = \frac{1}{mk}\sum_{i,j}^{m,k} X_i,{s_j} = \frac{1}{n}\sum_{i}^{n} X^s_i$$
and reasoning is the same as for $Var(\frac{1}{n}\sum_i^n X_i$). However, this is a shortcut and conceptually is different from considering the sampling distribution of the $\bar X$. Also, this approach will get nasty for $Var(\bar X)$ if we estimate $\bar X$ as the mean of means from sub-samples: $\bar X =\frac{1}{n} \sum^n_i \bar X_{s_i}$.
A: According to Bradley & Efron if you want to estimate the standard error $se(\hat\theta)$ by the sample standard deviation of the $B$ replications you will take $$\hat{se}_{B} = \sum_{b=1}^{B}\left([\hat\theta^{\star}(b)-\hat\theta^{\star}(\cdot)]^2 / (B-1)\right)^{1/2}$$
where $\hat\theta^{\star}(\cdot) =\sum_{b=1}^{B}\hat\theta^{\star}(b)/B $
Doing so in R :

#Regression Data
#===============
n <- 27
x <- c(99, 152, 293, 155, 196, 53, 184, 171, 52, 376, 385, 402, 29, 76, 296, 151, 177, 209, 119, 188, 115, 88, 58, 49, 150, 107, 125)
y <- c(25.8, 20.5, 14.3, 23.2, 20.6, 31.1, 20.9, 20.9, 30.4, 16.3, 11.6, 11.8, 32.5, 32.0, 18.0, 24.1, 26.5, 25.8, 28.8, 22.0, 29.7, 28.9, 32.8, 32.5, 25.4, 31.7, 28.5)


#OLS Regression estimates
#========================
ones <- rep(1,n)
DM <- cbind(ones,x)
betas.OLS <- solve(t(DM)%*%DM)%*%t(DM)%*%y
betas.OLS

ehat.OLS <- y-DM%*%betas.OLS
ehat.OLS
sigma2.hat <- t(ehat.OLS)%*%ehat.OLS/(n-2)
sigma2.hat
sigma.hat <- sqrt(sigma2.hat)
sigma.hat
var.betas.OLS <- drop(sigma2.hat)*solve(t(DM)%*%DM)
var.betas.OLS
std.betas.OLS <- sqrt(diag(var.betas.OLS))
std.betas.OLS

betas.OLS
std.betas.OLS

#Nonparametric Bootstrap
#=======================
B<-10000
sample.betas<-NULL

for (i in 1:B) 
{
  Boot.residual.sample <- sample(ehat.OLS,size=n,replace=T)
  Y.sample <- DM%*%betas.OLS+Boot.residual.sample
  Boot.betas <- solve(t(DM)%*%DM)%*%t(DM)%*%Y.sample
  sample.betas <- cbind(sample.betas,Boot.betas)
}
sampling.mean <- apply(sample.betas,1,mean)
sampling.mean
sampling.var <- apply(sample.betas,1,var)
sampling.s.e <- sqrt(sampling.var)
sampling.s.e



#Parametric Bootstrap
#====================
B<-10000
sample.betas<-NULL

for (i in 1:B) 
{
  Boot.residual.sample <- rnorm(n,0,sigma.hat)
  Y.sample <- DM%*%betas.OLS+Boot.residual.sample
  Boot.betas <- solve(t(DM)%*%DM)%*%t(DM)%*%Y.sample
  sample.betas <- cbind(sample.betas,Boot.betas)
}
sampling.mean <- apply(sample.betas,1,mean)
sampling.mean
sampling.var <- apply(sample.betas,1,var)
sampling.s.e <- sqrt(sampling.var)
sampling.s.e
```

