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.
105,481 characters · 22 sections · 0 citation commands
A More Robust t-Test
\setcounter{page}{1}
The usual t-test for inference about the mean of a population from an i.i.d. sample is a key building block of statistics and econometrics. Not only does it have many direct applications, but also many other standard forms of inference reduce to the application of a t-test applied to a suitably defined population. For example, consider a linear regression with scalar regressor, $Y_{i}=X_{i}\beta +\varepsilon _{i}$, $\mathbb{E} [X_{i}\varepsilon _{i}]=0$. A test of the null hypothesis $H_{0}:\beta =\beta _{0}$ reduces to a test of $\mathbb{E}[W_{i}]=0$ for $ W_{i}=(Y_{i}-X_{i}\beta _{0})X_{i}$, and the usual t-statistic computed from the i.i.d. sample $W_{i}$ amounts to a specific version of the usual heteroskedasticity robust test suggested by \citeasnoun{White80}. Under clustering that allows for arbitrary correlations between $\varepsilon _{j}$ for all $ j\in \mathcal{C}_{i}$, $i=1,\ldots ,n$, the effective observations become $ W_{i}=\sum_{j\in \mathcal{C}_{i}}(Y_{j}-X_{j}\beta _{0})X_{j}$. In the presence of additional controls $Y_{i}=X_{i}\beta +Z_{i}^{\prime }\gamma +\varepsilon _{i}$, the equivalence to the inference for the mean problem holds approximately after projecting $Y_{i}$ and $X_{i}$ off $Z_{i}$. This further extends to instrumental variable regression and parameters estimated by GMM.
The asymptotic validity of standard t-statistic based inference relies on two arguments. First, the law of large numbers implies that the variance estimator in the denominator has negligible estimation error. Second, a central limit theorem applied to the numerator yields approximate normality. Underlying populations with heavy tails are a threat to both. Even if the second moment exists, so that t-statistic based inference is asymptotically justified, large samples might be required before these approximations become accurate.
The effective sample size in empirical work is often considerably smaller than the raw number of observations, and not all that large. This can arise because researchers are interested in inference for smaller subgroups, or because nonparametric kernel estimators are employed that effectively depend only on relatively few observations, or because the relevant variation only stems from a small fraction of observations, such as in studies about rare events. It is also very common for standard errors to be clustered, reducing the effective number of independent observations to the number of clusters, which tends to be only moderately large. What is more, in many applications clusters are of fairly heterogeneous size (think of the 50 states of the U.S., say). Even if none of the variables under study are heavy-tailed, a substantial portion of the parameter sampling variation will then stem from the randomness of the large clusters, inducing a form of heavy-tailedness in the resulting $W_{i}$ variables. See Section (ref) for an illustration.
This paper develops an alternative to the t-test that performs more reliably when the underlying population has potentially heavy tails. The focus is exclusively on the case of moderately heavy tails, that is, the first two moments of $W_{i}$ exist, so that asymptotically, standard t-statistic based inference is valid. The aim is to devise an inference method that does not overreject if the underlying population has moderately heavy tails, without losing much in terms of efficiency if the underlying population has light tails. The theoretical development only concerns the canonical inference for the mean problem, but we show how to adapt the procedure to also obtain more reliable inference about scalar parameters estimated by GMM, including the case with clustering. In such more general contexts, the quality of standard inference is poor if the induced $W_{i}$ has heavy tails, for which it is neither necessary nor sufficient that the variables under study are individually heavy-tailed.
To describe the key idea, consider the hypothesis test of $H_{0}:\mathbb{E} [W_{i}]=0$ against $H_{a}:\mathbb{E}[W_{i}]\neq 0$ based on an i.i.d. sample $W_{i}$, $i=1,\ldots ,n$, from a population $W$ with cumulative distribution function $F$. For expositional ease, suppose that $F$ has a thin left tail, but a potentially heavy right tail. For some given $k$, let $\mathbf{W} ^{R}=(W_{1}^{R},W_{2}^{R},\cdots ,W_{k}^{R})^{\prime }$ be the $k$ largest order statistics, with $W_{1}^{R}$ the sample maximum. Conditional on $ \mathbf{W}^{R}$, the remaining \textquotedblleft small\textquotedblright\ observations $W_{i}^{s}$, $i=1,\ldots ,n-k$ are i.i.d. draws from the truncated distribution with c.d.f. $F(w)/F(W_{k}^{R})$ for $w\leq W_{k}^{R}$ . The mean of this truncated distribution under $H_{0}$ is no longer zero, however, but is given by $-m(W_{k}^{R})<0$, where
Note that $m(w)$ for $w$ large is determined by the properties of $F$ in its right tail.
The idea now is to apply three asymptotic approximations. First, invoke standard extreme value theory to obtain an approximation for the distribution of $\mathbf{W}^{R}$ in terms of a (joint)\ extreme value distribution governed by three parameters describing location, scale and shape. Second, apply the central limit theorem to the conditional i.i.d. sample of remaining observations $W_{i}^{s}$ from the truncated (and hence no longer heavy-tailed) distribution to argue that $ (n-k)^{-1}\sum_{i=1}^{n-k}W_{i}^{s}$ is approximately normal with mean $ -m(W_{k}^{R})$ under $H_{0}$ (and arbitrarily different mean under the alternative $H_{a}$). Third, by the same arguments that justify extreme value theory, obtain an approximation to $m(w)$ in terms of the three parameters that govern the distributional approximation of $\mathbf{W}^{R}$.
These approximations lead to a parametric approximate joint model of $ k+1$ statistics: $\mathbf{W}^{R}$ is jointly extreme value, and $ (n-k)^{-1}\sum_{i=1}^{n-k}W_{i}^{s}$ is normally distributed with a mean that, under $H_{0}$, depends on $W_{k}^{R}$, and the parameters of the extreme value distribution. For given $k$, this is a small sample nonstandard parametric testing problem, and one can construct tests that are exactly valid under the approximate parametric model. Specifically, we apply computational techniques similar to those developed in \citeasnoun{Elliott15} to determine a powerful that is of level $\alpha $ in this parametric model. Once the test is applied to the original mean testing problem, it is no longer of level $\alpha $ by construction. But the explicit modelling of the potentially moderately heavy tail via extreme value theory might improve performance over the usual t-test.
The main theoretical result of this paper corroborates this conjecture by considering higher order improvements for populations with finite variance, but that do not possess a third moment, and for which extreme value theory applies. We consider asymptotics in which $k$ is a fixed number that does not vary as a function of $n$. In this way, the asymptotics reflect that moderately large samples only contain limited information about the tail properties of the underlying population. We show that the approximation error of the parametric model for $k$ fixed induces an error in the rejection probability in the mean testing problem that converges to zero faster than the error in rejection probability of the usual t-test. In that sense, the new approach yields a refinement over the usual t-test and provides theoretical support for the usefulness of the new perspective.
A natural alternative to obtain more accurate approximations is to consider the bootstrap. \citeasnoun{Bloznelis03} show that the percentile-t bootstrap provides a refinement whenever the underlying population has at least three moments. A second, apparently new theoretical result shows that the bootstrap does not provide a refinement when the underlying population has between two and three moments.
The new method readily generalizes to the case where both tails are potentially heavy. The approximate parametric model then consists of $2k+1$ statistics, with $k$ joint extreme value observations from the left tail governed by three parameters, $k$ extreme value observations from the right tail governed by their own three parameters, and the conditionally normal average of the middle observations. Since in most applications, there are no compelling reasons to assume any constraints between the properties of the left and right tail, the approximate parametric problem is thus indexed by a six dimensional nuisance parameter. We use a version of the the algorithm of \citeasnoun{Elliott15} to numerically determine a powerful test in this parametric problem for selected values of $k$. The large nuisance parameter space turns this into a major computational challenge. Once a powerful valid test has been determined, however, applying it in practice is entirely straightforward and does not pose any significant computational burden. This includes its use to obtain more reliable inference about scalar parameters estimated by GMM with potentially clustered errors; see Section (ref) for details.
Our preferred default method uses $k=8$ and is appropriate when the sample consists of at least 50 independent clusters or observations.\footnote{ We also provide an alternative, even more robust test for $k=4$ that is applicable to samples with as few as $25$ indepedent clusters or observation. } The tests were determined for various significance levels, enabling the construction of confidence intervals at the 90%, 95% and 99% level via test inversion, and the computation of p-values. Corresponding tables and STATA code is provided in the replication files.
Monte Carlo simulations show that the new approach leads to much better size control in moderately large samples compared to existing methods, at fairly small cost in terms of average confidence interval length for thin-tailed populations. This is true in the canonical inference for the mean case, as predicted by the theory, but also when comparing two means, and for inference about regression coefficients under clustering. In one design, the clusters are Metropolitan Statistical Areas, which are fairly heterogeneous in size. As discussed above, this heterogeneity induces the resulting $W_{i}$ to be quite heavy-tailed, which leads to poor performance of standard cluster robust inference. A moderately large number of heterogenous clusters (say, no more than 100 or 200) is quite common in empirical work, making the new approach a potentially attractive alternative in such settings.
The remainder of the paper is organized as follows. The next section discusses the related literature and provides a review of known results from extreme value theory and the approximate distribution of t-statistics. Section 3 contains the new theoretical results in the inference for the mean problem. Section 4 provides details on the construction of the new test. Small sample simulation results are reported in Section 5. Section 6 concludes.
The new method \textquotedblleft robustifies\textquotedblright\ the usual t-test in the sense of providing more reliable inference under moderately heavy tails. The classic robustness literature (see, for instance, \citeasnoun {Huber96} for an overview) is based on a very different notion; in this literature, it is assumed that the observations are contaminated, with a small fraction not stemming from the population of interest. This in turn raises the question what type of estimands can still be reliably learned about, and how to do so efficiently. In contrast, the estimands considered here are defined relative to the distribution that generated the data. Which of these views is appropriate depends on the application, and in particular whether relatively extreme observations are part of the population of interest, or rather are induced by measurement errors. It is also possible to combine the approaches, such as applying the new test in a regression of winsorized variables (with the estimand then defined relative to a winsorized population).
Even under the assumption that the data is entirely uncontaminated, informative inference requires assumptions beyond the existence of moments: The classic impossibility result of \citeasnoun{Bahadur56} shows that one cannot learn about the population mean from i.i.d. samples of any size, even if all moments are assumed to exist. The substantial assumption pursued here is that the population tails are such that extreme value theory provides reasonable approximations. This effectively amounts to an assumption that the tails of the underlying distribution are approximately (generalized) Pareto. Given the theoretical prevalence and empirical success of extreme value theory for learning about the tail of distributions (for overviews and references, see, for instance, \citeasnoun{Embrechts97} or \citeasnoun{DeHaan07}), this seems a reasonably general starting point, especially given that some assumption must be made. What is more, the approximate Pareto tail is only imposed in the extreme tail with approximate mass of $k/n$ for $k$ fixed, which is enough to ensure that the largest (and smallest) $k$ observations are governed by extreme value theory.
\citeasnoun{Muller17} pursue this \textquotedblleft fixed-$k$\textquotedblright\ approach for the purpose of inference about tail properties, such as extreme quantiles. In contrast, the remaining literature on the modelling of tails considers asymptotics where $k=k_{n}$ diverges with the sample size. In large samples, $k_{n}$ diverging asymptotics allow for consistent estimation of tail properties, at least pointwise for a fixed population. In practice, though, the approximations generated from $k_{n}$ diverging asymptotics are not very useful for, say, samples of size $n=50$ or $n=100$, as there are only a handful of observations that can usefully be thought of as stemming from the tail, so that any approximation that invokes \textquotedblleft consistency\textquotedblright\ of tail property estimators becomes misleading.
The separate analysis of the largest and remaining terms of a sum of independent random variables goes back to at least \citeasnoun{Csorgo88}; also see \citeasnoun{Zaliapin05}, \citeasnoun{Kratz2014} and \citeasnoun{Mueller19}. The relatively closest precursors to this work are Peng (2001, 2004)\nocite{Peng01}\nocite {Peng04}\ and \citeasnoun{Johansson03}. These authors are concerned with inference about the mean from an i.i.d. sample under very heavy tails, that is, the underlying population has less than two moments. For such populations, the usual t-statistic does not converge to a normal distribution. Peng (2001, 2004) and \citeasnoun{Johansson03} suggest estimating the contribution of the two tails to the overall mean by consistently estimating the tail Pareto parameters using the smallest and largest $k_{n}$ observations, with $k_{n}$ diverging, and combining those estimates with the estimate of the mean of the remaining middle observations.
Another approach to overcome the \citeasnoun{Bahadur56} impossibility result is to assume bounded support, with known bounds. \citeasnoun{Romano00}, \citeasnoun{Schlag07} and \citeasnoun{Gossner13} derive corresponding methods.
Let $W_{1}^{R}\geq W_{2}^{R}\geq ...\geq W_{k}^{R}$ denote the largest $k$ order statistics from an i.i.d. sample from a population with distribution $ F $. Suppose the right tail of $F$ is approximately Pareto in the sense that for some scale parameter $\sigma >0$ and tail index $\xi >0$
so that the second moment of $W$ exists if and only if $\xi <1/2$. Then $W$ is in the maximum domain of attraction of the Fr\'{e}chet limit law
where $X_{1}^{-1/\xi }\sim E_{1}$ with $E_{1}$ an exponentially distributed random variable.
As is well known (see, for instance, Theorem 2.8.2 of \citeasnoun{Galambos78}), ( (ref)) implies that extreme value theory also holds jointly for the first $k\,$order statistics
The distribution of $\mathbf{X}$ satisfies $\{X_{j}^{-1/\xi }\}_{j=1}^{k}\sim \{\sum_{l=1}^{j}E_{l}\}_{j=1}^{k}$, where $E_{l}$ are i.i.d. exponential random variables.
Since the new theoretical results of this paper concern rates of convergence, a suitable strengthening of the approximate Pareto tail assumption ((ref))\ is needed. \citeasnoun{Falk04} define the $\delta $ -neighborhood of the Pareto distribution with index $\xi $ as follows.
Theorem 5.5.5 of \citeasnoun{Reiss89} shows that under Condition 1, ((ref) ) provides accurate approximations in the sense that
for $\delta \leq 1$, where the supremum is taken over all Borel sets $ B\subseteq \mathbb{R} ^{k}$.
Many heavy-tailed distributions satisfy Condition 1: for the right tail of a student-t distribution with $\nu $ degrees of freedom, $\xi =1/\nu $ and $ \delta =2\xi $, for the tail of a Fr\'{e}chet or generalized extreme value distribution with parameter $\alpha $, $\xi =1/\alpha $ and $\delta =1,$ and for an exact Pareto tail, $\delta $ may be chosen arbitrarily large. But there also exist heavy-tailed distributions in the domain of attraction of a Fr\'{e}chet limit law that do not satisfy Condition 1, such as density of the form ((ref)) with $h(x)=1/\log (1+x)$, for example. Under some additional regularity conditions, Theorem 3.2 of \citeasnoun{Falk93} shows Condition 1 to be necessary to obtain an error rate of extreme value approximations of order $n^{-\delta }$ for $\delta >0$. Roughly speaking, Condition 1 thus formalizes the assumption that extreme value theory provides accurate approximations.
Let $T_{n}$ be the usual t-statistic computed from an i.i.d. sample $ W_{1},\ldots ,W_{n}$, where $W_{i}\sim W$. If $\mathbb{E}[W]=0$ and $\mathbb{ E}[W^{2}]<\infty $, then $T_{n}\Rightarrow \mathcal{N}(0,1).$ A seminal paper by \citeasnoun{Bentkus93} establishes a bound on the rate of this convergence which does not require the third moment of $W$ to exist. In particular, \citeasnoun{Bentkus93} show that for some $C>0$ that does not depend on $F$, and $\mathbb{E}[W^{2}]=1$,
where $\Phi (t)=\mathbb{P}(Z<t)$, $Z\sim \mathcal{N}(0,1)$. (The explicit claim of uniformity with respect to $F$ is made in \citeasnoun{Bentkus96}.) This result is key for the new theoretical results of this paper.
To put ((ref)) into perspective, recall that the classic Berry-Esseen bound shows a $n^{-1/2}$ rate of the approximation $ Z_{n}=n^{-1/2}\sum_{i=1}^{n}W_{i}/\sqrt{\mathbb{E}[W^{2}]}\Rightarrow \mathcal{N}(0,1)$, with a constant that involves the third moment of $W$. If $\mathbb{E}[|W|^{3}]$ does not exist, then the classic Berry-Esseen bound is inapplicable. Theorem 5 of \citeasnoun{Petrov75} provides the analogue of ((ref)) with $T_{n}$ replaced by $Z_{n}$. The contribution of \citeasnoun {Bentkus93} is thus to show that the random norming inherent in the t-statistic does not affect Petrov's (1975) result.
Subsequent research by \citeasnoun{Hall04} provides a sharp bound on the rate of convergence: If $\mathbb{E}[W^{2}]<\infty $, their results imply that
with $A_{n}=\mathbf{1}[|W|\leq n^{1/2}]$ is bounded away from zero and infinity uniformly in $n$.
A final relevant result from the literature concerns the bootstrap approximation to the distribution of the t-statistic. Let $\mathbf{W} =(W_{1},\ldots ,W_{n})$, and let $T_{n}^{\ast }$ be a bootstrap draw of $ T_{n}$ from the demeaned empirical distribution of $W_{i}$, conditional on $ \mathbf{W}$. \citeasnoun{Bloznelis03} show that if $F$ is non-lattice and $\mathbb{ E}[|W|^{3}]<\infty $, then
while, for $\mathbb{E}[W^{3}]\neq 0$, $\lim \inf_{n\rightarrow \infty }n^{1/2}\sup_{t}|\mathbb{P}(T_{n}<t)-\Phi (t)|>0$. In other words, as long as $W$ has finite non-zero third moment, the error in the bootstrap approximation to the distribution of the t-statistic is of smaller order than the normal approximation, and the bootstrap provides a refinement over the usual t-test.
To ease exposition, we focus in this section on the case where the left tail of $W$ is light, as in the introduction. The analogous results also hold when both tails are moderately heavy with tail index smaller than $1/2$; we provide an analogue of Theorem (ref) in Appendix (ref) .
Since $3(1/2-\xi )>1/(2\xi )-1$ for $1/3<\xi <1/2$, the triangle inequality implies that $\sup_{t}|\mathbb{P}(T_{n}^{\ast }<t|\mathbf{W})-\mathbb{P} (T_{n}<t)|=O_{p}(n^{1-1/(2\xi )})$, so Theorem (ref) shows that the bootstrap does not provide a refinement if the underlying population has between two and three moments, at least as long as the population has an approximate Pareto tail. This result is apparently new, but it is not difficult to prove. From Markov's inequality, $\int_{-\infty }^{0}|w|^{3}dF(w)<\infty $ implies that also $|W|$ has a Pareto tail with index $1/3<\xi <1/2$ in the sense of ((ref)). Part (a) now simply follows from evaluating the sharp bound on the rate of convergence in ((ref)). Part (b) follows from applying the \citeasnoun{Bentkus93} bound ((ref)) to the empirical distribution of $\bar{W}_{i}=W_{i}-n^{-1} \sum_{j=1}^{n}W_{j}$: By ((ref)), $n^{-\xi }\max_{i}|W_{i}|$ converges in distribution, so $\max_{i}|\bar{W}_{i}|=O_{p}(n^{\xi })$. Since $\xi <1/2$, this implies $n^{-1}\sum_{i=1}^{n}\bar{W}_{i}^{2}\mathbf{1}[| \bar{W}_{i}|>\sqrt{n}]\overset{p}{\rightarrow }0$. Furthermore, $|W_{i}|^{3}$ has a Pareto tail of index $3\xi >1$. Thus $n^{-3\xi }\sum_{i=1}^{n}|\bar{W} _{i}|^{3}$ converges in distribution to a stable distribution (see, for instance, \citeasnoun{LePage81}, who elucidate the connection between extreme value theory and stable limit laws), so that $n^{-3/2}\sum_{i=1}^{n}|\bar{W} _{i}|^{3}=O_{p}(n^{3\xi -3/2})$, and the result follows.
The existence of three moments, corresponding to a tail index of $\xi <1/3$, is necessary to obtain the first term of an Edgeworth expansion that underlies the proof of \citeasnoun{Bloznelis03}. More intuitively, recall that under $\xi <1/3$, the Berry-Esseen bound shows that the central limit theorem has an approximation quality of order $n^{-1/2}$. Now under ((ref)), $\mathbb{P}(W_{1}^{R}>\sigma \sqrt{n})$ is of order $(1-n^{-1/(2\xi )})^{n}\approx n^{1-1/(2\xi )}$. Thus, for $\xi >1/3$, the largest observation is of order $\sqrt{n}$ with a probability that is an order of magnitude larger than $n^{-1/2}$. Non-normal observations of order $\sqrt{n}$ are not negligible in the central limit theorem, so the rare large values of $W_{1}^{R}$ under $\xi >1/3$ are responsible for a deterioration of the central limit theorem approximation compared to the $\xi <1/3$ case (cf. \citeasnoun{Hall04}). But from ((ref)) $W_{1}^{R}$ is of order $n^{\xi }$ in nearly all samples, so the bootstrap approximation misses this effect, and systematically underestimates the heaviness of the tail.
We first discuss the approximate parametric problem in more detail. Under the Pareto tail assumption ((ref)), we find from a straightforward calculation that for large $w$, $m(w)=-\mathbb{E}[W|W\leq w]\approx \sigma ^{1/\xi }w^{1-1/\xi }/(1-\xi ).$ Let $s_{n}^{2}=(n-k)^{-1} \sum_{i=1}^{n-k}(W_{i}^{s}-\bar{W}^{s})^{2}$ be the usual variance estimator from the $n-k$ smallest observations. With $k$ fixed, $s_{n}^{2}$ still converges in probability to the unconditional variance of $W$, $s_{n}^{2} \overset{p}{\rightarrow }\limfunc{Var}[W]$. Since the ultimate test we derive is scale invariant, it is without loss of generality to normalize $ \limfunc{Var}[W]=1$. From the convergence to the joint extreme value distribution in ((ref)), $n^{-\xi }\mathbf{W}^{R}\overset{a}{\sim } \sigma \mathbf{X}$, where we write $\overset{a}{\sim }$ for \textquotedblleft is approximately distributed as.\textquotedblright\ Furthermore, under local alternatives $\mathbb{E}[W]=n^{-1/2}\mu ,$ the t-statistic
computed from $\{W_{i}^{s}\}_{i=1}^{n-k}$ is approximately normal with mean $ \mu -n^{-1/2}m(W_{k}^{R})\approx \mu -n^{-1/2}\sigma ^{1/\xi }(W_{k}^{R})^{1-1/\xi }/(1-\xi )$. Combining these two approximations yields
with $\eta _{n}=\sigma n^{-(1/2-\xi )}$ and $Z\sim \mathcal{N}(0,1)$ independent of $\mathbf{X}$.\ The first $k$ elements of $\mathbf{Y}_{n}$ are the largest $k$ order statistics divided by the denominator of the $k+1$ element $T_{n}^{s}$, so that $\mathbf{Y}_{n}$ is invariant to changes in scale $\{W_{i}\}_{i=1}^{n}\rightarrow \{cW_{i}\}_{i=1}^{n}$ for $c>0$. The approximate parametric model on the right-hand side of ((ref)) treats these as jointly extreme value with scale $\eta _{n}$ and tail index $\xi $, and conditionally normally distributed with some (negative) mean that is a function of $X_{k}$ and the parameters $\eta _{n}$ and tail index $\xi $ under $\mu =0$.
As discussed in the introduction, the central idea of this paper is to use the parametric model $\mathbf{Y}_{n}^{\ast }$ to determine a level $\alpha $ test $\varphi : \mathbb{R} ^{k+1}\mapsto \{0,1\}$ of $H_{0}:\mu =0$ that satisfies $\mathbb{E}[\varphi ( \mathbf{Y}_{n}^{\ast })]\leq \alpha $ by construction for all $\xi <1/2$, at least for all $n\geq n_{0}$ and some appropriate upper bounds on $\sigma $. We discuss the construction of such tests in the next section. Any such test $\varphi $ may then be applied to the left-hand side of ((ref)), $ \varphi (\mathbf{Y}_{n})$, to test $H_{0}:\mathbb{E}[W]=0$ from the observations $W_{1},\ldots ,W_{n}$.
Our main theoretical result is the following.
Recall from Theorem (ref) (a) above that the rate of convergence of the normal approximation to the distribution of the t-statistic is $ n^{1/(2\xi )-1}$. Since for $\frac{1+k}{1+3k}<\xi <1/2$, $r_{k}(\xi )>1/(2\xi )-1$, the theorem shows that the difference in the rejection rates of $\varphi $ in the parametric model $\mathbb{E}[\varphi (\mathbf{Y} _{n}^{\ast })]$ and in the original inference for the mean problem $\mathbb{E }[\varphi (\mathbf{Y}_{n})]$ is of smaller order. In this sense, the new approximation provides a refinement for underlying populations that have between two and three moments.
The \citeasnoun{Bentkus93} bound ((ref)) implies that conditional on $ \mathbf{W}^{R}$, $T_{n}^{s}$ is well approximated by a standard normal distribution, since the $W_{i}^{s}$ form an i.i.d. sample from a truncated distribution with less heavy tails compared to the original population. Furthermore, under Condition 1, it follows from ((ref)) that the distribution of $\mathbf{W}^{R}/\sqrt{(n-k)}$ is well approximated by the distribution of $\eta _{n}\mathbf{X}$. The difficulty in the proof of Theorem (ref) arises from the presence of $s_{n}^{2}$ in the scale normalization of $\mathbf{W}^{R}$ in $\mathbf{Y}_{n}$. While it is easy to show that $s_{n}^{2}\overset{p}{\rightarrow }\limfunc{Var}[W]=1$, the proof of Theorem (ref) requires this convergence to be sufficiently fast, and this complication leads to the presence of $k$ in the rate $r_{k}$ (intuitively, larger $k$ lead to more truncation, so $s_{n}^{2}$ is estimated from a distribution with a lighter tail), and the technical requirements on the form of $\varphi $.
Note that the usual full sample t-statistic of $H_{0}:\mathbb{E}[W]=0$, $ T_{n}$, is approximated in terms of $\mathbf{Y}_{n}=(Y_{1n},\ldots ,Y_{kn},T_{n}^{s})^{\prime }$ by
up to a $O_{p}(n^{-1/2})$ term. Under ((ref)), the distribution of $ \tilde{T}_{n}$ is approximated by
Application of Theorem (ref) with $\varphi (\mathbf{Y}_{n}^{\ast })= \mathbf{1}[T(\mathbf{Y}_{n}^{\ast })<t]$ shows that this approximation has a faster rate of convergence compared to the usual standard normal approximation. \citeasnoun{Mueller19} shows that one can combine extreme value theory to improve the rates of approximation to sums of i.i.d. random variables compared to the central limit theorem under $\xi >1/3$. On implication of Theorem (ref) is thus a corresponding result for the case of self-normalized sums ((ref))\ and ((ref)).
In principle, one could use this implication also to construct an alternative test $\varphi $ that simply amounts to a t-test with appropriately increased critical value to ensure size control in the approximate model, $\mathbb{E}[\varphi (\mathbf{Y}_{n}^{\ast })]\leq \alpha $ . This is woefully inefficient, however, since the much larger critical value is only needed for samples where $\xi $ and $\eta _{n}$ are large, which would defeat the objective of obtaining a test that remains close to efficient for populations with thin tails.
To obtain accurate approximations in small samples also for potentially thin-tailed distributions, it makes sense to extend the parametric approximation to populations with an approximate generalized Pareto tail. The c.d.f. $F$ of such populations satisfies
for all $w$ close to the upper bound of the support of $F$, and here and in the following, expressions of the form $(1+\xi x)^{-1/\xi }$ are understood to equal $e^{-x}$ for $\xi =0$. The Pareto tail assumption ((ref)) of Section (ref) is recovered as a special case for $\xi >0$ with $\nu =1/\xi $ and $\sigma $ rescaled by $\xi $.\footnote{ To avoid notational clutter, this section redefines some of the notation previously introduced in Sections 2 and 3 as appropriate for the more general model.}
Assumption ((ref))\ accommodates infinite support thin-tailed distributions, such as the exponential distribution, with $\xi =0$, as well as distributions with finite upper bound on their support, such as the uniform distribution with $\xi =-1$. From the seminal work of \citeasnoun {Balkema74} and \citeasnoun{Pickands75} (also see Theorem 5.1.1 of \citeasnoun{Reiss89} ), it follows that under an appropriate formalization of ((ref)), there exist real sequences $a_{n}$ and $\kappa _{n}$ such that
is (jointly) generalized extreme value distributed, so that $\{(\xi X_{j}+1)^{-1/\xi }\}_{j=1}^{k}\sim \{\sum_{l=1}^{j}E_{l}\}_{j=1}^{k}$ with $ E_{l}$ i.i.d. exponential random variables. If $F$ is exactly generalized Pareto in the sense of ((ref)), then Corollary 1.6.9 of \citeasnoun{Reiss89} \ implies
with $a_{n}=\sigma n^{\xi }$ and $\xi \kappa _{n}=1+n^{-\xi }(\xi \nu -1)$, so that $\sum_{l=1}^{n+1}E_{l}/n\approx 1$ is the only approximation involved in ((ref)). If only the right tail of $F$ of mass $p_{R}>0$ is exactly generalized Pareto, then ((ref)) holds conditionally on the event $\sum_{l=1}^{j}E_{l}/\sum_{l=1}^{n+1}E_{l}\leq np_{R}$, whose probability is larger than 99% for $k=8$ and all $p_{R}\geq 16/n$, $n\geq 50$.
Under ((ref)) and ((ref)), from the same logic that led to ( (ref)), we obtain the approximate model
where $\mathbf{e}$ is a $k\times 1$ vector of ones, $\eta _{n}=n^{-1/2}a_{n}$ and
With $\kappa _{n}\rightarrow 1/\xi $ for $\xi >0$, it is tempting to employ the additional approximation $\kappa _{n}=1/\xi $ to eliminate the location parameter in ((ref)), and this is implicitly applied in standard extreme value theory as reviewed in Section (ref). However, unless $ n $ is very large, this leads to a considerably deterioration of the approximation in ((ref)), and hence ((ref)), so we do not do so in the following.
For practical implementations it is important to allow for the possibility that both tails are potentially moderately heavy. This is straightforward under an assumption that also the left-tail of $F$ is approximately generalized Pareto in the sense of ((ref)): Let $\mathbf{W}^{L}$ be the set of smallest $k$ order statistics. Further let $W_{i}^{m}$ be the $ n-2k$ \textquotedblleft middle\textquotedblright\ order statistics $ k+1,\ldots ,n-k-1$, and let $s_{n}^{2}$ be the sample variance of $W_{i}^{m}$ . Then in analogy to ((ref)),
where $\mathbf{X}^{L}$ and $\mathbf{X}^{R}$ are independent and generalized extreme value distributed with tail index $\xi ^{L}$ and $\xi ^{R}$, respectively, and independent of $Z\sim \mathcal{N}(0,1)$.
The scale and location parameters $\eta _{n}$ and $\kappa _{n}$ in this generalized model depend on the known sample size $n$. But they also depend on the tail parameters of the underlying population: Recall that $\eta _{n}=n^{-1/2}a_{n}=\sigma n^{\xi -1/2}$ in ((ref)). With $ \sigma $ unknown, this product can in principle take on any positive value, even with $(n,\xi )$ known, and the same holds for the parameter $\kappa _{n} $. For this reason, we will now drop the index $n$ in the nuisance parameter $\theta =(\kappa ^{L},\eta ^{L},\xi ^{L},\kappa ^{R},\eta ^{R},\xi ^{R})\in \Theta _{0}$ and in the $2k+1$ dimensional observation $\mathbf{Y} ^{\ast }\mathbf{=Y}_{n}^{\ast }$ from the approximate parametric model in ( (ref)). In this notation, the problem becomes the construction of a powerful test $\varphi (\mathbf{Y}^{\ast })$ of $H_{0}:\mu =0$ against $ H_{a}:\mu \neq 0$ that satisfies
where
and $\mathbf{Y}^{J}=(Y_{1}^{J},\ldots ,Y_{k}^{J})^{\prime }$ for $J\in \{L,R\}$.
From the representation of the joint generalized extreme value distribution in terms of i.i.d. exponentially distributed random variables, it follows that the density of $\mathbf{Y}^{\ast }$ is given by
where $\theta ^{J}=(\kappa ^{J},\eta ^{J},\xi ^{J})$, $\phi $ is the density of a standard normal, $M^{\ast }(\mathbf{y},\theta ^{S})=\eta m^{\ast }( \mathbf{y/}\eta -\mathbf{e}\kappa ,\kappa ,\xi )$, and $f_{T}$ is the \textquotedblleft tail\textquotedblright\ density
with $\theta ^{S}=(\kappa ,\eta ,\xi )$ the parameter of a \textquotedblleft single tail\textquotedblright\ and $x_{i}=y_{i}/\eta -\kappa $ in obvious notation.
Allowing for arbitrary values of the location and scale parameters in the testing problem ((ref)) is not fruitful: An unreasonably large nuisance parameter space $\Theta _{0}$ leads to excessively conservative inference, and it renders the computational determination of powerful tests prohibitively difficult. With that in mind, in the default construction, we consider a nuisance parameter space $\Theta _{0}$ that is partially motivated by a desire to obtain good size control in samples from a demeaned Pareto population when $n\geq n_{0}=50$. In the description of $\Theta _{0}$ , we refer to the extreme value approximation extended to the most extreme $ n_{0}$ observations, $\{W_{i}^{J}\}_{i=1}^{n_{0}}\overset{a}{\sim } \{Y_{i}^{J}\}_{i=1}^{n_{0}}$ for $J\in \{L,R\}$. As noted in ((ref)), this remains an good approximation for an exact generalized Pareto population even if $n=n_{0}$.
In particular, for $J\in \{L,R\}$, we impose
Restriction (a) imposes that the tails are such that at least two moments exists. Restriction (b) says that any potential tail shift is only inward relative to the non-demeaned Pareto default. Note that arbitrarily large inward shifts are incompatible with the population having mean zero. Restriction (c) puts a corresponding lower bound on the inward shift: For the right tail, it requires that the sum of the largest $n_{0}$ observations still has positive mean. To motivate restriction (d), note that the normalization by $s_{n}$ implies that the sum of squared demeaned middle observations cannot be larger than unity. Ignoring the demeaning, taking expectations and approximating the distribution of these observations by again extending the extreme value distribution yields restriction (d) with a right-hand side of unity. We relax the upper bound to equal 2 to accommodate approximating errors in this argument.
We further impose cross restrictions between the two tails:
Restriction (e) amounts to an assumption that the two tails don't overlap. Under an extended tail assumption up to the most extreme $n_{0}/2$ observations, the middle observations take on values between $ -Y_{n_{0}/2}^{L}$ and $Y_{n_{0}/2}^{R},$ leading to restrictions (f)-(g) under the null hypothesis of the overall mean being zero. Finally, restriction (h) is the analogous version of restriction (d) for each tail.
While restriction (c) involves the extreme value approximation for the most extreme $n_{0}$ observations, note that this approximation is only used to motivate a lower bound on $\kappa ^{J}$, and for no other purpose. Consider, for instance, a sample of size $n=n_{0}=50$ from a mean-zero population with a Pareto right tail and a uniform left tail, with overall continuous density. Since the uniform distribution is relatively more spread out compared to the left-tail of a demeaned Pareto distribution, the right tail is shifted outward compared to a demeaned Pareto distribution. Thus, restriction (c) is satisfied for this population, and as long as $k$ is smaller than $n_{0}/2=25$, the approximate parametric model ((ref) ) can still be a good approximation
At the same time, one might argue that if the sample size $n$ is much larger than $n_{0}$, this default parameter space $\Theta _{0}$ is artificially large, and more powerful inference could be obtained by suitably reducing it. Note, however, that for any sample size $n$, the tails could be as large as they are in a sample of size $n_{0}=50$. For instance, consider a sample of size $n=5000$ from a population that is a mixture between a point mass at zero and a demeaned Pareto distribution, with 99% mass on the point mass at zero. Then only approximately 50 observations in the sample will be non-zero, and those follow the demeaned Pareto distribution, so $\Theta _{0}$ is again appropriate, and mechanical reduction of $\Theta _{0}$ as a function of $n$ leads to a poorly performing test in this problem.
Ultimately, inference about the mean requires a substantial assumption about the tails, and the stronger the assumptions, the more powerful the potential inference. The restriction $\Theta _{0}$ as described here is one such choice, and as will be shown below, it yields informative inference while maintaining a high degree of robustness under moderately heavy tails.
Our approach is a variant of the algorithms in \citeasnoun{Elliott15}, denoted by EMW in the following, and \citeasnoun{MuellerWatsonLRCov}; see \citeasnoun{Mueller20} for a detailed survey. This approach yields a likelihood ratio-type test of the form
where the $M$ values of $\theta _{i}\in \Theta _{0}$ and associated positive weights $\lambda _{i}$ are iteratively determined so that the discrete mixture of $\theta $ taking on the values $\theta _{i}$ with probability $ \lambda _{i}/\sum_{j=1}^{M}\lambda _{j}$ forms an approximate least favorable distribution for testing $H_{0}$ against the alternative $H_{a}:$ \textquotedblleft the density of $\mathbf{Y}^{\ast }$ is $f_{a}$ \textquotedblright . Here $f_{a}$ is chosen to equal $f_{a}(\mathbf{y}^{\ast })=\int f(\mathbf{y}^{\ast }|\theta ,\mu )dF_{a}(\theta ,\mu )$ for some weighting function $F_{a}$ that determines against what kind of alternatives the resulting test is designed to be particularly powerful.
We choose $F_{a}(\theta ,\mu )$ to be an improper weighting function with density that is proportional to
so that the implied density on $\mu $, $\kappa ^{L}$ and $\kappa ^{R}$ is flat. This choice is numerically convenient, as it leads to the product form $f_{a}(\mathbf{y}^{\ast })=f_{a}^{S}(\mathbf{y}^{L})f_{a}^{S}(\mathbf{y} ^{R}) $ with $f_{a}^{S}(\mathbf{y)=}\int_{-1/2}^{1/2}f_{a|\xi }^{S}(\mathbf{y }|\xi )d\xi $ and $f_{a|\xi }^{S}$ proportional to the density of the scale and location maximal invariant considered in \citeasnoun{Muller17}. By the same arguments as employed there, $f_{a|\xi }^{S}$ can be obtained by one dimensional Gaussian quadrature, and we approximate $f_{a}^{S}$ by an average of those over a grid of values for $\xi \in \lbrack -1/2,1/2]$.
The lower bound of $-1/2$ on $(\xi ^{L},\xi ^{R})$ in ((ref)) plays no important role, since for values of $\xi ^{J}$ that imply a thin tail, with very high probability the test is constrained to be of a form that does not involve $f_{a}$, as discussed next.
A key ingredient in the algorithm of EMW is an importance sampling estimate of the null rejection probability $\limfunc{RP}(\theta )=\mathbb{E}_{\theta }[\varphi ^{c}(\mathbf{Y}^{\ast })]=\int \varphi ^{c}(\mathbf{y}^{\ast })f( \mathbf{y}^{\ast }|\theta ,0)d\mathbf{y}^{\ast }$ of a candidate test $ \varphi ^{c}: \mathbb{R} ^{2k+1}\mapsto \{0,1\}$ under $\theta ,$
where $\mathbf{Y}_{(l)}^{\ast }$, $l=1,\ldots ,N$ are i.i.d. draws from the proposal density $\bar{f}$ (so that by the LLN, $\widehat{\limfunc{RP}} (\theta )\rightarrow \mathbb{E}_{\bar{f}}[\varphi ^{c}(\mathbf{Y}^{\ast })f( \mathbf{Y}^{\ast }|\theta ,0)/\bar{f}(\mathbf{Y}^{\ast })]=\limfunc{RP} (\theta )$ in obvious notation), where an appropriate $\bar{f}$ may be obtained by the algorithm in \citeasnoun{MuellerWatsonLRCov}. Clearly, the larger $ \Theta _{0}$, the larger the number of importance sampling draws $N$ needs to be for $\widehat{\limfunc{RP}}$ to be of satisfactory accuracy uniformly in $\theta \in \Theta _{0}$.
The nuisance parameter space $\Theta _{0}$ of the last section is unbounded: the restrictions there did not put any lower bound on the scale parameters $ \eta ^{J}$ or the shape parameters $\xi ^{J}$, $J\in \{L,R\}$. Since the distribution of $\mathbf{Y}^{\ast }$ is highly informative about the scale of the tails, it is not possible to obtain uniformly accurate approximations via $\widehat{\limfunc{RP}}$ over $\Theta _{0}$, even with arbitrary computational resources. It is therefore necessary to choose the test $ \varphi $ in a way that does not require a computational check of $\mathbb{E} _{\theta }[\varphi (\mathbf{Y}^{\ast })]\leq \alpha $ over the entirety of $ \Theta _{0}$.
The solution to this challenge suggested in \citeasnoun{Elliott15} is to switch to a default test with known size control under $\Theta _{00}\subset \Theta _{0}$, where the switching rule is such that the default test is employed with probability very close to one whenever $\mathbf{Y}^{\ast }$ is generated from $\Theta _{00}$. If for $J\in \{L,R\}$, $\xi ^{J}$ is very small or $\kappa ^{J}$ is very large, then the resulting observation $Y$$^{J}$ is highly compressed in the sense that $Y_{1}^{J}/Y_{k}^{J}$ is positive and not much larger than one. Also, if $Y_{1}^{J}$ (and thus the entire vector $\mathbf{Y}^{J}$) is small, then even if the tail is heavy in the sense of $Y_{1}^{J}/Y_{k}^{J}$ being large, the tail still only makes a minor contribution to the overall variation of the data. We operationalize this by introducing the switching index
for positive values of $\rho _{r}$, $\rho _{1}$ close to zero, so that $\chi (\mathbf{Y}^{J})=0$ implies that either $Y_{1}^{J}$ is small or $ Y_{1}^{J}/Y_{k}^{J}$ is close to unity. If $\chi (\mathbf{Y}^{J})=0$ for one tail, but not the other, then the problem is heuristically close to knowing that only one of the tails is moderately heavy. For example, suppose $\chi ( \mathbf{Y}^{L})=0$, so the left-tail seems thin. Under approximation ((ref)), the sum of all observations that are not in the right tail equals $\tilde{Y}_{0}^{L}=Y_{0}-\sum_{i=1}^{k}Y_{i}^{L}$, with corresponding approximate variance equal to $\tilde{V}^{L}=1+\sum_{i=1}^{k}(Y_{i}^{L})^{2}$ . It hence makes sense to switch to a \textquotedblleft single tail\textquotedblright\ test $\varphi ^{S}: \mathbb{R} ^{k+2}\mapsto \{0,1\}$ that treats $\mathbf{Y}^{R}$ as the extreme observations from the potentially heavy tail, and $\tilde{Y}_{0}^{L}$ to be approximately normal with mean $-M^{\ast }(\mathbf{Y}^{R},\theta ^{R})$ and variance $\tilde{V}^{L}$. In analogy to ((ref)), such a test is of the form
where $\theta _{i}^{S}\in \mathbb{R} ^{3}$ and $\lambda _{i}^{S}>0$ form again a numerically determined approximate least favorable distribution, and
The corresponding test for a thin right tail is given by $\varphi ^{S}( \mathbf{y}^{L}\mathbf{,y}^{R},-y_{0})$.
If both tails seem thin, then one would expect that the distribution of the analogue to the full-sample t-statistic (cf. ((ref)) from Section (ref))
to be reasonably well approximated by a standard normal distribution, especially if $\sum_{i=1}^{k}(Y_{i}^{R})^{2}+\sum_{i=1}^{k}(Y_{i}^{L})^{2}$ is small.
These considerations, a numerical analysis, and the sequential structure of the eventual algorithm presented in Section (ref) below motivate the restriction of tests $\varphi $ to reject, $\varphi (\mathbf{Y} ^{\ast })=1$, only if all of the following four conditions hold:
where
The additional term $\exp [5\chi (\mathbf{y}^{L})]$ in ((ref)) compared to ((ref)) ensures that conditions 2 and 3 are not binding whenever the corresponding switching index $\chi (\mathbf{y}^{L})$ is large. Its continuity in $\mathbf{y}^{L}$ avoids the sharp change of the form of the rejection region as a function of $\mathbf{y}^{\ast }$ that would be induced by a simpler hard threshold rule $\varphi _{\chi }^{S}(\mathbf{y}^{R} \mathbf{,y}^{L},y_{0})=\mathbf{1}[\chi (\mathbf{y}^{L})=0]\varphi ^{S}( \mathbf{y}^{R}\mathbf{,y}^{L},y_{0})$.
Condition 1 implies that $\varphi $ never rejects if the analogue $T(\mathbf{ Y}^{\ast })$ of the usual t-statistic does not reject; in that sense, we seek to \textquotedblleft robustify\textquotedblright\ the usual t-statistic to obtain better size control. Condition 1 has the additional appeal that sums of the form ((ref)) then effectively only involve $\mathbf{Y} _{(l)}^{\ast }$ for which $|T(\mathbf{Y}_{(l)}^{\ast })|\geq \func{cv}_{T}( \mathbf{Y}_{(l)}^{\ast })$, with an associated gain in computing speed.
We stress that the definition of a \textquotedblleft thin tail\textquotedblright\ in ((ref)), the approximate normality of ((ref))\ and so forth are purely heuristic and do not enter the evaluation of $\mathbb{E}_{\theta }[\varphi (\mathbf{Y}^{\ast })]$ by the algorithm; this probability is always computed from the distribution ((ref)) of $ \mathbf{Y}^{\ast }$. The heuristics merely motivate the particular form of $ \varphi $ just described. As discussed, it is not possible to numerically check that $\mathbb{E}_{\theta }[\varphi (\mathbf{Y}^{\ast })]\leq \alpha $ for all $\theta \in \Theta _{0}$. But we employ extensive numerical analysis to ensure that $\mathbb{E}_{\theta }[\varphi (\mathbf{Y}^{\ast })]\leq \alpha $ over a very large set, including values of $\theta $ that lead to the events $\chi (\mathbf{Y}^{J})=0$ for $J\in \{L,R\}$ with probability close to zero, close to one or in between. The simple form that $ \varphi $ takes on with very high probability in the remainder of the parameter space makes it plausible that $\mathbb{E}_{\theta }[\varphi ( \mathbf{Y}^{\ast })]\leq \alpha $ for all $\theta \in \Theta _{0}$, or at the least, very nearly so.
Even though the switching rule of the last subsection reduces the numerically relevant parameter space to a bounded set, this set still turns out to be so large that a very large number $N$ of importance sampling draws are necessary to obtain adequate approximations. The computationally expensive part in the evaluation of $\widehat{\limfunc{RP}}(\theta )$ in ( (ref)) for different $\theta $ is the evaluation of $f(\mathbf{Y} _{(l)}^{\ast }|\theta ,0)$ (since all $\bar{f}(\mathbf{Y}_{(l)}^{\ast })$ can be computed once and stored).
These evaluations can be dramatically sped up by recombining two \textquotedblleft single tails\textquotedblright\ in different combinations: For a given $\theta ^{S}=(\kappa ,\eta ,\xi )$, let $\mathbf{Y}^{e}\in \mathbb{R} ^{k+1}$ be an \textquotedblleft extended\textquotedblright\ single tail with distribution
where $\mathbf{X}$ is distributed as as in ((ref)), independent of $ Z\sim \mathcal{N}(0,1)$. Denote the density of $\mathbf{Y}^{e}$ by $f^{e}( \mathbf{y}^{e}|\theta ^{S})$. Given two independent vectors $\mathbf{Y} _{(1)}^{e}$ and $\mathbf{Y}_{(2)}^{e}$ distributed according to $\theta _{1}^{S}=\theta ^{L}$ and $\theta _{2}^{S}=\theta ^{R}$, respectively, note that their combination into the \textquotedblleft both tails\textquotedblright\ observation $(\mathbf{Y}_{(1)}^{S\prime },\mathbf{Y} _{(2)}^{S\prime },Y_{0,(1)}^{e}-Y_{0,(2)}^{e})^{\prime }\in \mathbb{R} ^{2k+1}$ has the same distribution as $\mathbf{Y}^{\ast }$ in ((ref)), since the difference of two independent normals of variance 1/2 is again standard normal. Thus, with $\mathbf{Y}_{(l)}^{e}$ i.i.d. draws from a suitable proposal density $\bar{f}^{e}$, one obtains the alternative estimator
that recombines each extended single tail with $K$ different other extended single tails, for a total of $KN$ importance draws. Yet evaluation of ((ref)) only requires a simple product of the $(K+N)$ values $f^{e}( \mathbf{Y}_{(l)}^{e}|\theta ^{S})$ for $\theta ^{S}\in \{\theta ^{L},\theta ^{R}\}$. We let $K=128$ and $N=640,000$ for a total of nearly 82 million importance sampling draws.
The overall algorithm proceeds in four stages. To describe these stages, let $\Theta _{0}^{S}\subset \mathbb{R} ^{3}$ be the set of parameters satisfying the constraints (a)-(d) of Section (ref) on one tail. Let $\Theta _{s}^{S}\subset \Theta _{0}^{S}$ be such that for $\theta ^{S}\in \Theta _{s}^{S}$, the event that the switching index is zero, $\chi (\mathbf{Y}^{J})=0$, happens with at least 90% probability, and $\Theta _{ss}^{S}\subset \Theta _{s}^{S}$ be such that $\chi (\mathbf{Y}^{J})=0$ with probability of exactly 90%. Note that $ \Theta _{s}^{S}$ and $\Theta _{ss}^{S}$ depend on $(\rho _{r},\rho _{1})$ in ((ref)).
Note that the parameter set under consideration becomes consecutively larger in Steps 1-3, and the form imposed on tests described in Section (ref) ensures that any potential remaining overrejections of a stage can be corrected by the subsequent stage, which increases the numerical stability of the algorithm. Null rejection probabilities are estimated throughout with the importance sampling estimator of Section (ref). This estimator has an importance sampling standard error (appropriately adjusted for the dependence in ((ref))) of no more than $0.05\%$, $0.15\%$ and $0.2\%$ for $\alpha =1\%$, $5\%$, 10%, respectively.
We apply this algorithm to the default nuisance parameter space with $ n_{0}=50$ of Section (ref) for $k=8$ and $\alpha \in \{0.002,0.004,\ldots $, $0.008,0.01,0.02,\ldots $, $0.10,0.12,\ldots $, $ 0.20,0.25,0.30,0.40,0.50\}$, and also for $k=4$ to the larger nuisance parameter space where $n_{0}=25$ in the constraints of Section (ref). We ensure that the 95% and 99% level confidence intervals obtained via test inversion\footnote{ In the rare samples where test inversion yields disconnected sets, we set the confidence interval equal to the smallest interval that contains all non-rejections.} always contain the 90% and 95% level intervals, respectively, and that the p-value is always coherent with the level of the reported confidence interval by adding the obvious additional constraints to the form of the tests for $\alpha \neq 0.05$. After trivial modifications that decrease their rejection probability by an arbitrarily small amount, these tests satisfy the condition of the two-tailed analogue of Theorem (ref); see Appendix (ref) for details. For comparison purposes, we also generate tests with $k\in \{4,12\}$ for $\alpha \in \{0.01,0.05\}$ in the default parameter space. For $k=8$ and a given level $ \alpha $, the computations take about one hour on a modern 24 core workstation in a Fortran implementation, and about 3 hours for $k=12$. Once the values for $\{\lambda _{i}^{S},\theta _{i}^{S}\}_{i=1}^{M_{S}}$ and $ \{\lambda _{i},\theta _{i}\}_{i=1}^{M}$ are determined in this fashion, the evaluation of the resulting test $\varphi $, as required in applications, is computationally trivial.
Suppose we estimate the parameter vector $\vartheta =(\beta ,\gamma ^{\prime })^{\prime }\in \mathbb{R} ^{q}$ by Hansen's (1982)\nocite{Hansen82} Generalized Method of Moments using the $r\times 1$ moment condition $\mathbb{E}[g(\vartheta ,z)]=0$ from data $z_{j}$, $j=1,\ldots ,n_{z}$ and $r\times r$ positive definite weighting function $\hat{\Psi}$. Suppose further that the data $z_{j}$ is i.i.d. across clusters defined by the partition $\{\mathcal{C} _{i}\}_{i=1}^{n}$ of $\{j:1\leq j\leq n\}$ (so that $\mathcal{C}_{i}=\{i\}$ and $n=n_{z}$ under i.i.d. sampling of $z_{j}$). Then as $n\rightarrow \infty $, under standard regularity conditions, $\hat{\vartheta}=(\hat{\beta} ,\hat{\gamma}^{\prime })^{\prime }$ satisfies
where $G_{i}=\sum_{j\in \mathcal{C}_{i}}g(\vartheta ,z_{j})$ are i.i.d., $ \hat{\Gamma}=n^{-1}\sum_{j=1}^{n_{z}}\partial g(\vartheta ,z_{j})/\partial \vartheta ^{\prime }|_{\vartheta =\hat{\theta}}\overset{p}{\rightarrow } \Gamma $ and $\hat{\Psi}\overset{p}{\rightarrow }\Psi $ with $\Gamma $ and $ \Psi $ non-stochastic, so that the large sample variability of $\hat{ \vartheta}$ is entirely driven by the average of i.i.d. observations $G_{i}$ . Correspondingly, the standard GMM hypothesis test of $H_{0}:\beta =\beta _{0}$ is numerically equivalent to the usual t-test of $H_{0}:\beta =\beta _{0}$ computed from the $n$ observations
where $\hat{G}_{i}=\sum_{j\in \mathcal{C}_{i}}g(\hat{\vartheta},z_{j})$ and $ \iota _{1}$ is the $q\times 1$ vector $(1,0,\ldots ,0)^{\prime }$.
Thus, to the extent that $G_{i}$ follows a distribution with moderately heavy tails, one would expect that small sample inference is improved by applying the new test to the observations $\{\hat{W}_{i}\}_{i=1}^{n}$. A corresponding analytical refinement result analogous to Theorem (ref) is beyond the scope of this paper.
This section presents six sets of small sample results: two for inference about the mean from an i.i.d. sample, two for the difference of population means from two independent samples, and two for a regression coefficient with clustered standard errors. In all three cases, the data is either generated from analytical distributions, or from draws with replacement from a large data set. We focus on tests of nominal 5% level in the main text; results for 1% level tests are reported in the appendix and exhibit broadly similar patterns.
We initially compare our default test with $k=8$ (\textquotedblleft new default\textquotedblright ) with three standard tests for the population mean: standard t-statistic based inference with critical value from a student-t distribution with $n-1$ degrees of freedom \textquotedblleft t-stat\textquotedblright ; the percentile-t bootstrap based on the absolute value of the t-statistic \textquotedblleft sym-boot\textquotedblright ; and the percentile-t bootstrap based on the signed t-statistic \textquotedblleft asym-boot \textquotedblright . The data is generated from one of seven populations: the standard normal distribution N(0,1), the log-normal distribution LogN, the F-distribution with 4 degrees of freedom in the numerator and 5 in the denominator F(4,5), the student-t distribution with 3 degrees of freedom t(3), an equal probability mixture between a N(0,1) and LogN distribution Mix1, and a 95 / 5 mixture between a N(0,1/25) and a LogN distribution Mix2. All population distributions are normalized to have mean zero and unit variance; the corresponding densities are plotted in Figure 1. Technically, only the Pareto distribution, the t-distribution and the F-distribution exhibit heavy Pareto like tails in the sense of ((ref)) with tail indices $\xi =0.4$, $\xi =1/3$ and $\xi =0.4$, respectively, but as a practical matter, also the log-normal and the two mixture distributions are right-skewed enough to make small sample inference challenging.
Table 1 reports null rejection probabilities, along with the average length of the resulting confidence interval, expressed as a multiple of the average length of the infeasible confidence interval that is based on the t-statistic, but applies the size adjusted critical value. As can be seen from Table 1, the new method comes much closer to controlling size under moderately heavy-tailed distributions. For the thin-tailed normal population, the new method only leads to 10% longer intervals for $n=50$, and essentially no excessive length for $n\in \{100,500\}$. For other populations, the intervals of the new method are often much longer than those from other methods; but since the other methods do not come close to controlling size, that comparison is not meaningful (entries in bold indicate where tests are close to valid with a null rejection probability below 6%). Remarkably, for $n=50$, the new method yields shorter confidence intervals than the size corrected t-statistic for some populations while still controlling size. The explicit modelling of the tails can also yield efficiency gains, since under a Pareto-like tail, the sample mean is not the efficient estimator of the population mean.
An exception to the good performance of the new method is the student-t population with three degrees of freedom. Even though it has fairly heavy tails, with the third moment not existing, its symmetry enables t-stat and sym-boot to control size at much less cost to average length compared to the new method.\footnote{ The analytical result by \citeasnoun{Bakirov05} shows that the usual 5% level t-test remains small sample valid under arbitrary scale mixtures of normals, which includes all t-distributions.}
Table 2 compares different versions of the new method across the same set of seven populations. We consider $k\in \{4,8,12\}$ for the default parameter space with $n_{0}=50$, and also include the even more robust test with $k=4$ constructed from the parameter space in Section (ref) with $ n_{0}=25$. For $n=25,$ only the test with $n_{0}=25$ comes close to controlling size for non-thin tailed populations (and the test with $k=12$ cannot be applied at all, since there is only a single \textquotedblleft middle\textquotedblright\ observation). For larger $n$, the tests with $k=4$ are even more successful in controlling size compared to the default method, but at a non-negligible cost in terms of longer confidence intervals. In contrast, the test for $k=12$ does not yield an additional substantial reduction in average length, and has worse size control for $n=50$. These results underlie our choice of the test with $k=8$ and $n_{0}=50$ as the default, and in the following, we exclusively focus on this variant.
A potential objection to this first set of Monte Carlo results is that the underlying populations have smooth tails, which might overstate the effectiveness of the new method \textquotedblleft in practice\textquotedblright . To address this concern, consider a population that is equal to the (discrete) distribution from a large economic data set. We use the income data of 2016 mortgage applicants as reported by banks under the Home Mortgage Disclosure Act (HMDA). From this database of more than 16 million applications, we create subpopulations that condition on U.S. state and the gender of the applicant, as well as the purpose of the mortgage (home purchase, home improvement or refinancing) and whether or not the unit is owner-occupied. We eliminate all records with missing data, and only retain subpopulations with at least 5000 observations. For each of the resulting 300 subpopulations, we compare the performance of alternative methods for inference about the mean, based on i.i.d. samples of size $n$ (that is, sampling is with replacement).
Panel (a) of Figure 2 plots the cumulative distribution function of the null rejection probabilities over the 300 subpopulations for each tests considered in Table 1, estimated from 20,000 draws from each subpopulation. Nominally, all mass should be to the left of the 5% line, but the traditional tests don't comes close. For instance, for $n=100,$ the usual t-statistic has null rejection probability of less than 10% for only approximately 40% of the 300 subpopulation. In comparison, the new test controls size much more successfully.
Panel (b) of Figure 2 plots the cumulative distribution function of the average length of the confidence intervals, relative to the average length of the size corrected t-statistic based interval. For $n=50$, the new method not only controls size better than the bootstrap tests, but it also leads to confidence intervals that are typically shorter on average. In fact, they are substantially shorter than what is obtained from the infeasible size corrected interval. For $n=\{100,500\}$, this is no longer the case and the better size control of the new method comes at the cost of somewhat longer confidence intervals.
One might argue that in the HMDA example, one could avoid the complications of the heavy right tail of the income distribution by considering the logarithm of the applicants' income. But, of course, there is no robust way to transform a confidence interval for the population mean of log-income into a valid confidence interval for the population mean income. What is more, in many contexts, the policy relevant parameter is the population mean (and not, say, the median) of some potentially heavy-tailed distribution: think of health care costs, or flood damage, or asset returns.
Our second set of Monte Carlo experiment concerns inference about the difference of two population means $\mathbb{E}[W^{\func{I}}]-\mathbb{E}[W^{ \func{II}}]$ based on two independent equal-sized i.i.d. samples $ W_{i}^{j}\sim W^{j}$, $i=1,\ldots ,n/2$, $j\in \{$I$,$II$\}$. Casting this in terms of a linear regression and applying the general mapping ((ref))\ yields
where $\bar{W}^{j}=(n/2)^{-1}\sum_{i=1}^{n/2}W_{i}^{j}$ are the sample means for $j\in \{$I$,$II$\}$.
We initially generate data according to
for $i=1,\ldots ,n/2$, where $\varepsilon _{i}^{j}\sim iid\mathcal{N} (0,1/10) $ across $i$ and $j\in \{$I$,$II$\}$, and $\nu _{i}$ is distributed according to one of the distributions of Table 1. Inference about $\mathbb{E} [W^{\func{I}}]-\mathbb{E}[W^{\func{II}}]$ can then be thought of as inference about the average treatment effect $\mathbb{E}[\nu _{i}]$, with the design amounting to a large but highly heterogeneous additive treatment effect.
Table 3 compares the new method to standard t-statistic based inference and a symmetric and asymmetric percentile-t bootstrap, where now the bootstrap samples combine $n/2$ randomly selected observations with replacement from each of the two samples. In this exercise the design with $\nu _{i}\sim \mathcal{N}(0,1)$ leads to a much longer confidence interval from the new method with $n=50$. The reason is that with $\varepsilon _{i}^{j}\sim \mathcal{N}(0,1/10)$ in ((ref)), $W_{i}^{\func{I}}$ has much larger variance than $W_{i}^{\func{II}}$. The distribution of $\hat{W}_{i}$ is thus approximately equal to a 50-50 mixture of two normal distributions with very different variances, which is heavier tailed than a normal distribution. At the same time, for asymmetric $\nu _{i}$, standard methods do not control size well, while the new method does so much more successfully.
As a second exercise, we generate $W_{i}^{j}$ as $n/2$ i.i.d. draws of two randomly selected subpopulations of the HMDA data set considered in the last section. Inference about $\mathbb{E}[W^{\func{I}}]-\mathbb{E}[W^{\func{II}}]$ then corresponds to inference about the average treatment effect if the treatment induces a change from the distribution of income in one subpopulation to the distribution in another---maybe a plausible calibration for an intervention that affects individuals' incomes. Figure 3 reports the performance of the inference methods of Table 3 for 200 randomly selected pairs of subpopulations, in analogy to Figure 2 above. We find that also in this exercise, standard methods fail to produce reliable inference, while the new method is substantially more successful at controlling size.
A third set of Monte Carlo experiments explores the performance of the new method for inference in a clustered linear regression
with conditionally mean zero $u_{it}$, so that there are $T_{i}$ observations in cluster $i$. Viewing linear regression as a special case of GMM inference, we obtain from the development of Section (ref) and the Frisch-Waugh Theorem that
where $\hat{u}_{it}$ and $\hat{\beta}$ are the OLS\ estimates of $u_{it}$ and $\beta $, and $\hat{X}_{it}$ are the residuals of a OLS regression of $ X_{it}$ on $Z_{it}$. We consider four tests of $H_{0}:\beta =\beta _{0}$: The t-statistic implemented by STATA, which is nearly identical to a standard t-test applied to $\hat{W}_{i}$, except for degree of freedom corrections; the suggestion of \citeasnoun{Imbens16} to account for a potentially small number of heterogeneous clusters \textquotedblleft Im-Ko \textquotedblright\ (we consider the variant that involves the data dependent degree of freedom adjustment $K_{IK}$ in their notation); the wild cluster bootstrap that imposes the null hypothesis suggested by \citeasnoun {Cameron08} \textquotedblleft CGM\textquotedblright ; and the new default test applied to $\hat{W}_{i}$ \textquotedblleft new default \textquotedblright .
We initially consider data generated from model ((ref)) where
$\nu _{i}$ is i.i.d. mean-zero with a distribution that is one of the seven populations considered in Table 1, one element of $Z_{it}$ is a constant, and $X_{it},$ the 5 non-constant elements of $Z_{it}$, and $\varepsilon _{it} $ are independent standard normal. We set $T_{i}=T=10$ for all clusters. The presence of $\nu _{i}$ induces heteroskedastic correlations within each cluster of observations $\{Y_{it}\}_{t=1}^{T}$.
Table 4 reports the results. As in the inference about the mean problem, the new method is seen to control size much more successfully compared to the other methods, although at a cost in average confidence interval length that is more pronounced than in Table 1 for the thin-tailed $\nu _{i}\sim \mathcal{N}(0,1)$. Intuitively, the product of two independent normals $\nu _{i}X_{it}$ has considerably heavier tails than a normal distribution, but it is still symmetric.
In the final Monte Carlo exercise we again consider a discrete population from a large economic data set. In particular, we consider a sample of all employed workers aged 18-65 from the 2018 merged outgoing rotation group sample of the Current Population Survey (CPS). We let the dependent variable $Y_{it}$ be the logarithm of wages, and pick the regressor of interest $ X_{it}$ and the 5 non-constant controls $Z_{it}$ as a random subset of potential regressor including gender, race, age and dummies for Hispanic, non-white, married, public sector employer, union membership and whether hours or the wage was imputed. The resulting coefficient $\beta $ on $X_{it}$ in the regression using the entire 145,838 individuals in the database is the population coefficient. We cluster at the level of 308 Metropolitan Statistical Areas (MSAs).\footnote{ For the purposes of this exercise, we treat as additional MSAs the part of each U.S. state outside of any CBSA area.} The four different methods of Table 4 are then employed to conduct inference about $\beta $ based on a sample consisting of all individuals that reside in $n$ randomly selected MSAs, where the MSAs are drawn with replacement. By construction the clusters are thus i.i.d. and the population regression coefficient is equal to $\beta $.
Figure 4 depicts the results over 200 populations generated in this manner, analogous to Figure 3, for $n\in \{50,100,200\}$. (We consider $n=200$ rather than $n=500$ for the largest sample size to avoid that with high probability, samples contain many identical clusters.) In this design none of the methods come close to perfectly controlling size. Still, the new method is substantially more successful, albeit at the cost of considerably longer average confidence intervals for $n\in \{100,200\}$.
The poor performance of the standard methods might come as a surprise given that none of the variables in the CPS exercise are heavy-tailed, and the number of clusters is not particularly small. Approximately (cf. equation ( (ref))), the variability of the OLS estimator $\hat{\beta}$ is driven by the average of the $n$ i.i.d. random variables
where $u_{it}$ is the population regression error and $\tilde{X}_{it}$ is the residual of a population regression of $X_{it}$ on $Z_{it}$. The distribution of $G_{i}$ may be heavy-tailed because (i) $u_{it}$ has a heavy-tailed component, as in ((ref)) above; (ii) the joint distribution of $(\tilde{X}_{it},u_{it})$ is such that $\tilde{X}_{it}u_{it}$ is heavy-tailed; (iii) $\tilde{X}_{it}$ is heavy-tailed; (iv) $T_{i}$ is heterogeneous across $i$, so that large clusters with big $T_{i}$ lead to $ G_{i}$ with high variance; or a combination of these effects. MSAs are highly heterogeneous in their size: the largest contains 6,163 individuals, and the smallest only 42. Effect (iv) is thus clearly present, and the suggestion by \citeasnoun{Imbens16} is designed to accommodate effects (iii) and (iv). But as reported in Table 4, if $G_{i}$ is heavy-tailed due to effect (i), then the adjustment of \citeasnoun{Imbens16} does not help much. The CPS\ design seems to exhibit all four effects to some degree, making correct inference quite challenging, and the new method relatively most successful at controlling size.
Whenever researchers compare a t-statistic to the usual standard normal critical value they effectively assume that the central limit theorem provides a reasonable approximation. This is true when conducting inference for the mean from an i.i.d. sample, but it holds more generally for linear regression, GMM inference, and so forth. As is well understood, the central limit theorem requires that the contribution of each term to the overall variation is small. To some extent, this is empirically testable: one can simply compare the absolute values of each (demeaned) term with the sample standard deviation. The normal approximation then surely becomes suspect if the largest absolute term is, say, equal to half of a standard deviation.
One may view the new test suggested here as a formalization of this notion: the extreme terms are set apart, and if they are large, then the test automatically becomes more conservative. What is more, even if the sample realization from an underlying population with a heavy tail fails to generate a very large term, it still leaves a tell-tale sign in the large spacings between the largest terms. Correspondingly, the test also becomes more conservative if the largest observations are far apart from each other, even if the largest one isn't all that large---the new method seeks to infer the likelihood of a potential large outlier based on the spacings of the extreme terms. These adjustments are disciplined by an assumption of Pareto-like tails. But as was found in the small sample simulations, they help generate more reliable inference also when the underlying population is more loosely characterized by a moderately heavy tail.
It would be desirable to extend the new method also to $F$-type tests of null hypotheses that restrict more than one parameter. Such an extension is far from straightforward, though: extreme value theory for the joint extremes of vector-valued observations does not yield a tightly parametrized approximate model, and the approach pursued here to determine the appropriate adjustments is already very computationally challenging.