TL;DR The article that you refer to makes things look more worse than they actually are. Their bootstrapping procedure is not a good way to apply bootstrapping. In the case of OLS there shouldn't be big problems with high dimensionality if the sample size is large. If you can not get correct results with OLS, where a correct confidence interval can be easily computed analytically, then something must be wrong with the implementation of the bootstrapping method.
It is good though to be reminded that the residuals are not the same as the errors and that we can use simulations with OLS to test (potentially wrong) implementations of bootstrapping.
Simple reproduction of the article results
The article that you refer to is performing simulations of errors by bootstrapping/resampling of the residuals. Below is a simple example that reproduces this.
The model is a linear regression with $n=500$ samples (or 250 pairs) and $p=125$ parameters. The distributions that are plotted here are just for the first parameter estimate $\hat{\beta}_1$.
Discrepancy in estimated sample variance
The third image, resampling the true errors, gives a correct indication of the sample distribution of the coefficient.
The first and second images, resampling all residuals, or resampling the pairs, have distributions with a different variance. They lead to errors in the estimates of standard errors and confidence intervals.
The reason for the discrepancy is that bootstrapping only works when the bootstrapped samples are a good representation of the true distribution. This is not the case when $p/n$ is large.
resampling residuals The bootstrap samples are created by simulating errors by sampling from the residuals, however the variance of the residuals is lower than the variance of the errors $$\text{Var}(r_i) \approx \left(1-\frac{p}{n}\right) \text{Var}(\epsilon)$$
pairwise resampling in the case of pairwise resampling the distribution is effectively a scaled binomial distribution. The variance will be
$$\text{Var}(r_{i,paired}) \approx \frac{1}{2} \text{Var}(\epsilon)$$
The factor $0.5$ stems from the Bernoulli variable having variance $0.25$ and the scaling is the difference of the two pairs which has variance $2\text{Var}(\epsilon)$
Discussion
You should only use bootstrapping when the sampling is done from a distribution that represents the population distribution. This is not the case with paired resampling, which is sampling a distribution with half the variance, and with the residual resampling, which is sampling a distribution with smaller variance if $p/n$ is large.
The bootstrapping is often performed when a distribution of is difficult to compute. This is either the case when 1) the assumptions about the error distribution are false, or the error distribution is unknown 2) the propagation of errors is difficult to compute.
For ordinary linear regression, the second case is not an issue. The statistic is a linear sum of the data and it's sampling distribution will often approximate a normal distribution. With different cost functions the behaviour might not be too far. The problem is just to estimate the variance, and the residuals are often a good indication for this. But, one has to apply the right corrections.
The problem is more difficult in the situation where the error distribution has large tails and the variance is not easily estimated with a small sample. In this case the typical remedy is to simply gather more data. Potentially one could do an advanced semi-parametric bootstrapping by combining the residuals with a normal distribution that relates to the residuals being the errors with the estimate subtracted (the estimate being some correlated normal distribution).
Plot of reproduction
set.seed(1)
n = 500 # data samples
p = 125 # parameters
m = 1000 # times resampling
### create paired data
X = matrix(rnorm(n*p/2),ncol =p)
X = rbind(X,X)
Y = rnorm(n)
solve(t(X)%*%X)[1,1] ### this is the theoretic variance
### compute main model
mod = lm(Y~X+0)
### variables used for resampling
Y_m = predict(mod)
res = mod$residuals
err = Y
### perform resampling of residuals
b_residuals = sapply(1:m, FUN = function(i) {
Y_s = Y_m + sample(res,n)
lm(Y_s~X+0)$coefficients[1]
})
### perform resampling of errors
b_errors = sapply(1:m, FUN = function(i) {
Y_s = Y_m + sample(err,n)
lm(Y_s~X+0)$coefficients[1]
})
### perform paired resampling
b_paired = sapply(1:m, FUN = function(i) {
selection = rep(1:(n/2),2)+rbinom(n,1,0.5)*n/2
Y_s = Y_m + res[selection]
lm(Y_s~X+0)$coefficients[1]
})
### plot histograms
layout(matrix(1:3,3))
hist(b_residuals, breaks = seq(-0.5,0.5,0.02),
freq = 0, ylim = c(0,10), main = "resampling of residuals", xlab = expression(beta[1]))
lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)
hist(b_paired, breaks = seq(-0.5,0.5,0.02),
freq = 0, ylim = c(0,10), main = "resampling of residual pairs", xlab = expression(beta[1]))
lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)
hist(b_errors, breaks = seq(-0.5,0.5,0.02),
freq = 0, ylim = c(0,10), main = "resampling of true errors", xlab = expression(beta[1]))
lines(mod$coefficients[1]*c(1,1),c(0,10), lty = 2, col = 2)
var(b_residuals)/var(b_errors)
var(b_paired)/var(b_errors)
