Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
45,939 characters · 6 sections · 0 citation commands
From Replications to Revelations: Heteroskedasticity-Robust Inference
A considerable body of literature has proposed and recommended different specifications for heteroskedasticity-robust inference. For instance, based on Monte Carlo evidence, Long and Ervin (2000) strongly recommend for sample sizes below 250 observations HC3 standard errors rather than Stata's robust default: HC1.\footnote{Similar recommendations to adopt more robust inference methods for smaller samples have been made e.g. by MacKinnon and White (1985), who introduced HC1, HC2 and HC3, Chesher and Jewitt (1987), Chesher and Austin (1991), and Cattaneo, Jansson, and Newey (2018).}
Is this recommendation widely followed in empirical practice? No, the opposite holds true. We analyse the code of regression commands found in the Stata scripts of 4420\ reproduction packages from leading economic journals. 40,571 regression commands specify heteroskedasticity-robust standard errors and 98.1% stick with HC1 standard errors. Appendix A provides more details.
Would inference results substantially differ if the recommendation were more widely followed? We proceed with a subset of 608\ regressions from {155\ } reproduction packages that have fewer than 1000 observations and can be reproduced using the current version of the toolbox repbox.\footnote{The repbox tool chain consists of a series of R packages that I am developing to facilitate and semi-automate methodological meta-studies. This paper constitutes my first application of that toolchain. The toolchain automatically reproduces Stata code found in reproduction packages, performing essential tasks such as automatic file path correction. Additionally, it systematically stores information from regression commands and the underlying data sets, enabling the replication and modification of these regressions in both Stata and R. The general endeavor is complex and remains a work in progress. Currently, repbox does not yet robustly work for all reproduction packages that are theoretically reproducible. Thorough development and documentation will take a lot more time. Another part of repbox and area of ongoing work, which is not yet used in this paper, is the automatic mapping of regressions from reproduction packages to the regression tables displayed in the corresponding articles. Thankfully, the Deutsche Forschungsgemeinschaft (DFG) supports future work on repbox as part of the larger SocEnRep project where it will also benefit from input from colleagues from computer sciences and social sciences.} Appendix B provides details on the sample selection.
For Table 1 we explore significance tests for the null hypothesis that a true regression coefficient $\beta_{k}$ is zero. We run each significance test twice: once with HC1 standard errors and once with HC3 standard errors. HC3 standard errors and p-values are always larger than their HC1 counterparts. Table 1 only considers the subsample of the tests for which the HC1 p-value is below 5% and shows for different intervals of that p-value the fraction of tests that are no longer significant at a 5% level if one uses HC3 standard errors instead.
The results show that the switch from HC1 to HC3 has substantial effects on statistical significance. For instance, fewer than 10% of HC1-based tests with p-values between 4% and 5% remain significant when using HC3. This result holds even for regressions with sample sizes ranging from 250 to 1,000 observations.
While these results are insightful, several important questions remain. For instance, to what extent are HC1-based p-values excessively low, and to what extent are HC3-based p-values overly conservative? Beyond HC3 standard errors, alternative methods for improving robust inference have been proposed. These include HC4 standard errors, introduced by Cribari-Neto (2004); HC2 standard errors with alternative calculations of degrees of freedom for the t-test, suggested by Bell and McCaffrey (2002) and Imbens and Kolesar (2016); Hansen's (2024) modified jackknife estimator, which also extend the degrees of freedom adjustments proposed by Bell and McCaffrey (2002); and wild bootstrap inference, as proposed e.g. by Wu (1986) and Roodman et al. (2019).
Do some of these methods perform better than others across typical situations encountered in economic analyses? How heterogeneous is a method's performance across different situations? Do applied researchers encounter scenarios where these established approaches systematically fail, and novel methods offer meaningful improvements?
To address these questions, we employ a large-scale Monte Carlo study. Monte Carlo simulations are a common tool in research on robust inference, but they typically examine only a small set of regression specifications, usually not based on real world data sets.\footnote{A notable exception is Young (2022), who conducts a large-scale Monte Carlo study based on instrumental variable regressions extracted from 30 reproduction packages of economic articles. Already, in his previous article, Young (2019), hand-collected reproducible regressions from a large set of 53 reproduction packages, but did not yet base the Monte Carlo simulations on those regressions. } We perform Monte Carlo studies based on a large set of 608\ OLS regressions originally conducted in 155\ different reproduction packages of published economic articles. Ideally, our approach offers insights that are representative of the situations typically encountered by empirical researchers. Additionally, this broad scope enables us to investigate heterogeneity by examining how the performance of robust inference methods varies across the different original regression specifications.
In a nutshell, for each original regression of the form
we specify a custom data generating process
with the same $n\times K$ matrix of explanatory variables $X$ as in the original regression.
The true coefficients, $\beta$, are set to zero. In line with the original researchers' assumptions, error terms are heteroskedastic. Concretely, we assume $\varepsilon_{i}\sim N(0,\sigma_{i}^{2})$ for each observation $i=1,\ldots,n$. The specification of the standard deviations, $\sigma_{i}$, of the error terms, $\varepsilon_{i}$, is a detailed process performed separately for each original regression. Initially, multiple candidate FGLS specifications, constructed using random forests, are estimated and calibrated. Subsequently, one candidate is selected by comparing the moments of the original OLS residuals with those of the OLS residuals obtained from Monte Carlo simulations of the various candidates. Details are provided in Appendix C.
For each original regression, we draw $M=10,000$ Monte Carlo samples and compute the p-values for a t-test of the null hypothesis $\beta_{k}=0$ for up to 25 coefficients, $\beta_{k}$, per regression. Coefficients of fixed effects dummies are not tested. Each of the 3280\ tested coefficients from the 608\ original regressions constitutes a distinct test situation, indexed by $s$.
For each test situation, we compare different test specifications $\tau\in\{\text{IID, HC1, HC2, \ensuremath{\ldots}}\}$, which vary by the type of standard error and the specification of the degrees of freedom used in the t-distribution. Mathematical background on the different specifications is provided in Section 2.
Let $p_{\tau,s}(m)$ denote the realized p-value for Monte Carlo sample $m=1,\ldots,M$ in specification $\tau$ and test situation $s$. Our analysis focuses on the 5% significance level. The simulated rejection rate is defined as the proportion of Monte Carlo samples for which the p-value is below 0.05:
where $I(\cdot)$ is the indicator function. Since the null hypothesis is true in all test situations, p-values should be uniformly distributed under a correctly specified t-test. Consequently, the ideal value of the rejection rate $\pi_{\tau,s}^{0.05}$ is 0.05.
We measure deviations from this ideal value using the excess and lack of the rejection rate, defined as:
While an excessive rejection rate increases the risk of false discoveries, lack in rejection rates can lead to under-powered significance tests. Excess is generally regarded as more problematic than an equally high lack. However, opinions may differ regarding acceptable levels of excess and the degree of lack one is willing to tolerate for a given reduction in excess.
Figure 1 shows for each specification $\tau$ the average excess and lack and their distribution across all 3280\ test situations. Consistent with conventional wisdom, average excess decreases when moving in order from HC1, HC2, HC4, to HC3 standard errors, while average lack correspondingly increases.
More surprisingly, in our sample of regressions with no more than 1000 observations, both HC1 and HC2 yield on average more excessive rejection rates than inference based on i.i.d. standard errors, which is consistent only under homoskedasticity.\footnote{Thus, adding the robust option to a Stata regress command, such that HC1 standard errors are used, may, in smaller sample sizes, misleadingly suggest that the resulting standard errors and test results are more conservative than without the robust option.} While HC3 is most conservative in the sense of having the largest average lack in rejection rates, it is not the specification with the lowest average excess.
To partially order specifications, we say that one specification outperforms another on average if it has a lower weighted sum of average excess and average lack if the weight on excess is at least as large as the weight on lack.
We study four specifications employing alternative degrees-of-freedom adjustments. HC2-BM is based on the methods of Bell and McCaffrey (2002) and Imbens and Kolesar (2016), JK-H refers to the jackknife estimator proposed by Hansen (2024), and HC1-PL and HC2-PL implement a novel partial leverage-based degrees-of-freedom adjustment introduced in this paper. Each of these methods outperforms, on average, the HC1 to HC4 specifications that rely on the default $n-K$ degrees of freedom in the t-tests. This result highlights the importance of correctly specifying degrees of freedom.
Among all methods considered, JK-H and HC2-PL outperform all others on average, with no clear ranking between them. Both exhibit very small average excess in their rejection rates of 0.06% and 0.1%, respectively. For comparison, a simulation study assuming uniformly distributed p-values yields an average excess in the rejection rates of 0.09%, which is attributable to noise in the Monte Carlo simulations.
Each pane of Figure 1 also plots the rejection rates $\pi_{\tau,s}^{0.05}$ of all 3280\ test situations, arranged in increasing order with the quantile level on the x-axis. Blue and red shaded areas represent rejection rates $\pi_{\tau,s}^{0.05}$ associated with positive lack and excess, respectively. The average lack and excess for specification $\tau$ correspond to the total size of the respective blue and red areas.
There is substantial heterogeneity. No specification $\tau$ exhibits uniformly excessive or uniformly lacking rejection rates across all test situations. As shown on the left-hand side of each pane, all specifications have rejection rates close to zero in some test situations. Conversely, except for JK-H, HC1-PL and HC2-PL, all specifications include a few test situations with extremely high rejection rates exceeding 50%.
Figure 2 plots the rejection rates against the sample size of each test situation. While highly excessive rejection rates are more likely for smaller sample sizes, we find them also in test situations with moderate sample sizes.
The next section provides mathematical background on robust inference and deeper insights into the role of partial leverages and leverages. We will also pin down the subset of test situations that have highly excessive rejection rates for specifications HC2-HC4 and HC2-BM and propose remedies. Section 3 then briefly concludes with additional observations, like the extension of some insights to cluster robust standard errors.
Due to the significant computational demands, wild bootstrap methods are analyzed for only a subset of test situation, with the corresponding results presented in Appendix D. While wild bootstrap specifications outperform the conventional HC1 and HC2 specifications, they are outperformed by the HC2-BM, JK-H, HC1-PL, and HC2-PL specifications that utilize customized degrees of freedom.
This section provides theoretical background on heteroskedasticity-robust inference and yields deeper insights into the role of partial leverages and leverages. Consider the linear regression model:
where $y$ is the dependent variable, $X$ is the matrix of explanatory variables, $\beta$ is the vector of coefficients. The error terms $\varepsilon_{i}$ are independently distributed with mean zero and variance $\sigma_{i}^{2}$.
Using the Frisch-Waugh-Lovell (FWL) theorem, coefficient $\beta_{k}$ can be estimated with the simpler model:
where $\tilde{y}_{k}$ and $\tilde{x}_{k}$ are the residuals from regressing $y$ and $x_{k}$, respectively, on all other explanatory variables in $X$ except for $x_{k}$. Estimation of model ((ref)) yields the same OLS residuals $\hat{\varepsilon}$ and estimator $\hat{\beta}_{k}$ as the original model ((ref)).\footnote{The FWL representation also enhances computational efficiency of our Monte Carlo simulations. While certain elements, in particular, the leverages $h$ (introduced further below), cannot be derived from the FWL representation, these elements need to be computed only once per original regression. Since calculation of a full-blown $K\times K$ variance-covariance matrix can be avoided, the FWL representation considerably accelerates computations that must be repeated for each Monte Carlo sample.}
The partial leverage for observation $i$ with respect to explanatory variable $x_{k}$ is defined as:
where $\tilde{x}_{k,i}$ is the $i$-th element of $\tilde{x}_{k}$. Partial leverages satisfy $0\leq\tilde{h}_{k,i}\leq1$ and $\sum_{i=1}^{n}\tilde{h}_{k,i}=1$.\footnote{Furthermore, for $n\geq2$, we have $\tilde{h}_{k,i}<1$ if the original regression model includes a constant.}
The FWL representation yields a simple formula for the true variance of $\hat{\beta}_{k}$:
The impact of error variance $\sigma_{i}^{2}$ on the variance of $\hat{\beta}_{k}$ is thus proportional to the partial leverage $\tilde{h}_{k,i}$.
Heteroskedasticity-robust variance estimators for $\hat{\beta}_{k}$ of types HC0 to HC4 can all be expressed in the general form:
with
where $\hat{\varepsilon}_{i}$ is the OLS residual for observation $i$, and $\alpha_{i}^{\tau}$ is an adjustment factor that depends on the specification $\tau$.
Intuitively, the term $\hat{\sigma}_{i}^{\tau}$ estimates the unknown standard deviation $\sigma_{i}$ of $\varepsilon_{i}$. Of course, being based on a single residual, $\hat{\sigma}_{i}^{\tau}$ will never consistently estimate $\sigma_{i}$. Yet, an important insight of Eicker (1967), Huber (1967), and White (1980) is that, usually, one can still consistently estimate the (co-)variances of $\hat{\beta}$. Nowadays, their proposed estimator is often referred as HC0: it sets $\alpha_{i}^{HC0}=1$. While consistent, HC0 can exhibit severe bias in small samples and is rarely used in practice. MacKinnon and White (1985) introduced three variants, HC1, HC2 and HC3, with improved small-sample properties.
The HC1 adjustment factor is:
This correction mirrors the degrees of freedom adjustment used in the unbiased estimator $\hat{\sigma}^{2}=\frac{1}{n-K}\sum_{i=1}^{n}\hat{\varepsilon}_{i}^{2}$ for the error variance $\sigma^{2}=E[\varepsilon_{i}^{2}]$ under homoskedasticity.
The HC2 adjustment factor is given by:
where $h_{i}$ is the leverage of observation $i$, defined as the $i$-th diagonal element of the hat matrix:
Leverages satisfy $0\leq h_{i}\leq1$ and $\sum_{i=1}^{n}h_{i}=K$. HC2 is motivated by the property that, under homoskedasticity, $E[\hat{\varepsilon}_{i}^{2}]=(1-h_{i})\sigma^{2}$. Furthermore, given certain regularity conditions Pustejovsky and Tipton (2017) show that HC2 constitutes an unbiased variance estimator.
The HC3 adjustment factor is also based on leverages:
As shown in Hansen (2022), HC3-adjusted residuals $\sqrt{\alpha_{i}^{HC3}}\hat{\varepsilon}_{i}$ are equivalent to the leave-one-out prediction error $y_{i}-x_{i}\hat{\beta}^{(i)}$, where $\hat{\beta}^{(i)}$ denotes the OLS estimator obtained from the regression excluding observation $i$. The HC3 estimator can also be interpreted as a jackknife estimator, satisfying\footnote{See MacKinnon et al. (2023) for a derivation. MacKinnon and White (1985) introduced the jackknife formulation \[ \hat{V}_{k}^{JK}=\frac{N-1}{N}\sum_{i=1}^{N}(\hat{\beta}_{k}^{(i)}-\bar{\beta}_{k}^{JK})^{2} \] as the HC3 variance estimator where $\bar{\beta}_{k}^{JK}$ denotes the mean of the leave-one-out estimators $\hat{\beta}\halfnote_{k}^{(i)}$. The current HC3 formulation was popularized by Davidson and MacKinnon (1993).}
Hansen (2024) adapts the jackknife estimator ((ref)) by employing generalized Moore-Penrose inverses in cases an $\hat{\beta}^{(i)}$ cannot otherwise be computed. The HC4 estimator, introduced by Cribari-Neto (2004), aims to better handle cases with high leverages by modifying the adjustment factor to:
where $\delta_{i}=\min\left\{ 4,\frac{nh_{i}}{K}\right\} $.
In our specifications HC1 to HC4 we follow Stata and other software packages and use $n-K$ degrees of freedom in the $t$-tests. Yet, it is evident from ((ref)) that if partial leverages are concentrated in a few observations, the standard errors are typically imprecisely estimated as they depend strongly on the residuals of those few observations. Choosing simply $n-K$ degrees of freedom in the $t$-test does not account for such imprecisions arising from highly concentrated partial leverages.
The idea behind our newly proposed methods HC1-PL and HC2-PL is to base the degree of freedoms in the $t$-test on the concentration of partial leverages. A widely used concentration measure in competition policy is the Herfindahl-Hirschman index, defined as the sum of the squared market shares of all competitors in a market. Similar to market shares, partial leverages are non-negative and add up to 1. Thus the concentration of partial leverages can be measured by the sum of squared partial leverages. We denote the inverse of this measure as:
We refer to $\tilde{n}_{k}$ as the partial-leverage-adjusted sample size. It satisfies $1\leq\tilde{n}_{k}\leq n$. If all observations have the same partial leverage then $\tilde{n}_{k}=n$, while $\tilde{n}_{k}=1$ in the limit case that the partial leverage becomes concentrated in a single observation.
For HC1-PL and HC2-PL we specify the degrees of freedom in the $t$-test as $\tilde{n}_{k}-1$ and use HC1 and HC2, respectively, as standard errors. Appendix E provides a justification for this partial-leverage-based degrees of freedom by deriving it from a Satterthwaite approximation, similar to the approach of Bell and McCaffrey (2002).
To understand why we specify the degrees of freedom as $\tilde{n}_{k}-1$ instead of $\tilde{n}_{k}$, consider the following. As long as the original sample contains at least two observations, we have $\tilde{n}_{k}>1$. Since a $t$-distribution can be defined for any fractional degrees of freedom strictly greater than zero, this adjustment ensures properly defined degrees of freedom whenever $n\geq2$. In the degenerate case of a single observation, no standard error can be computed, and the corresponding $t$-distribution with zero degrees of freedom becomes degenerate. Using $\tilde{n}_{k}-1$ as the degrees of freedom allows for continuous convergence to this degenerate case as $\tilde{n}_{k}\to1$.
Figure (ref) provides additional the empirical support for our proposal. For the standard HC1-HC4 specifications, as well as for specification HC2-BM: all test situations with highly excessive rejection rates exhibit very small partial leverage adjusted sample sizes. By imposing very wide t-distributions in those cases, HC1-PL and HC2-PL effectively mitigate those excessive rejection rates. Hansen's (2004) jackknife estimator also robustly handles cases with small partial leverage.
The insights discussed thus far do not yet provide a complete picture. We now demonstrate that excessive rejection rates are particularly pronounced in a subset of regressions involving observations with a leverage $h_{i}$ of one and high partial leverages. We refer to a leverage value of one as full leverage. For cases with full leverage, we propose a modification to the currently implemented methods for computing standard errors.
Indeed, 15.6% of the regressions in our replication sample include at least one observation with full leverage. Full leverage occurs whenever one or several explanatory variables can perfectly predict the corresponding observation $i$ such that $\hat{\varepsilon}_{i}=0$. In that case $\hat{\varepsilon}_{i}$ provides no information about $\sigma_{i}$. Our definitions above for the HC2, HC3, and HC4 variance estimators are then not well-defined as $\hat{\sigma}_{i}^{\tau}$ involves the mathematically indeterminate expression $0/0$.
Fortunately, we don't need to predict $\sigma_{i}$ if $x_{k}$ is not required to perfectly predict observation $i$, i.e. if a modified regression that omits explanatory variable $x_{k}$ also has full leverage for observation $i$. We then have $\tilde{x}_{k,i}=0$ and also a zero partial leverage $\tilde{h}_{k,i}=0$. This means $\sigma_{i}$ does not affect the true variance of $\hat{\beta}_{k}$ and the general variance estimator $\hat{V}_{k}^{\tau}$ specified in ((ref)) does not depend on $\hat{\sigma}_{i}^{\tau}$.
In contrast, if observation $i$ has full leverage and a positive partial leverage $\tilde{h}_{k,i}>0$, explanatory variable $x_{k}$ is essential to perfectly predict observation $i$: $\sigma_{i}$ then influences the actual variance of $\hat{\beta}_{k}$ and variance estimators $\hat{V}_{k}^{\tau}$ depend on which value is assumed for $\hat{\sigma}_{i}^{\tau}$. The larger the partial leverage $\tilde{h}_{k,i}$, the more relevant becomes the assumed value of $\hat{\sigma}_{i}^{\tau}$ for the variance estimator.
There is not yet a consistent treatment of cases with full leverage across statistical software packages. The sandwich package in R (Zeileis, 2004) returns an NaN value for HC2 to HC4 standard errors if an observation has full leverage.
A similar non-invertibility problem arises in the context of cluster-robust standard errors. Pustejovsky and Tipton (2017) theoretically justify using a generalized Moore-Penrose inverse, which Kolesar (2023) and the corresponding R package dfadjust adapt to HC2 heteroskedasticity-robust standard errors. The Moore-Penrose inverse of $1-h_{i}$ for $h_{i}=1$ is zero, i.e. one sets $\hat{\sigma}_{i}^{\tau}=0$ if $i$ has full leverage. Poetscher and Preinerstorfer (2023) adopt this approach for HC3 and HC4 as well.
Stata's approach is not yet well-documented. But experiments and personal communications suggest that for observations with full leverage, Stata evaluates $\alpha_{i}^{\tau}\hat{\varepsilon}_{i}^{2}=0/0$ numerically for HC3, which yields results determined by rounding errors. For HC2 computation, as well as HC2-BM computation using Stata 18's new option vce(hc2, dfadjust), Stata seems to set $\hat{\sigma}_{i}^{\tau}=0$ for observations with full leverage.
Being in line with existing conventions, we have set $\hat{\sigma}_{i}^{\tau}=0$ for observations with full leverage in our Monte Carlo simulations reported in Section 1. It is important to note that the theoretical justification, that under this approach HC2 is then an unbiased estimator of the true variance, only applies for the case that every observation with full leverage has zero partial leverage. Yet, in those cases also any other finite value for $\hat{\sigma}_{i}$ yields the same variance estimator. For cases where the choice of $\hat{\sigma}_{i}$ matters, setting $\hat{\sigma}_{i}=0$ actually seems a recipe for underestimating standard errors.
Table (ref) shows the performance of different inference specifications for different subset of test situations. Columns (1) and (2) show results for tests without full leverage: HC2-BM and HC2-PL perform both very well. Columns (3) and (4) study regressions that have observations with full leverage, but all of them have zero partial leverage. Those cases constitute 95.6% of all test situations with full leverage. Again both HC2-BM and HC2-PL perform both very well. Columns (5) and (6) show results for the very small number of 31 test situations where observations with full leverage are accompanied by positive partial leverages. In those cases, all specifications except for JK-H, HC1-PL and HC2-PL exhibit extremely excessive rejection rates, with average excess well above 20%.
We suggest selecting a different value for $\hat{\sigma}_{i}^{\tau}$ than 0 for cases with full leverage: one natural candidate is the homoskedastic standard error estimate $\text{\ensuremath{\hat{\sigma}}}$. As shown in Columns (7) and (8), this adjustment reduces average excess of rejection rates to more reasonable levels. Of course, we cannot rule out the possibility that, in some real-world data-generating processes, observations with full leverage systematically exhibit higher standard errors $\sigma_{i}$ than the homoskedastic estimate $\hat{\sigma}$. Nonetheless, setting $\hat{\sigma}_{i}^{\tau}=\hat{\sigma}$ always yields more conservative results compared to setting $\hat{\sigma}_{i}^{\tau}=0$.
Interestingly, these 31 particular test situations are more diverse than one might suspect. For example, only in 19 cases $x_{k}$ is the obvious dummy variable that is equal to $1$ only for the observation with full leverage and thus perfectly predicts it. If the regressions had no control variables that are correlated with $x_{k}$ such a dummy variable would exhibit a very high partial leverage for the observation $i$ that is perfectly predicted by $x_{k}$. Yet, with control variables, the structure of the FWL residual $\tilde{x}_{k}$ may substantially differ from $x_{k}$ and partial leverages may be far less concentrated. In our sample one such dummy variable even has a partial leverage of only 0.6% at the perfectly predicted observation and and a large partial-leverage adjusted sample size of $\tilde{n}_{k}=93.5$. For these test situations, empirical rejection rates remain close to 5% for both HC2-PL and HC2-BM, regardless of whether $\hat{\sigma}_{i}^{\tau}$ is set to $0$ or $\hat{\sigma}$.
However, most of the 31 test situations exhibit highly concentrated partial leverages. Specifically, for 20 test situations, we observe $\tilde{n}_{k}<2$. This explains the very conservative rejection rates of HC1-PL and HC2-PL reported in columns (5) to (8). Those conservative rejection rates reflect high uncertainty about the actual value of the standard error. Technically, by fixing for observations with full leverage $\hat{\sigma}_{i}^{\tau}$ to $0$ or $\hat{\sigma}$, the actual variance of the standard error can be much smaller than if $\hat{\sigma}_{i}^{\tau}$ is determined by the noisy OLS residual $\hat{\varepsilon}_{i}$. It is nonetheless prudent to treat these cases as if we are highly uncertain as we don't know how well the fixed $\hat{\sigma}_{i}^{\tau}$ approximates the unknown true $\sigma_{i}$. Also Hansen's jacknife estimator JK-H has very conservative rejection rates in those test situations.
We conclude briefly with our key recommendations and additional remarks.
Overall, we recommend using Hansen's (2024) jackknife estimator or JK-H, or the novel HC2-PL specification for heteroskedasticity-robust inference, with HC2-BM being in most cases an equally good alternative. A limitation of the HC2-PL specification, unlike JK-H and HC2-BM, is that it cannot be readily extended to hypothesis tests involving multiple coefficients. Computing the jackknife estimator JK-H likely requires a longer runtime compared to HC2-PL and HC2-BM; however, for a single regression, this is typically not a concern on modern computers.
We further recommend the following modification for all specifications but JK-H. For observations $i$ that have full leverage, set $\hat{\sigma}_{i}^{\tau}$ to a more conservative value than 0, such as the homoskedastic error estimate $\hat{\sigma}$. Additionally, software packages could report the percentage of the standard error computed using the homoskedastic error estimate by summing the partial leverages of all observations $i$ with full leverage.
MacKinnon et al. (2023b) propose reporting various diagnostic statistics to evaluate the robustness of inference, including the distribution of partial leverages. As compact concentrations measures, one could report the partial leverage-adjusted sample sizes, $\tilde{n}_{k}$, potentially accompanied by a warning for very low values.
All methods can be easily adapted to compute confidence intervals. Moreover, for methods that compute customized degrees of freedom $\nu_{k}^{\tau}$, one could follow the suggestion of Imbens and Kolesar (2016), and report adjusted standard errors by multiplying the original standard errors with $q_{\nu_{k}^{\tau}}(0.975)/q_{n-K}(0.975)$ where $q_{\nu}$ is the quantile function of the t-distribution with $\nu$ degrees of freedom.
In our Monte Carlo studies, HC1-PL was relatively close behind the performance of HC2-PL and JK-H. The reghdfe Stata command (Correia, 2017) and the fixest R package (Berge, 2018), which are widely used for fixed effects regressions, provide heteroskedasticity-robust standard errors based solely on HC1. This is likely because extending the performance gains from fixed-effects absorption to the computation of hat values, which are required for HC2-, HC3-, or HC4-based standard errors, or to JK-H estimation, is non-trivial.\footnote{Personal communication and the inclusion of HC2 standard errors with absorbed fixed effects in the areg function of Stata 18 suggest that StataCorp might have developed a yet unpublished, performant method for this computation.} In contrast, absorbed fixed effects do not pose any challenges for computing partial leverages and HC1 standard errors. Therefore, HC1-PL might be a promising alternative specification for regressions with absorbed fixed effects.
This paper does not study robust inference methods proposed by Cattaneo et al. (2018) and by Poetscher & Preinerstorfer (2023). I still need to better understand those methods and also need to figure out whether a fast implementation for our large scale Monte Carlo study is possible. If not included in a future version of this paper, hopefully once the repbox toolbox is fully developed and well documented, it will facilitate such studies by more skilled researchers.
Variants of JK-H and HC2-BM have also been developed for cluster-robust inference. Also HC1-PL and HC2-PL can be readily extended to cluster-robust inference. The partial leverage of cluster $g$ is defined as the sum of the partial leverages of all observations within cluster $g$. The inverse Herfindahl-Hirschman index is then computed by treating each cluster as a single observation, yielding a partial-leverage-adjusted number of clusters, $\tilde{G}_{k}$. A comprehensive Monte Carlo assessment for cluster-robust inference is planned for a separate paper.