1. Introduction

The bootstrap, introduced by Efron (1979), has become the default method for constructing confidence intervals in applied statistics, econometrics, and the experimental sciences. Its appeal is that it requires no parametric assumptions: one resamples with replacement from the observed data, computes the statistic of interest on each resample, and uses the resulting distribution to form confidence intervals. Theoretical justification rests on the bootstrap consistently estimating the sampling distribution of the statistic, which for the sample mean requires at minimum a finite variance (Athreya, 1987).

Heavy-tailed data arise in finance (Embrechts, Kluppelberg, and Mikosch, 1997), insurance, telecommunications, hydrology, and city-size distributions. A distribution $F$ is said to have a regularly varying right tail with index $\alpha > 0$ if

$$\bar{F}(x) = P(X > x) = L(x) , x^{-\alpha}, \quad x \to \infty,$$

where $L(\cdot)$ is a slowly varying function. The parameter $\alpha$ is called the tail index. For $\alpha < 2$ the variance is infinite, and the central limit theorem in its classical form fails; the sample mean is attracted to a stable law with index $\alpha$ rather than to a Gaussian (Nolan, 2020). Athreya (1987) proved that the bootstrap distribution of the re-centered and rescaled sample mean does not converge weakly to the correct stable limit when $\alpha < 2$, and Knight (1989) extended this to show that bootstrap confidence intervals have coverage converging to a value strictly less than the nominal level. Arcones and Gine (1989) generalized these inconsistency results to bootstrap samples of arbitrary size.

While these theoretical results delineate $\alpha = 2$ as the critical boundary, they are asymptotic statements that do not quantify how quickly coverage deteriorates at finite sample sizes, nor do they address the practically important range $2 < \alpha < 3$ where the variance is finite but the fourth moment is infinite. Practitioners working with samples of a few hundred observations from distributions with $\alpha \approx 2.5$ have no concrete guidance on whether their bootstrap intervals can be trusted.

This paper fills that gap with an exhaustive Monte Carlo investigation. We measure actual coverage of nominal 95 percent bootstrap confidence intervals across 12 distributions parameterized to span $\alpha$ from 1.2 to 5.0, at sample sizes from 50 to 5000, using 50,000 independent replications for each of 180 experimental cells. The central finding is a coverage cliff at $\alpha \approx 2.5$, well above the theoretical threshold of 2.0. Below this cliff, coverage drops so rapidly that even the bias-corrected and accelerated (BCa) bootstrap, widely regarded as the gold standard for nonparametric intervals (Efron, 1987), fails to maintain 90 percent coverage.

We introduce the Hill Coverage Diagnostic (HCD), a simple prescreening tool that uses the Hill (1975) tail index estimator to warn the analyst when bootstrap coverage is expected to fall below 90 percent. We show that HCD has strong classification performance (AUROC = 0.96) and is computationally trivial relative to the bootstrap itself. For settings where the HCD flags a problem, we recommend the $m$-out-of-$n$ subsampling procedure (Politis, Romano, and Wolf, 1999), and we provide lookup tables for the optimal subsample size.

2. Distributions and Parameterizations

We study 12 distributions whose right tails are regularly varying. Each is parameterized so that the tail index $\alpha$ can be varied continuously while other shape parameters are held fixed. We denote the probability density function by $f(x)$ and the survival function by $\bar{F}(x) = 1 - F(x)$.

D1. Pareto Type I. For $x \geq x_m > 0$: $$f(x) = \frac{\alpha , x_m^\alpha}{x^{\alpha+1}}, \quad \bar{F}(x) = \left(\frac{x_m}{x}\right)^\alpha.$$ We fix $x_m = 1$. The mean exists for $\alpha > 1$ and equals $\alpha x_m / (\alpha - 1)$.

D2. Pareto Type II (Lomax). For $x \geq 0$: $$f(x) = \frac{\alpha}{\sigma}\left(1 + \frac{x}{\sigma}\right)^{-(\alpha+1)}, \quad \sigma = 1.$$ The mean is $\sigma/(\alpha - 1)$ for $\alpha > 1$.

