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,801 characters · 10 sections · 38 citation commands
Fast Test Inversion for Resampling Methods
The usual method to produce confidence intervals in randomization-based inference is grid search garthewaite1996. This method exploits the duality between hypothesis testing and confidence intervals. Given a level $\alpha$ test, the confidence interval is the set of $\beta_0$ values for which we fail to reject the null hypothesis. While there are times when grid search is appropriate, it should not be our primary method as it has two major drawbacks.
First, it is computationally inefficient hansenprob. For example, in d2024robust, the authors develop a randomization test and apply their method to construct confidence intervals for two different examples. In the first example, they construct numerous confidence intervals using different methods by inverting “the tests of $\beta = \beta_0$ for $\beta_0 \in \{-1.7, -1.69, \dots, 0.3\}$,” resulting in 201 tests. Their method takes 5 seconds for this simple example, which is a considerable amount of time given that there are only 51 observations. In a second example, d2024robust “invert the tests of $\beta = \beta_0$ for $\beta_0 \in \{-40, -39.99, \dots, 39.99, 40\}$,” resulting in a staggering 8,001 hypothesis tests. This computational burden becomes even more pronounced for the randomization test developed by diciccio:2017aa, which requires approximately 17 hours to complete.\footnote{The supplementary materials estimates that approximately 33 hours are needed altogether to construct the confidence intervals in Table 3, so 17 hours may be conservative.}; this method takes much longer in part because they need to construct a confidence region in $\mathbb{R}^2$ and project onto one dimension, meaning that they need two dimensional grid search. Clearly, there is a need for a faster method to create confidence intervals and regions.
Second, the confidence intervals produced by grid search are conservative. This conservatism arises because the sequence of $\beta_0$ values increases in steps of $\epsilon$, leading to undercoverage if we include only the convex hull of $\beta_0$ values for which we fail to reject. Conversely, if we include adjacent $\beta_0$ values where rejection occurred, the confidence interval will be conservative. One possible alternative is to use bisection search algorithms paired with linear interpolation, but these usually result in approximations garthewaite1996.
We develop a method to construct fast and exact confidence intervals by exploiting the structure of the test and randomization statistics. We define exact confidence intervals as confidence intervals that only include the values of $\beta_0$ for which we fail to reject the null hypothesis. To the best of our knowledge, CaiCanayKimShaikh:2023aa and guan2024 present the first exact confidence intervals that use this approach. We extend their work by considering a much broader class of test statistics, including those where the variance estimator depends on the null hypothesis.
As a proof of concept, we primarily use the tests from Pouliot2024. He develops a test that can handle general design matrices and provides exact inference under homoskedasticity in small samples and asymptotically exact inference under general model misspecification. We will extend the principles applied to Pouliot2024 to other randomization tests—including that developed by d2024robust and diciccio:2017aa—in Section (ref) and the Appendix. Since we generally adopt the notation from each source, readers should consider the notation from each test and its analysis independently unless we explicitly indicate otherwise.
We encourage readers to read Pouliot2024, but to provide some background and perspective, consider the linear regression model $$ \mathbf{Y} = \mathbf{1}_n \beta_0 + \mathbf{X}_1 \beta_1 + \mathbf{X}_2 \beta_2 + \varepsilon, $$ where $\mathbf{Y}, \mathbf{1}_n, \mathbf{X}_1, \varepsilon \in \mathbb{R}^{n \times 1}$; $\mathbf{X}_2 \in \mathbb{R}^{n \times p}$; $\beta_0, \beta_1 \in \mathbb{R}$; and $\beta_2 \in \mathbb{R}^{p \times 1}$.
To test $H_0:\beta_1 = \beta_1^0$ against $H_1:\beta_1 < \beta_1^0$ (see Appendix (ref) for the two-sided hypothesis test), we use a randomization test by comparing the test statistic against the randomization statistics. The randomization statistic is defined as $$ T_g(\beta_1^0) \equiv \frac{(\mathcal{Q}_1 \mathbf{X}_1)^\top g (\mathbf{Y} - \mathbf{X}_1 \beta_1^0)}{\hat{\sigma}_g}, $$ where $$ \hat{\sigma}_g^2 \equiv \frac{1}{n}\sum_{i=1}^{n} (\mathcal{Q}_1 \mathbf{X}_1)_i^2 (\mathcal{Q}_{2} g \mathbf{Y})_i^2; $$ $g$ is a block permutation matrix in $\mathbf{G}$, which is the set of all block permutation matrices for a prespecified number of blocks; $\mathcal{Q}_1$ is the annihilator matrix that projects onto the orthogonal complement of the span of $\left\{g\mathbf{X}_2\right\}_{g \in \mathbf{G}}$; and $\mathcal{Q}_2$ is the annihilator matrix that projects onto the orthogonal complement of the span of $\mathbf{1}$, $\left\{g\mathbf{X}_1\right\}_{g \in \mathbf{G}}$, and $\left\{g\mathbf{X}_2\right\}_{g \in \mathbf{G}}$. Let $M = |\mathbf{G}|$, the total number of block permutations $g$ (this definition applies throughout the paper). The test statistic corresponds to $g = \text{Id}$, the identity matrix.
The full procedure to perform the hypothesis test is detailed in Algorithm (ref). We can also obtain a $p$-value for a given $\beta_1^0$ value by calculating the proportion of randomization statistics that are greater or equal to the test statistic. This is essentially comparing the rank of the test statistic with the randomization distribution.
As noted above, the traditional method to construct confidence intervals is grid search, but this method is computationally inefficient and results in conservative confidence intervals. There are two other challenges with this method here. First, we don't know where the center of the confidence interval is, so we don't know the range of $\beta_1^0$ values to search over. Traditionally, we would create a grid of $\beta_1^0$ values to the left and right of $\hat{\beta}_1$ because the confidence intervals are centered at $\hat{\beta}_1$, but $\hat{\beta}_1$ may not be the center of the interval here.\footnote{Even creating a grid of $\beta_1^0$ values from $\hat{\beta}_1 - 4\text{s.e.}(\hat{\beta}_1)$ to $\hat{\beta}_1 + 4\text{s.e.}(\hat{\beta}_1)$ will not necessarily solve this issue because the confidence intervals can be wide for poor orderings of the data that reduce power (i.e. $\mathcal{Q}_1$ orthogonalizes almost all variation in $\mathbf{X}_1$.} Second, we don't know when to stop searching. The confidence interval from randomization tests can have “holes”—the interval could be a union of two disjoint sets. This arises because the $p$-value function may not be monotonic, so it could be the case that the $p$-value initially decreases as $\beta_1^0$ gets further from the truth, but later increases.
The solution to these woes is to pair the structure of the test statistic with algebra and geometry to find the exact $p$-value function from which we can derive precise confidence intervals for any choice of $\alpha$. To see this, note that the randomization statistics can be reformulated as: $$
$$ where $m_g = -(\mathcal{Q}_1 \mathbf{X}_1)^\top g \mathbf{X}_1 / \hat{\sigma}_g$ is the slope and $b_g = (\mathcal{Q}_1 \mathbf{X}_1)^\top g \mathbf{Y} / \hat{\sigma}_g$ is the intercept for the randomization statistic with block permutation $g$. From this, we see that randomization statistics (including the test statistic) are lines. This means that each of the randomization statistics $T_g(\beta_1^0)$ (excluding the test statistic) will be exactly equal to $T_\text{Id}(\beta_1^0)$ at a unique $\beta_1^0$ value because two lines with differing slopes must intersect at a single point.\footnote{Two lines can intersect at one, zero or an infinite number of points, but the latter two require the slopes to be equal, which has zero probability.} We denote these $\beta_1^0$ values for which the test statistic line intersects the randomization statistic line as $\dot{\beta}_1^{0, g}$. Finding these $\dot{\beta}_1^{0, g}$ values only requires a bit of algebra; we just need to set the test statistic equal to the randomization statistic and solve for $\beta_1^0$: $$
$$ where $m_{Id} \neq m_g$. With $\dot{\beta}_1^{0, g}$ and the slope of the test and randomization statistics, we can easily deduce the relationship between $T_{Id}(\beta_1^0)$ and any randomization statistic $T_g(\beta_1^0)$ for all $\beta_1^0 \in \mathbb{R}$.\footnote{If $m_{Id} = m_g$, then $b_{Id} > b_g \Leftrightarrow T_{Id} > T_g$.} For example, notice that if $m_{Id} < m_g$, then $T_{\text{Id}} > T_g$ for all $\beta_1^0 < \dot{\beta}_1^{0, g}$, and $T_{\text{Id}} < T_g$ for all $\beta_1^0 > \dot{\beta}_1^{0, g}$. These $\dot{\beta}_1^{0, g}$ values are critical because it is only at these specific $\beta_1^0$ values that the rank of the test statistic changes in relation to the randomization distribution. At all other $\beta_1^0$ values, the test and randomization statistics are changing in value, but the rank of the test statistic does not change in relation to the randomization distribution, meaning that the $p$-value does not change. See Figure (ref) for a visualization.
The $p$-value function is $$ p(\beta_1^0) \equiv \frac{1}{M}\left(1 + \sum_{g \in \mathbf{G}_n \setminus \text{Id}} \mathbbm{1}[T_{\text{Id}}(\beta_1^0) \leq T_g(\beta_1^0)]\right). $$ Note that the $p$-value curve is a step function (see Figure (ref) for an example). We can find the curve using Algorithm (ref). The exact level $1-\alpha$ confidence interval is the set of $\beta_1^0$ values with $p$-value greater than $\alpha$, which we can find by filtering the rows of $F$, which is the output of Algorithm (ref).
We also make the observation here that $T_g(\beta_1^0) \in \mathbb{R}[\beta_1^0]$, a scalar polynomial of degree 1 over the field $\mathbb{R}$. This perspective does not necessarily unlock new insights in this simple example, but we nevertheless introduce it here as it will be helpful with more complex randomization statistics like the over-identified instrumental variable case discussed in the following section.
Pouliot2024 also defines a test for overidentified instrumental variable analysis. Suppose there are $k$ instruments and $\beta_1^0$ is a scalar. Then, the heteroskedasticity-robust Wald-like randomization statistic is $$ t_g(\beta_1^0) = T^{iv}_g(\beta_1^0)^T \hat{\Sigma}_g(\beta_1^0)^{-1} T^{iv}_g(\beta_1^0), $$ where $$
$$ $\mathbf{Y}, \mathbf{X}_1 \in \mathbb{R}^{n \times 1}$; $\mathbf{X}_2 \in \mathbb{R}^{n \times p}$; $\mathbf{Z}_2 \in \mathbb{R}^{n \times k}$; $\mathcal{Q}_1$ is defined in Section \ref{onesidedPouliot}; and $\mathcal{Q}_3$ is the annihilator matrix that projects onto the orthogonal complement of the span of $\mathbf{1}$, $\mathbf{X}_1$, and $\left\{g\mathbf{X}_2\right\}_{g \in \mathbf{G}}$. Note that unlike the previous example where the variance estimator $\hat{\sigma}_g^2$ was not a function of $\beta_1^0$, that is not the case here with $\hat{\Sigma}_g(\beta_1^0)$.\footnote{This statistic slightly differs from that in the paper as the variance estimator $\hat{\Sigma}_g$ is a function of $\beta_1^0$. This variance estimator slightly increases power.}
We would like to apply the same principles as before: transform the test and randomization statistics into some form where we can find the $\beta_1^0$ values for which the test statistic intersects the randomization statistic. The issue here is that our covariance matrix estimator—that we need to invert—is sandwiched between two outer "bread" matrices. We cannot set the test and randomization statistics equal to each other and multiply $\hat{\Sigma}_g(\beta_1^0)$ and $\hat{\Sigma}_\text{Id}(\beta_1^0)$ on both sides to remove $\hat{\Sigma}_g(\beta_1^0)^{-1}$ and $\hat{\Sigma}_\text{Id}(\beta_1^0)^{-1}$. We solve this problem by recognizing that $\hat{\Sigma}_g(\beta_1^0)$ is in $\mathbb{R}^{k \times k}[\beta_1^0]$: it is a symmetric polynomial matrix where each element is a quadratic function in $\beta_1^0$. To see this, observe that $\hat{\Sigma}_g(\beta_1^0)$ can be reformulated as $$ \hat{\Sigma}_g(\beta_1^0) =
, $$ where for $1 \leq v,w \leq k$, $$
$$ Here, diag$[\bullet]$ denotes the matrix formed from vector $\bullet$, and $\circ$ denotes the Hadamard (element-wise) product (these definitions apply throughout the paper). Then, $$
$$
$\hat{\Sigma}_g(\beta_1^0)$ is invertible if the determinant is a non-zero polynomial kaczorek2007polynomial. This means that a sufficient condition for $\hat{\Sigma}_g(\beta_1^0)$ to be invertible is if it is invertible when $\beta_1^0 = 0$ (i.e. the determinant of $\hat{\Sigma}_g(0) \neq 0$), which should occur with probability one. We can invert $\hat{\Sigma}_g(\beta_1^0)$ using cofactor expansion. Thus, our Wald-like statistic can be rewritten as $$
$$ where $N_g(\beta_1^0) = T^{iv}_g(\beta_1^0)^T adj[\hat{\Sigma}_g(\beta_1^0)] T^{iv}_g(\beta_1^0)$, and $D_g(\beta_1^0) = det[\hat{\Sigma}_g(\beta_1^0)]$. Here, $adj[\bullet]$ and $det[\bullet]$ denote the adjugate and determinant of $\bullet$ respectively (these definitions apply throughout the paper). Note that $T^{iv}_g(\beta_1^0) \in \mathbb{R}^{k \times 1}[\beta_1^0]$, and $adj[\hat{\Sigma}_g(\beta_1^0)] \in \mathbb{R}^{k \times k}[\beta_1^0]$, so $N_g(\beta_1^0)$ and $D_g(\beta_1^0)$ are scalars in $\mathbb{R}[\beta_1^0]$.
These Wald-like statistics $t_g(\beta_1^0)$ differ from the statistics in the previous section because $t_g(\beta_1^0)$ is a rational function. Regardless, we can still find the intersection values algebraically by finding the real roots of $D_g(\beta_1^0)N_{\text{Id}}(\beta_1^0) - D_{\text{Id}}(\beta_1^0)N_g(\beta_1^0)$ = 0; we just need to ensure that $D_{\text{Id}}(\beta_1^0)$ and $D_g(\beta_1^0)$ are non-zero for these solutions. After finding the intersection values, we can use an algorithm similar to Algorithm (ref) to find the exact $p$-value curve from which we can easily extract the level $1-\alpha$ confidence interval(s); filter the output to only include rows for which the $p$-value is greater than $\alpha$, then the confidence interval(s) are the regions defined by $\beta_1^{\text{start}}$ and $\beta_1^{\text{end}}$.
Thus far we have only considered tests where the coefficient of interest $\beta_1^0$ is a scalar, but economists are often interested in constructing confidence regions for multiple coefficients and projecting this confidence region onto each dimension to create projected confidence intervals. We present two methods to create such intervals—fast grid search and optimization. Fast grid search works generally, but does not yield exact projected confidence intervals. We can use optimization to find the exact projected confidence intervals, but this requires assumptions that we discuss later.
Consider the overidentified instrumental variable statistic from Pouliot2024, but where $\beta_1^0 \equiv [\beta^0_{1,1}, \ \beta^0_{1,2}]^T$ $\in \mathbb{R}^{2 \times 1}$ and we annihilate all $g \mathbf{X}_1$ variation when estimating $\hat{\varepsilon}$ so that $\hat{\Sigma}_g$ is not a function of $\beta_1^0$.\footnote{This is the precise Wald-like randomization test in Pouliot2024.} In this case, the Wald-like statistic is $$
$$ where $$
$$ and $$ \hat{\Sigma}_g = \frac{1}{n}\sum_{i=1}^{n} (\mathcal{Q}_1 \mathbf{Z})_i (\mathcal{Q}_1 \mathbf{Z})_i^T \left(\mathcal{Q}_2 g \mathbf{Y}\right)_i^2. $$ Here, $\mathbf{X}_1 \in \mathbb{R}^{n \times 2}$; and $\mathbf{Y}$, $\mathbf{X}_2$, $\mathbf{Z}$, $\mathcal{Q}_1$, and $\mathcal{Q}_2$ are defined in Section~\ref{onesidedPouliot}. Similar to before, $\hat{\Sigma}_g$ can be rewritten as $$
$$ where $$
$$ Notice here that $\hat{\Sigma}_g$ is a real valued matrix that we can invert using cofactor expansion (assuming it is invertible). Then $t_g(\beta_1^0)$ can be rewritten $$
$$ where $A_g = g^T \mathcal{Q}_1 \mathbf{Z} \hat{\Sigma}_g^{-1} \mathbf{Z}^T \mathcal{Q}_1 g$ and $\Omega_g$ is the center matrix in the third equality. Note that $t_g(\beta_1^0)$ can define a conic section: it is a polynomial of degree two in $\beta_{1,1}^0$ and $\beta_{1,2}^0$. Like above, we want to consider the values of $[\beta^0_{1,1}, \ \beta^0_{1,2}]^T$ for which $t_{\text{Id}}(\beta_1^0) = t_g(\beta_1^0)$ because it is at these $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ values that the rank of $t_{\text{Id}}(\beta_1^0)$ changes relative to the distribution of $t_g(\beta_1^0)$, which changes the $p$-value. To find this set of $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ values, we adopt a similar procedure to that above: $$
$$ This is the matrix representation of a conic section (it can be degenerate such that the solution could be a point, two parallel lines, etc.). We can classify the general type of conic section (circle, ellipse, parabola, or hyperbola) by its discriminant $-4\Delta$, where $\Delta$ is the determinant of the upper-left $2\times2$ submatrix of $\Omega_{Id} - \Omega_{g}$. The $4$ is irrelevant in this context, so the discriminant is equivalent to $\mathbf{X}^T_{1,1} (A_Id - A_g) \left(\mathbf{X}_{1,2} \mathbf{X}^T_{1,1} - \mathbf{X}_{1,1} \mathbf{X}^T_{1,2}\right) (A_Id - A_g) \mathbf{X}_{1,2}$. The details for the classification can be found in Table \ref{conic}. This conic section partitions the $\mathbb{R}^2$ space so that depending on which side $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ is relative to the conic section, $t_{Id}(\beta_{1,1}^0, \beta_{1,2}^0)$ will be greater or smaller than $t_{g}(\beta_{1,1}^0, \beta_{1,2}^0)$.
The $p$-value for any given value of $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ is $$ p(\beta_{1,1}^0, \beta_{1,2}^0) \equiv \frac{1}{M} \left(1 + \sum_{g \in \mathbf{G} \setminus \text{Id}} \mathbbm{1}\left[ t_{\text{Id}}(\beta_{1,1}^0, \beta_{1,2}^0) \leq t_g(\beta_{1,1}^0, \beta_{1,2}^0) \right] \right). $$ The projected $p$-value function is equivalent to finding the max $p$-value over all other coefficients. For example, the projected $p$-value function for $\beta_{1,1}^0$ is $$ p(\beta_{1,1}^0) = \max_{\beta_{1,2}^0} p(\beta_{1,1}^0, \beta_{1,2}^0). $$
Unlike the scalar $\beta_1^0$ case, finding the exact confidence region for the vector $\beta_1^0$ case is difficult. It involves calculating the $p$-values in each region of the partitioned $\mathbb{R}^2$ space, then choosing the regions for which the $p$-value is greater than $\alpha$. As mentioned above, we present two methods.
The first method is fast grid search, but compared with traditional grid search methods, which would calculate the test and randomization statistics from start for each value of $[\beta_{1,1}^0, \beta_{1,2}^0]^T$, our method is an order faster because we only need to pay the fixed cost of calculating the coefficients of the second degree polynomial function once. After calculating these values, the marginal cost to calculate the test and randomization statistics for values of $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ is very small. Calculating the $p$-value for any $[\beta_{1,1}^0, \beta_{1,2}^0]^T$ amounts to checking the proportion of $t_g(\beta_{1,1}^0, \beta_{1,2}^0)$ values greater than or equal to $t_{\text{Id}}(\beta_{1,1}^0, \beta_{1,2}^0)$. We fully outline the method in Algorithm (ref). To extract the level $1-\alpha$ valid confidence region(s), consider the values of $\beta_1^0$ for which the $p$-value is greater than $\alpha$, in addition to adjacent values of $\beta_1^0$ where rejection occurred.
The second method involved optimization, but it requires that all the conic sections be ellipses. This seems to be a reasonable assumption in practice. While this method is easily extended to Exact Wald-tests of arbitrary dimension, we present this method for the two-dimensional case.
Note that when $-\Delta_g < 0$, then $f_g(\beta_{1,1}^0, \beta_{1,2}^0) =
(\Omega_{Id} - \Omega_g)
^T \leq 0 $ is a convex elliptical disk. We plot three such disks for visualization purposes.
In this toy example, we see 8 distinct regions (since for each ellipse, we can either be inside or outside of it). To construct an exact projected $p$-value curve, one method is create a axis-aligned bounding box for each region so that we can easy project each region onto the $\beta{1,1}^0$ and $\beta_{1,2}^0$ axes.
One could theoretically construct a bounding box for each region by solving a quadratically constrained problem. For instance, the optimization problem find the bounding box for region is $$
$$ Clearly, however, the quadratic constraints are not convex, so finding the maximum value of $\beta_{1,1}^0$ subject to the constraints will not be easy. To solve this problem, one could remove the constraint that $f_3(\beta_{1,1}^0, \beta_{1,2}^0) > 0$. Then, instead of finding the bounding box for region $x$, we could easily find the bounding box for region $x \cup y$. While this is a larger bounding box, this has no effect on the projected $p$-value curve because we overwrite the projected $p$-value defined by $x \cup y$ with the projected $p$-value defined by $y$.
If we were interested in the projected confidence intervals for a particular level of significance, notice that there is no need to compute the entire projected $p$-value curve for $\beta_{1,1}^0$ and $\beta_{1,2}^0$. Only the portion of the $p$-value curve for which the $p$-value is greater than our level of significance, as that is the interval for which we fail to reject the null hypothesis. Therefore, we can modify the optimization problem to be a mixed integer quadratically constrained problem.
$$
$$
Note that we added a constraint that the solution needs to be in a certain number of ellipses, that corresponds to our level of significance $\alpha$. These optimization programs can be solved using modern solvers.
The purpose of this paper to develop a general procedure to calculate the exact confidence intervals for permutation tests. To that end, we demonstrate that this method also extends to other tests in this section.
d2024robust also develop a heteroskedasticity-robust, Wald-like test. Their test, however, only handles discrete nuisance covariates, so consider the setup in Section (ref) where $\mathbf{X}_2$ is composed of discrete covariates. We are primarily focused on the structure of the test and randomization statistic, so read d2024robust for details about the inner workings of the test. The randomization statistic that they define is $$ \mathcal{W}^{\pi} = g(W, (y - X\beta_0)_{\pi}), $$ where $$ g(W, v) = v' \Tilde{X}(\Tilde{X}' \Sigma(v) \Tilde{X})^{-1} \Tilde{X}' v; $$ $W = [X, Z]$; $\Tilde{X}$ is a matrix derived from $X$; $v \in \mathbb{R}^{n \times 1}$; $\Sigma(v) \in \mathbb{R}^{n \times n}$ is a diagonal matrix with $(i,i)$ element equal to $(Dv)_i^2$; $D \in \mathbb{R}^{n \times n}$; and $\pi \in \mathbb{R}^{n \times n}$ is a permutation matrix. Like with the tests in Pouliot2024, the test statistic corresponds to the case when $\pi = \text{Id}$. Notice that $$
$$ From this perspective, it is easy to see that the structure of the test statistic is similar to the over-identified instrumental-variable test statistic in Section \ref{waldPouliot} when $\beta_0$ is a scalar: $D$ and $\Tilde{X}$ are analogous to $\mathcal{Q}_3$ and $\mathcal{Q}_1 \mathbf{Z}$ respectively. We can follow the procedure outlined in Section \ref{waldPouliot} or \ref{jointValidPouliot} depending on whether $\beta^0$ is a scalar or vector to find the exact or approximate $p$-value curve and level $1-\alpha$ confidence region respectively.
Another test that we touched upon in Section (ref) is that of diciccio:2017aa. diciccio:2017aa develop a partial correlation test for specific components. Their test is not exact in finite samples, but is asymptotically valid. d2024robust compare their method against the first randomization test in Section 3 of diciccio:2017aa. That randomization statistic is $$
$$ where $\hat{\beta}_{\pi}$ is the regression coefficients obtained from regressing $\pi Y$ on $X$; $\hat{\Sigma}^T_{\mathbf{XX}} = \frac{1}{n}\sum_i^n X^T_i X_i$ and $\hat{\Omega}_{\pi}(\beta^0) = \frac{1}{n} \sum_i^n ((\pi Y)_i - X^T_i \beta^0)^2 X_i X^T_i$. Like with the other randomization tests, the test statistic corresponds to $\pi = Id$.
For simplicity, consider the case where $\beta^0 \equiv [\beta_1^0, \beta_2^0 ]^T \in \mathbb{R}^2$. Notice that $\hat{\Omega}(\beta^0)$ is similar in form to $\hat{\Sigma}_g(\beta_1^0)$ in the Wald statistic discussed in Section (ref). Then $\hat{\Omega}_{\pi}(\beta^0)$ can be rewritten as $$
$$ where $$
$$ and $\Sigma_{vw} = diag[\mathbf{X}_v \circ \mathbf{X}_w]$. Note that $\hat{\Omega}(\beta^0) \in \mathbb{R}^{2 \times 2}[\beta_1^0]$ and $(\hat{\beta} - \beta^0) \in \mathbb{R}^{2 \times 1}[\beta_1^0]$, so $U_{\pi}(\mathbf{X}, \mathbf{Y}) \in \mathbb{R}[\beta_1^0]$. Since $U_{\pi}(\mathbf{X}, \mathbf{Y})$ is just a polynomial, we can apply the same principles in Section \ref{jointValidPouliot} to construct jointly valid confidence regions. We will need to use grid search here, but we stress that our method is an order of magnitude faster than traditional grid search because we only calculate the coefficients of the polynomial for every permutation once rather than for each hypothesized value of $\beta^0$. Consequently, calculating the $p$-values and subsequent level $1-\alpha$ confidence interval for some region of $\beta^0$ will be much faster than traditional grid search.
The traditional method to produce confidence intervals in randomization-based inference in grid search, which is computationally expensive and leads to conservative confidence intervals as evidenced in d2024robust. We show that many test and randomization statistics, including those where the variance estimator depends on the null hypothesis, are rational functions (if the variance estimator does not depend on the null hypothesis, the statistics reduce to a polynomial). We propose a method to use this structure to produce exact confidence intervals for a scalar $\beta_1^0$ for a broad class of test statistics. For cases where the null hypothesis $\beta_1^0$ is a vector, we use grid search; however, we show that by storing the coefficients of the function, our method is an order of magnitude faster than traditional grid search methods which require computing the test and randomization statistics from scratch for each null hypothesis value.