Bootstrap bias of the sample correlation I am working on a problem asking
"Provide a Bootstrap bias of $\hat{r}$ using R = 1000 bootstrap estimators." However, I am stuck on what even to do. As far as I know, the formula for finding bootstrap bias is
$Bias_{boot}(\hat{r}) = E[\hat{r}^\ast] - \hat{r} $
Please correct me if I am wrong. So from here, I used simple bootstrap function as
boot.fn = function (data, index) {
    X = data[index]
    Y = data[index]
    sample_cor_hat = cor(X, Y)
  return(sample_cor_hat)
}

b = boot(sample_distribution, boot.fn, R = 1000)
b

Then I would get
Bootstrap Statistics :
    original       bias     std. error
t1*        1 1.110223e-16 7.719711e-17

Is my formula for bias correct? am I using the correct way for the function?
update: here's the given data
sample_size = 20
sample_meanvector = c(2, 3)
sample_covariance_matrix = matrix(c(2, 0.4, 0.4, 1), ncol = 2)
sample_distribution = mvrnorm(n = sample_size, mu = sample_meanvector, Sigma = sample_covariance_matrix)
X = matrix(sample_distribution[,1])
Y = matrix(sample_distribution[,2])


It  seems you are using a simple quantile style of
bootstrap. Perhaps it would be better to use bias-corrected
quantile bootstrap as below.
Roughly, each bootstrap
re-sample of the sample correlation is compared to the
observed sample correlation, and the distances d.re between
them are found. Then the boostrap CI is based on quantiles .025 and .975 of the $B$ d.res.
For my fictitious normal data x and y, this bootstrap CI and the classical CI for population $\rho$ are essentially
both $(-0.14, 0.25.).$
# Fictitious normal data and classical CI
set.seed(2022)
x = rnorm(100, 50, 6)
y = rnorm(100, 50, 6)
r.obs = cor(x,y)
cor.test(x,y)$conf.int
[1] -0.1409677  0.2506400
attr(,"conf.level")
[1] 0.95
plot(x, y, pch=19)


# 95% Nonparametric bootstrap CI
set.seed(1234)
B= 2000;  seq = 1:100;  d.re=numeric(B)
for (i in 1:B)  {
 seq.re = sample(seq, 100,rep=T)
 d.re[i] = cor(x[seq.re],y[seq.re]) - r.obs
 }
LR = quantile(d.re, c(.975,.025))
r.obs - LR 
     97.5%       2.5% 
-0.1387333  0.2473008 

Then for non normal data,
the same bootstrap method provides a reasonable
95% CI for $\rho.$
set.seed(608)
x = rnorm(100, 50, 6) + 1:100
y = rnorm(100, 50, 6) + 1:100
r.obs = cor(x,y);  r.obs
[1] 0.9559053
plot(x, y, pch=19)
shapiro.test(x)$p.val
[1] 0.002321988       # Not normal
shapiro.test(y)$p.val
[1] 0.01022498        # Not normal

Because x and y are not normal one cannot rely
on the CI from cor.test; see end note. However, my bootstrap procedure should give useful CIs for samples as large as 100.

set.seed(1234)
B= 2000;  seq = 1:100;  d.re=numeric(B)
for (i in 1:B)  {
 seq.re = sample(seq, 100,rep=T)
 d.re[i] = cor(x[seq.re],y[seq.re]) - r.obs
 }
LR = quantile(d.re, c(.975,.025))
r.obs - LR
    97.5%      2.5% 
0.9412820 0.9751636 

Note from R documentation for cor.test: "If method is 'pearson', the test statistic is based on Pearson's product moment correlation coefficient cor(x, y) and follows a t distribution with length(x)-2 degrees of freedom if the samples follow independent normal distributions. [Emphasis added.]