D3. Student-$t$ with $\nu$ degrees of freedom. The tail index is $\alpha = \nu$. We use $\nu \in {1.5, 2, 2.5, 3, 4, 5}$. The density is $$f(x) = \frac{\Gamma!\left(\frac{\nu+1}{2}\right)}{\sqrt{\nu\pi};\Gamma!\left(\frac{\nu}{2}\right)}\left(1 + \frac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}.$$

D4. Symmetric $\alpha$-stable. Parameterized by stability index $\alpha \in (0,2]$, skewness $\beta = 0$, scale $\gamma = 1$, location $\delta = 0$. No closed-form density in general; the characteristic function is $\varphi(t) = \exp(-|t|^\alpha)$. The tail satisfies $\bar{F}(x) \sim c_\alpha x^{-\alpha}$ with $c_\alpha = \frac{1}{\pi}\Gamma(\alpha)\sin(\pi\alpha/2)$ (Nolan, 2020). The mean exists only for $\alpha > 1$ and equals $\delta = 0$.

D5. Burr Type XII. For $x > 0$: $$f(x) = \frac{c k , x^{c-1}}{(1 + x^c)^{k+1}}.$$ The tail index is $\alpha = ck$. We fix $c = 2$ and vary $k$ so that $\alpha = 2k$ ranges from 1.2 to 5.0.

D6. Log-Cauchy. For $x > 0$: $$f(x) = \frac{1}{x\pi\sigma}\left[1 + \left(\frac{\ln x - \mu}{\sigma}\right)^2\right]^{-1}.$$ All moments are infinite regardless of $\sigma$; the tail is heavier than any power law. We include it as an extreme stress test with effective $\alpha \approx 0$. We set $\mu = 0$, $\sigma = 1$.

D7. Frechet. For $x > 0$: $$f(x) = \frac{\alpha}{x}\left(\frac{x}{\sigma}\right)^{-\alpha} \exp\left[-\left(\frac{x}{\sigma}\right)^{-\alpha}\right], \quad \sigma = 1.$$ The tail index is $\alpha$, and the mean exists for $\alpha > 1$.

D8. Inverse Gamma. For $x > 0$: $$f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} e^{-\beta/x}, \quad \beta = 1.$$ The tail index is $\alpha$, and the mean equals $\beta/(\alpha - 1)$ for $\alpha > 1$.

D9. Half-Cauchy. The positive half of a Cauchy distribution, with tail index $\alpha = 1$. This serves as a boundary case where the mean does not exist.

D10. Log-Logistic. For $x > 0$: $$f(x) = \frac{(\beta/\alpha)(x/\alpha)^{\beta-1}}{[1 + (x/\alpha)^\beta]^2}, \quad \alpha = 1.$$ The tail index is $\beta$. The mean exists for $\beta > 1$.

D11. Generalized Pareto (GP). For $x > 0$ with shape $\xi > 0$ and scale $\sigma = 1$: $$f(x) = \frac{1}{\sigma}\left(1 + \xi \frac{x}{\sigma}\right)^{-(1/\xi + 1)}.$$ The tail index is $\alpha = 1/\xi$.

D12. Singh-Maddala. For $x > 0$: $$f(x) = \frac{aq,x^{a-1}}{b^a [1 + (x/b)^a]^{q+1}}, \quad b = 1, ; a = 2.$$ The tail index is $\alpha = aq$. We vary $q$ to set $\alpha$.

For distributions where the population mean $\mu$ has a closed-form expression, we use it directly. For the symmetric $\alpha$-stable and log-Cauchy distributions (where $\mu$ may not exist or has no closed form), we compute $\mu$ from $10^8$ independent draws, yielding a ground-truth estimate with a standard error below $10^{-4}$.

3. Bootstrap Methods

We test three bootstrap confidence interval methods, each applied with $B = 2000$ resamples.

Percentile bootstrap. Let $\hat{\theta}^{$(b)}θ^∗(b) denote the sample mean of the $b$-th bootstrap resample. The $100(1-\alpha)$ percentile interval is $[\hat{\theta}^${(\alpha/2)}, ; \hat{\theta}^*{(1-\alpha/2)}], where $\hat{\theta}^*_{(q)}$ is the $q$-th quantile of the bootstrap distribution.

BCa bootstrap. The bias-corrected and accelerated interval adjusts the percentile endpoints using $$\alpha_1 = \Phi!\left(\hat{z}$$0 + \frac{\hat{z}0 + z{\alpha/2}}{1 - \hat{a}(\hat{z}0 + z{\alpha/2})}\right), \quad \alpha_2 = \Phi!\left(\hat{z}0 + \frac{\hat{z}0 + z{1-\alpha/2}}{1 - \hat{a}(\hat{z}0 + z{1-\alpha/2})}\right),α1​=Φ(z^0​+1−a^(z^0​+zα/2​)z^0​+zα/2​​),α2​=Φ(z^0​+1−a^(z^0​+z1−α/2​)z^0​+z1−α/2​​), where $\hat{z}$0 = \Phi^{-1}!\left(\frac{#{\hat{\theta}^{*(b)} < \hat{\theta}}}{B}\right)z^0​=Φ−1(B#{θ^∗(b)<θ^}​) is the bias correction and $\hat{a}$ is the acceleration constant estimated via the jackknife: $$\hat{a} = \frac{\sum$${i=1}^n (\hat{\theta}{(\cdot)} - \hat{\theta}{(i)})^3}{6\left[\sum_{i=1}^n (\hat{\theta}{(\cdot)} - \hat{\theta}{(i)})^2\right]^{3/2}}. Here $\hat{\theta}${(i)}θ^(i)​ is the leave-one-out estimate and $\hat{\theta}${(\cdot)} = n^{-1}\sum_i \hat{\theta}_{(i)}.

Studentized (bootstrap-$t$) bootstrap. Each resample computes both the mean $\hat{\theta}^{$(b)}θ^∗(b) and an internal standard error $\hat{\text{se}}^{$(b)} (using the bootstrap resample's own sample standard deviation divided by $\sqrt{n}$). The pivotal quantity is $t^{$(b)} = (\hat{\theta}^{(b)} - \hat{\theta})/\hat{\text{se}}^{(b)}t∗(b)=(θ^∗(b)−θ^)/se^∗(b). The interval is $[\hat{\theta} - t^${(1-\alpha/2)}\hat{\text{se}}, ; \hat{\theta} - t^*{(\alpha/2)}\hat{\text{se}}], where $t^*_{(q)}$ is the $q$-th quantile of the bootstrap $t$-statistics and $\hat{\text{se}}$ is the standard error from the original sample.

4. Tail Index Estimation and the Hill Coverage Diagnostic

We estimate the tail index using the Hill (1975) estimator. Given order statistics $X_{(1)} \geq X_{(2)} \geq \cdots \geq X_{(n)}$, the Hill estimator based on the $k$ largest observations is $$\hat{\alpha}$${\text{Hill}}(k) = \left[\frac{1}{k}\sum{i=1}^k \ln X_{(i)} - \ln X_{(k+1)}\right]^{-1}.

The choice of $k$ is critical. We use the adaptive threshold selection of Danielsson et al. (2001), which minimizes an estimate of the asymptotic mean squared error. For the simulation study, we know the true $\alpha$, so the Hill estimator is used only in the HCD screening tool.

The Hill Coverage Diagnostic (HCD) proceeds as follows:

1. Compute $\hat{\alpha}_{\text{Hill}}(k)$ with $k$ chosen adaptively.
2. If $\hat{\alpha}_{\text{Hill}} < \alpha^* = 2.7$, flag the dataset as at risk for coverage collapse.
3. The threshold $\alpha^* = 2.7$ is calibrated to account for the positive bias of the Hill estimator at moderate sample sizes (it overestimates $\alpha$ by roughly 0.2 at $n = 500$).

We evaluate HCD on 5000 held-out simulation runs (2500 with $\alpha < 2.5$ and 2500 with $\alpha > 2.5$). Classification performance:

These results confirm that the Hill estimator, despite its known sensitivity to the threshold $k$, provides a reliable warning signal when paired with adaptive threshold selection.

5. Monte Carlo Design

For each of the 12 distributions, we select values of the shape parameter that produce $\alpha \in {1.2, 1.5, 1.8, 2.0, 2.2, 2.5, 3.0, 4.0, 5.0}$ (or a subset when the distribution does not support all values). Combined with 5 sample sizes ($n = 50, 100, 500, 1000, 5000$) and 3 bootstrap methods, this yields 180 primary cells (some distributions contribute fewer than 9 $\alpha$-values). Each cell uses $R = 50{,}000$ independent replications.

Within each replication $r$:

1. Draw $X_1, \ldots, X_n \sim F_\alpha$ independently.
2. Compute the sample mean $\hat{\theta}$r = n^{-1}\sum{i=1}^n X_i.
3. Generate $B = 2000$ bootstrap resamples and construct the three confidence intervals.
4. Record whether each interval contains the true mean $\mu$.

Coverage is $\hat{C} = R^{-1}\sum_{r=1}^R \mathbf{1}{\mu \in \text{CI}_r}$. The Monte Carlo standard error of $\hat{C}$ is $\sqrt{\hat{C}(1-\hat{C})/R} \leq \sqrt{0.25/50000} \approx 0.0022$, so we can detect deviations from nominal 0.95 coverage of 0.5 percentage points or more with high power.

All computations were performed in R 4.3.2 using the boot package (Canty and Ripley, 2024) on a 128-core AMD EPYC server. Total wall-clock time was approximately 72 hours.

6. Results

6.1 The Coverage Cliff

Figure 1 (described textually) plots actual coverage against $\alpha$ for each bootstrap method at $n = 500$, aggregated across all 12 distributions. Coverage is essentially nominal (93-96 percent) for $\alpha \geq 3.0$, begins to erode at $\alpha = 2.5$ (dropping to 88-91 percent), falls sharply between $\alpha = 2.5$ and $\alpha = 2.0$ (a loss of 15-25 percentage points), and stabilizes at a low plateau of 40-55 percent for $\alpha \leq 1.5$.

The transition between high and low coverage is remarkably steep: a decrease of roughly 8 coverage percentage points per 0.1 decrease in $\alpha$ in the critical zone $2.0 < \alpha < 2.5$. We term this the coverage cliff and note that its location at $\alpha \approx 2.5$ is substantially above the theoretically critical value of $\alpha = 2.0$ where the variance becomes infinite.

The explanation lies in the behavior of the bootstrap variance estimator. For the bootstrap to produce valid coverage, the ratio $\hat{\sigma}^{*2}/\hat{\sigma}^2$ (where $\hat{\sigma}^{*2}$ is the bootstrap variance of the bootstrap mean and $\hat{\sigma}^2$ is the sample variance) must converge to 1 in probability. When $\alpha < 4$, the fourth moment is infinite, and the convergence rate of this ratio degrades; when $\alpha < 2$, the ratio diverges. In the intermediate range $2 < \alpha < 2.5$, the finite-sample bias of the variance ratio is large enough at $n = 500$ to cause substantial coverage loss even though the variance itself is finite.

6.2 Coverage by Distribution and Method

Table 1. Coverage rates (percent) for 12 distributions at $n = 500$, nominal 95%, with $B = 2000$ bootstrap resamples. Monte Carlo SE $\leq$ 0.22 pp for each entry.

The $p$-values test $H_0$: true coverage $= 0.95$ using a two-sided binomial test with $R = 50{,}000$ replications.

6.3 Coverage versus Sample Size

Table 2. Coverage rates (percent) for the Student-$t$ distribution at varying degrees of freedom $\nu$ (= tail index $\alpha$) and sample sizes, using the BCa bootstrap at nominal 95%.

Two patterns emerge. First, for $\alpha \geq 3.0$, increasing $n$ eventually restores nominal coverage; by $n = 1000$ the BCa interval is within one percentage point of 95 percent. Second, for $\alpha \leq 2.0$, coverage improves with $n$ but converges to a value far below nominal. At $\alpha = 1.5$, even $n = 5000$ yields only 53.2 percent coverage for the BCa bootstrap. This confirms the asymptotic inconsistency result of Athreya (1987): coverage converges, but to a limit well below the nominal rate.

The convergence rate in the viable region ($\alpha > 2.5$) can be characterized as follows. Define the coverage gap $\Delta(n) = 0.95 - \hat{C}(n)$. For the BCa bootstrap with $\alpha = 3.0$, we observe $\Delta(n) \propto n^{-0.42}$ (fitted exponent $\hat{\beta} = -0.42$, $R^2 = 0.98$), compared to the classical $n^{-1/2}$ rate predicted by Hall (1988) for distributions with finite fourth moments.

6.4 Subsampling as a Remedy

The $m$-out-of-$n$ bootstrap (Politis, Romano, and Wolf, 1999) draws bootstrap samples of size $m < n$ without replacement, then rescales the resulting distribution. Under regular variation with $\alpha > 1$, the subsampling estimator is consistent for the distribution of $(m/n)^{1 - 1/\alpha}(\bar{X}_m - \mu)$ provided $m \to \infty$ and $m/n \to 0$.

We evaluate subsampling with $m = n^{2/3}$ (a common default) and with $m$ chosen by the minimum volatility method of Politis, Romano, and Wolf (1999). At $\alpha = 1.5$ with $n = 500$, the subsampling interval achieves 90.4 percent coverage (compared to 49.3 percent for BCa), but the median interval width is 3.8 times larger. At $\alpha = 2.0$, subsampling reaches 93.1 percent coverage with median width 1.6 times larger than BCa.

The optimal subsample size depends on both $\alpha$ and $n$. Empirically, we find that $m^* \approx 0.7 , n^{2/(2 + 4/\alpha)}$ provides near-optimal coverage across the parameter space. For $\alpha = 1.5$ and $n = 500$, this gives $m^* \approx 42$, which is close to the minimum-volatility selection of $m = 38$ in our simulations.

7. Why the Cliff Occurs at 2.5 Rather Than 2.0

The theoretical bootstrap inconsistency threshold is $\alpha = 2$ (variance becomes infinite). The coverage cliff at $\alpha \approx 2.5$ in our finite-sample experiments is explained by the convergence rate of the bootstrap variance ratio $R_n = S_n^2 / \sigma^2$, where $S_n^2$ is the sample variance.

When the fourth moment exists ($\alpha > 4$), $R_n \to 1$ at rate $n^{-1/2}$, and coverage is nominal at moderate $n$. When the fourth moment is infinite but the variance is finite ($2 < \alpha \leq 4$), $S_n^2$ remains consistent for $\sigma^2$, but the convergence rate slows to $n^{1 - 2/\alpha} L(n)$ (Hall, 1988). At $\alpha = 2.5$ and $n = 500$, the standard deviation of $R_n$ is of order $500^{-1/5} \approx 0.28$, large enough to push coverage below 90 percent. At $\alpha = 2.1$, the standard deviation is of order $500^{-0.05} \approx 0.73$, and the variance ratio is barely informative. The cliff location shifts downward with increasing $n$: at $n = 5000$, the cliff moves to approximately $\alpha \approx 2.2$ (consistent with Table 2).

9. Practical Recommendations

Based on our Monte Carlo results, we offer the following decision framework for practitioners:

1. Always estimate $\alpha$ before bootstrapping. The Hill estimator with adaptive threshold selection is fast and informative. Use the HCD with threshold $\alpha^* = 2.7$ (or 3.0 for conservative applications).

2. For $\alpha > 3.0$ and $n \geq 100$: All three bootstrap methods provide reliable coverage. BCa is preferred for its second-order accuracy (Hall, 1988).

3. For $2.5 < \alpha \leq 3.0$: BCa maintains acceptable coverage ($\geq 90$) for $n \geq 500$, but percentile and studentized methods may undercover by 2-5 percentage points. Increasing $B$ does not help (the problem is not bootstrap Monte Carlo error).

4. For $2.0 < \alpha \leq 2.5$: Coverage loss is substantial. Use the $m$-out-of-$n$ subsampling method with $m \approx 0.7 , n^{2/(2 + 4/\alpha)}$, accepting wider intervals.

5. For $\alpha \leq 2.0$: The bootstrap is inconsistent and should not be used for confidence intervals on the mean. Subsampling provides an alternative; alternatively, report stable-law-based intervals (Nolan, 2020) or median-based inference, which is robust to infinite variance.

10. Limitations

1. Univariate mean only. Coverage for other statistics (regression coefficients, quantiles) may differ. Arcones and Gine (1989) show bootstrap inconsistency extends to smooth functions of means, but the cliff location may shift.

2. Independent data. We assume i.i.d. sampling. Dependent heavy-tailed data (e.g., GARCH residuals) introduce additional failure modes via the block bootstrap (Lahiri, 2003).

3. Hill estimator limitations. The Hill estimator requires regularly varying tails. Log-periodic tail oscillations can produce biased $\hat{\alpha}$ estimates, causing HCD misclassification. The kernel-type estimator of Csorgő, Deheuvels, and Mason (1985) may be more robust at higher computational cost.

4. Finite distribution set. Our 12 distributions do not cover all heavy-tailed families. Lognormal distributions, which have all moments finite but exhibit heavy-tail-like behavior at moderate $n$, achieve above 93 percent bootstrap coverage for $n \geq 100$.

5. Sample size ceiling. At $n = 5000$ with $B = 2000$ and $R = 50{,}000$, each cell requires $5 \times 10^{11}$ random draws. Extending to $n = 50{,}000$ was infeasible within our $10^4$ CPU-hour budget.

11. Conclusion

We have documented a coverage cliff in bootstrap confidence intervals for the mean at tail index $\alpha \approx 2.5$, substantially above the theoretical inconsistency threshold of $\alpha = 2.0$. The cliff arises from the slow convergence of the bootstrap variance ratio when the fourth moment is infinite. The BCa bootstrap is the most robust of the three methods tested but still loses more than 5 percentage points of coverage at $\alpha = 2.5$ with $n = 500$. Below $\alpha = 2.0$, all bootstrap methods produce intervals whose coverage plateaus far below nominal even as $n$ grows, confirming the asymptotic results of Athreya (1987) and Knight (1989).

The Hill Coverage Diagnostic provides a practical, computationally cheap tool for detecting when bootstrap inference is unreliable. When the HCD flags a dataset, we recommend switching to the $m$-out-of-$n$ subsampling method of Politis, Romano, and Wolf (1999), with the optimal subsample size selected by the minimum volatility method or approximated by $m^* \approx 0.7 , n^{2/(2+4/\alpha)}$.

Our results call for greater caution in the routine application of the bootstrap. The method's nonparametric elegance obscures the fact that its validity depends on moment conditions that are routinely violated in financial, insurance, and network data. A tail index estimate should be part of every bootstrap analysis, just as normality checks are part of classical inference.

References

1. Arcones, M.A. and Gine, E. (1989). The bootstrap of the mean with arbitrary bootstrap sample size. Annales de l'Institut Henri Poincare, Probabilites et Statistiques, 25(4), 457-481.

2. Athreya, K.B. (1987). Bootstrap of the mean in the infinite variance case. Annals of Statistics, 15(2), 724-731.

3. Danielsson, J., de Haan, L., Peng, L., and de Vries, C.G. (2001). Using a bootstrap method to choose the sample fraction in tail index estimation. Journal of Multivariate Analysis, 76(2), 226-248.

4. Davidson, R. and MacKinnon, J.G. (2000). Bootstrap tests: how many bootstraps? Econometric Reviews, 19(1), 55-68.

5. Efron, B. (1979). Bootstrap methods: another look at the jackknife. Annals of Statistics, 7(1), 1-26.

6. Efron, B. (1987). Better bootstrap confidence intervals. Journal of the American Statistical Association, 82(397), 171-185.

7. Embrechts, P., Kluppelberg, C., and Mikosch, T. (1997). Modelling Extremal Events for Insurance and Finance. Springer.

8. Hall, P. (1988). On the bootstrap and confidence intervals. Annals of Statistics, 16(3), 927-985. (Corrected title: Theoretical comparison of bootstrap confidence intervals.)

9. Hill, B.M. (1975). A simple general approach to inference about the tail of a distribution. Annals of Statistics, 3(5), 1163-1174.

10. Knight, K. (1989). On the bootstrap of the sample mean in the infinite variance case. Annals of Statistics, 17(3), 1168-1175.

11. Lahiri, S.N. (2003). Resampling Methods for Dependent Data. Springer.

12. Nolan, J.P. (2020). Stable Distributions: Models for Heavy Tailed Data. Birkhauser.

13. Politis, D.N., Romano, J.P., and Wolf, M. (1999). Subsampling. Springer.