The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
77,827 characters
When Does Inexact Matching Ensure Balance and Inference without Adjustment?
\maketitle
\begin{abstract}
One-to-one matching without replacement is a classical approach to constructing comparable treated and control samples in the design of observational studies. It pairs each treated unit with a distinct control while minimizing a covariate distance objective. With continuous covariates, the matched pairs generally remain inexact, which contributes to bias in downstream analysis. Its key theoretical properties, such as the resulting imbalance between matched pairs and when it is negligible to support valid inference, remain unclear.
In this paper, we analyze one-to-one matching based on $d$-dimensional, continuous covariates with a quadratic covariate-distance objective.
First, we find that when $d\leq 3$, under standard conditions on the propensity score ensuring abundant control samples near each treated sample, the imbalance (difference between within-group averages) is root-$n$ negligible uniformly over the family of smooth functions with a common first- and second-order derivative bound. However, such balance is subject to a dimension restriction, as we construct examples in which the imbalance is root-$n$ non-negligible when $d=4$ and dominates root-$n$ rate when $d>4$.
Second, we show that when $d\leq 3$, the matched design allows valid Wald-type and bootstrap inference for the average treatment effect on the treated, distributional treatment effects, and quantile treatment effects. Thus, the same outcome-blind matched design supports various downstream inferences without having to tailor the design to the targets. Finally, paired randomization inference based on the matched design is asymptotically valid in the super-population sense for $d\leq 3$ but can fail when $d=4$. We corroborate the theoretical results with numerical experiments.
\end{abstract}
\section{Introduction}
Matching is widely used in the observational studies to construct comparable treated and control groups~\citep{rubin1973matching,rosenbaum1985constructing,rosenbaum2020design}.
One canonical form is one-to-one matching without replacement, which pairs each treated unit with a distinct control using their pre-treatment covariates~\citep{rosenbaum1989optimal} while minimizing certain objective for the distance between matched pairs.
In this process, the pairs are formed before the outcomes are examined, thereby separating the design and analysis of observational studies and promoting transparency and objectivity.
One criterion of a good match is the matched pairs be close enough---and ideally identical---in the pre-treatment covariates to help reduce bias from differences in observed covariates~\citep{rubin1973matching}. However, with continuous covariates, the matching generally remains inexact, and the differences in the covariates can leave bias in the subsequent analysis.
A second criterion is balance between the treated and matched control groups \citep{rosenbaum1989optimal}. Even when individual pairs are inexact, their differences may cancel across pairs, making the overall bias relevant in downstream inference smaller than the individual discrepancies suggest.
The key question is thus whether such bias is negligible relative to the sampling uncertainty so valid inference is possible.
Establishing this connection is challenging: because assigning a control to a treated unit can change the choices available to the others, the matching assignments need to be analyzed jointly across all samples. Although matching is intuitive, transparent, and simple to describe, the implications of inexact matching for the balance of the selected groups and the validity of downstream inference remain unclear.
Existing theory provides important results for matching estimators, including matching without replacement and inference after matching~\citep{AI2012,AS2022}. There, a central requirement for proving inferential guarantees is that the approximation error (pairwise distance) be sufficiently small. This requirement typically requires a large amount of control samples, yet does not necessarily follow from that alone: \cite{Savje2022} shows that the shortage of controls in a local region of the covariate space can result in inconsistency.
Even when the controls are abundant locally in the entire sample space, the magnitude of the matching bias still matters.
For optimal pair matching based on ordinary Mahalanobis distance, \cite{GR2023} establish the validity of an unadjusted paired randomization test with one continuous covariate and conjecture an extension to a dimension of three.
In this paper, we study optimal one-to-one matching with a quadratic covariate-distance objective commonly used in the literature~\citep{rubin1980bias,rosenbaum1985constructing,bind2019bridging,zhang2023matching}. We work under an i.i.d.~super-population model where the covariates are $d$-dimensional, bounded, continuously valued, and the dimension $d$ is viewed as fixed as $n\to \infty$.
We assume a favorable setting for matching in which the covariates admit a density with respect to the Lebesgue measure on $[0,1]^d$ that is bounded above and away from zero. The propensity score is smooth, bounded away from zero and strictly below $1/2$ to ensure that control samples are abundant around every treated sample, thereby separating the key analysis from the complications pointed out by~\cite{Savje2022}.
In this regime, the treated and control sample sizes grow at the same order, so we do not require the control reservoir to grow extremely large~\citep{AI2012,AS2022}.
\subsection{Preview of results}
Our main result concerns the imbalance of the matched sample over smooth functions. Let $\{T_i\}_{i=1}^{n_1}$ and $\{C_j\}_{j=1}^{n_0}$ be the treated and control covariates, respectively, both generated under a super-population model (Section~\ref{subsec:setup}), and let the learned matching be represented by an injection $\hat\sigma_n \colon \{1,\dots,n_1\}\to \{1,\dots,n_0\}$ which matches treated sample $i$ to control sample $\hat\sigma_n(i)$, $i=1,\dots,n_1$. In Section~\ref{sec:balance}, we show that, when $d\leq 3$, the minimizer $\hat\sigma_n$ of the quadratic matching objective obeys
\@\label{eq:intro_balance}
\sup_{g\in \cF_L} \bigg| \frac{1}{n_1}\sum_{i=1}^{n_1} \big\{g(T_i) - g(C_{\hat\sigma_n(i)}) \big\} \bigg| = o_P(n^{-1/2}),
\@
where $\cF_L$ is the class of twice-continuously-differentiable functions with uniformly bounded first- and second-order derivatives. Notably, the functions in this class are not fitted or included as constraints in the matching problem: the procedure simply minimizes the total quadratic distance. The result also allows the matrix defining the distance (such as that in the Mahalanobis distance) to be estimated in root-$n$ rate, and allows $\hat\sigma_n$ to be an approximate solution with a small optimization gap.
The uniform balance~\eqref{eq:intro_balance} for $d\leq 3$ appears to respect a sharp dimension restriction. We further construct an example in which the imbalance in one function $g$ is non-negligible at the root-$n$ scale for $d=4$, and dominates the root-$n$ scale when $d>4$, even though all the stated favorable conditions hold.
This result partially echoes the conjecture in~\cite{GR2023} albeit using the quadratic objective.
This balance result provides a common basis for various downstream inferences based on one matched design.
We study a variety of population-level inference tasks in Section~\ref{sec:wald_inference}.
Under suitable outcome conditions, for $d\leq 3$, we establish asymptotic normality and consistent variance estimation for the population average treatment effect on the treated (ATT), as well as valid pairwise bootstrap inference without recomputing the matching. To demonstrate the breadth of inferential guarantees, we show the matched design can analogously support pairwise bootstrap inference for joint confidence bands for the counterfactual distribution functions, distributional treatment effects, and quantile treatment effects.
In Section~\ref{sec:rand_inference}, we turn to the randomization-based inference commonly used after a matched design.
We find that for $d\leq 3$, the paired difference-in-means randomization test is asymptotically valid under the sharp and weak null in the super-population sampling sense.
A dimension boundary is present, as we show that for $d=4$, the randomization test for the weak null can violate the type-I error at any nominal level $1-\alpha$.
Together, these results justify the use of one outcome-blind matched design across several analyses for low-dimensional problems. Under the stated conditions, the matching and the inferential procedures studied in this paper do not require outcome regression or propensity-score estimation. While each inferential target requires distinct comparisons of the matched pairs (such as the paired averages in mean effect estimation and the empirical CDFs for distributional effect estimation), the uniform balance alone automatically controls the discrepancies in the variety of outcome-related functions relevant in each analysis, without knowing these functions in advance or involving them in constructing the matching. It remains in the design stage of the observational study but retains strong inferential guarantees when the feature dimension is small.
Finally, we corroborate our theory with numerical demonstrations in Section~\ref{sec:simu}. In simulations across ATT, distributional and quantile treatment effects, and randomization tests, we observe robust coverage of confidence intervals/bands when $d\leq 3$, systematic undercoverage at $d=4$, and deteriorating coverage for $d>4$. We include a detailed discussion on the related literature in Section~\ref{sec:literature}.
\subsection{Problem setup and notations}
\label{subsec:setup}
We assume a super-population setting where the full sample $\{(X_i,Z_i,Y_i(1),Y_i(0))\}_{i=1}^n$ are i.i.d.~from a common distribution $\PP$, and we observe the i.i.d.~triplets $\{(X_i,Y_i,Z_i)\}_{i=1}^n$, with binary treatment $Z_i\in \{0,1\}$, pre-treatment covariates $X_i\in [0,1]^d$, and outcome $Y_i\in \RR$.
Throughout the paper, we assume standard SUTVA and unconfoundedness, so that $\{Y(1),Y(0)\}\indep Z\given X$ under $\PP$, and $Y_i=Y_i(Z_i)$.
The propensity score is denoted as $e(x)=\PP(Z_i=1\given X_i=x)$.
In addition, we make the following assumptions on the covariate density and strictly upper bounded propensity scores to ensure a clean and favorable setting for matching.
\begin{assumption}\label{assump:obs-sample}
Assume the i.i.d.~super-population model with SUTVA and unconfoundedness, so that $\{Y(1),Y(0)\}\indep Z\given X$ under $\PP$, and $Y_i=Y_i(Z_i)$.
In addition, suppose $X$ has marginal density $f$ with respect to the Lebesgue measure on $\Omega=[0,1]^d$. Suppose $f(\cdot)$ and $e(\cdot)$ have $C^1$ extensions to
a neighborhood of $\Omega$. For fixed
constants $\underline f,\overline f,\eta_0,\eta_1$, we have
$0<\underline f\le f(x)\le\overline f<\infty$ and
$\eta_0\le e(x)\le\tfrac12-\eta_1$,
for all $x\in\Omega$.
\end{assumption}
The condition $\eta_0\le e(x)\le\tfrac12-\eta_1$ is also assumed in \cite{GR2023}, which ensures sufficiently abundant control samples in every local region to rule out the complication mentioned in~\cite{Savje2022}, so that it is eventually possible to find a control for
each treated observation.
Under this assumption, the event $\cI_n=\{2\le n_1\le n_0\}$ happens with probability tending to one.
\begin{remark}
The super-population model used in this paper follows that of~\cite{GR2023} but differs from some common alternatives in the literature, such as the design-based inference where $\{(X_i,Y_i(1),Y_i(0))\}_{i=1}^n$ are viewed as fixed and the sole source of randomness is the treatment assignment. In that regime, the randomness in the matching design arises from the treatment assignment only, and the fixed covariate values makes the analysis even more challenging. In contrast, the super-population model enables the analysis of large-sample behavior via concentration inequalities under transparent conditions. Our setting ensures $n_0$ and $n_1$ are of the same order, in contrast to the requirements in~\cite{AI2012} that $n_1^r/n_0=O(1)$ for some $r>d$ so the control reservoir grows much faster than the $n_1$.
\end{remark}
For notational convenience, we re-index the treated and control covariates as $\{X_i\}_{Z_i=1} =\{T_i\}_{i=1}^{n_1}$ and $\{X_i\}_{Z_i=0} = \{C_j\}_{j=1}^{n_0}$, where $n_1=\sum_{i=1}^n Z_i$ and $n_0=n-n_1$ are the treated and control group sizes.
We focus on the one-to-one matching problem with a quadratic objective.
Formally, consider an injection $\sigma\colon \{1,\dots,n_1\}\to \{1,\dots,n_0\}$, and let $\cA_{n_1,n_0}$ be the collection of all such injections. For a fixed symmetric positive-definite matrix $A_0\in \RR^{d\times d}$, we define the matching objective
\@\label{eq:def_matching_obj}
{Q}_n(\sigma;A_0) = \frac{1}{2n_1}\sum_{i=1}^{n_1} (T_i-C_{\sigma(i)})^\top A_0 (T_i-C_{\sigma(i)}),\quad \sigma\in \cA_{n_1,n_0},
\@
and the global optimal solution
\$
\sigma_n^\star \in \argmin_{\sigma\in \cA_{n_1,n_0}}Q_n(\sigma;A_0).
\$
Later on, we shall also extend the results so $A_0$ can be replaced by random, data-dependent matrices such as the inverse covariance matrix.
Write $Q_n^\star = \min_\sigma Q_n(\sigma;A_0)$ as the minimized matching objective.
In practice one may settle at an approximate solution $\hat\sigma_n$ with an objective gap $\Delta_n = Q_n(\hat\sigma_n;A_0) - Q_n^\star \geq 0$.
\begin{remark}
When $A_0=I$, the objective~\eqref{eq:def_matching_obj} is the total squared Euclidean distance. When $A_0$ is the (empirical) inverse covariance matrix, the objective is the total squared Mahalanobis distance. This differs from the Mahalanobis norm objective in \cite{GR2023}, but enables principled theoretical analysis and is broadly used in practice~\citep{bind2019bridging,zhang2023matching}.
\end{remark}
\paragraph{Notations.} We close the section by stating the notations used throughout this paper.
For a vector $v$, we denote it Euclidean norm by $\|v\|$. For a matrix
$M$, let
$\|M\|_{\mathrm{op}}=\sup_{\|v\|=1}\|Mv\|$
be its operator norm; when $M$ is symmetric,
$\lambda_{\min}(M)$ and $\lambda_{\max}(M)$ denote its smallest
and largest eigenvalues. We write $I_d$ for the $d\times d$
identity matrix. For a scalar function $g$ on a set $S$, let
$\|g\|_S=\sup_{x\in S}|g(x)|$, which is also written as $\|g\|_\infty$ when no confusion arises. The gradient and Hessian of $g$ are denoted by
$\nabla g$ and $\nabla^2 g$, respectively. For a square-integrable
real random variable $W$, we denote its $L_2$-norm as
$\|W\|_2=\{\EE(W^2)\}^{1/2}$.
We use $\ind\{\cdot\}$ for an indicator and
$x_+=\max\{x,0\}$ for the positive part.
For positive deterministic sequences $a_n,b_n$, write
$a_n\lesssim b_n$, or equivalently $a_n=O(b_n)$, if
$a_n\le Cb_n$ for some constant $C<\infty$ and all sufficiently
large $n$. Write $a_n\asymp b_n$ if both $a_n\lesssim b_n$ and
$b_n\lesssim a_n$, and $a_n=o(b_n)$ if $a_n/b_n\to0$.
For a positive sequence $a_n$, the notation
$U_n=O_P(a_n)$ means that $\|U_n\|/a_n$ is bounded in probability,
whereas $U_n=o_P(a_n)$ means that $\|U_n\|/a_n\to 0$ in probability;
the norm is chosen according to whether $U_n$ is a scalar,
vector, matrix, or function. We use $\to_P$ for convergence
in probability and $\Rightarrow$ for weak
convergence. Unless stated otherwise, all limits are taken as
$n\to\infty$ with $d$ fixed. Constants denoted by $c$ or $C$
may change from line to line and do not depend on the sample sizes;
their dependence on other quantities is specified where needed.
\section{Uniform balance from optimal quadratic matching}
\label{sec:balance}
We are interested in the balance between the matched groups.
For a function $g\colon \cX\to \RR$, define
\@\label{eq:def_imbalance_g}
B_{n,\sigma}(g) = \frac{1}{n_1}\sum_{i=1}^{n_1} \big\{ g(T_i) - g(C_{\sigma(i)})\big\}.
\@
For a fixed $L\in (0,\infty)$, we define $\cF_L$ as the set of functions with a $C^2$ extension on $\Omega$ and whose function values and all first-order and second-order partial derivatives are bounded in absolute value by $L$ on $\Omega$.
\subsection{Uniform balance for \texorpdfstring{$d\leq 3$}{d <=3}}
\label{subsec:uniform_balance}
Our first result is the uniform balance ensured by the matched design in Theorem~\ref{thm:observational}, whose proof is in Appendix~\ref{app:sec_proof_uniform_balance}.
For completeness, we define the corresponding optimization gap and balance measure to zero on $\cI_n^c$, which does not affect the high-probability result.
\begin{theorem}[Uniform balance]
\label{thm:observational}
Suppose $d\le3$, and Assumption~\ref{assump:obs-sample} holds.
Let $\hat{A}$ be a data-dependent symmetric matrix such that $\norm{\hat{A}-A_0}_{\textnormal{op}} = O_P(1/\sqrt{n})$ for some fixed symmetric positive-definite matrix $A_0$.
On the event $\cI_n$, let
$\hat\sigma_n\in\mathcal A_{n_1,n_0}$ be a measurable injection with optimization gap
$\xi_n = Q_n(\hat\sigma_n;\hat{A})
-\min_{\sigma\in\mathcal A_{n_1,n_0}}Q_n(\sigma;\hat{A}) = o_P(1/n)$.
Then for every
fixed $L\in(0,\infty)$, as $n\to \infty$,
\@\label{eq:observational-balance}
\sqrt n\sup_{g\in\mathcal F_L}
|B_{n,\hat\sigma_n}(g)| = o_P(1).
\@
\end{theorem}
The convergence in probability stated in~\eqref{eq:observational-balance} concerns the random imbalance as a function of the covariates. Such randomness arises from the super-population model, where i.i.d.~sampling of the covariates provides the probabilistic basis and needed structure for controlling the imbalance with high probability.
The condition $\|\hat{A}-A_0\|_{\textnormal{op}}$ allows the use of Mahalanobis distance where $\hat{A}= \hat\Sigma^{-1}$ for the empirical covariance matrix $\hat\Sigma$ which converges to the population covariance matrix in the required rate if tha latter has lower-bounded eigenvalues.
Theorem~\ref{thm:observational} implies that the matched control sample approximates the treated group in the sense of empirical measure. Namely, on $\mathcal I_n=\{2\le n_1\le n_0\}$, consider the empirical measures $\hat{P}_1 = n_1^{-1}\sum_{i=1}^{n_1} \delta_{T_i}$ and $\hat{Q}_0 = n_1^{-1} \sum_{i=1}^{n_1} \delta_{C_{\hat\sigma_n(i)}}$. It states that $\sqrt{n}\|\hat{P}_1-\hat{Q}_0\|_{\cF_L} = o_P(1)$, where we define $\|\mu\|_{\cF_L} = \sup_{g\in \cF_L} |\int g \ud \mu|$.
\begin{remark}[Fixed-function moment bounds]
The imbalance control in~\eqref{eq:observational-balance} is stated in probability. The proof of Theorem~\ref{thm:observational} (especially the intermediate result Theorem~\ref{thm:balance}) also yields stronger moment bounds for any fixed function $g\in \cF_L$
and exact optimizer under a fixed symmetric positive definite matrix $A_0$.
For all sufficiently large $n$, recalling that $\sigma_n^\star$ is the exact optimizer, we have
\[
\left|\EE [B_{n,\sigma_n^\star}(g)]\right|
\le C_g\Gamma_n,
\qquad
\EE\!\left[B_{n,\sigma_n^\star}(g)^2\right]
\le C_g\Big(\Gamma_n^2+\frac{\Gamma_n}{\sqrt n}\Big),
\qquad
\Gamma_n\asymp n^{-2/d}(\log n)^{1+2/d},
\]
where the constant $C_g$ does not depend on $n$.
Thus, for $d\leq 3$, the imbalance $B_{n,\hat\sigma_n}(g)$ is root-$n$ negligible in the $L^2$-norm. See Proposition~\ref{prop:observational-moments} for a formal statement.
\end{remark}
We now discuss our analysis techniques. A common way in the literature to control the approximation error left by matching is
to bound each pair's contribution in absolute value. For a smooth
function $g$,
\@\label{eq:first_order}
|B_{n,\sigma}(g)|
\le
\frac1{n_1}\sum_{i=1}^{n_1}
|g(T_i)-g(C_{\sigma(i)})|
\le
\bigg(\sup_{x\in\Omega}\|\nabla g(x)\|\bigg)
\frac1{n_1}\sum_{i=1}^{n_1}\|T_i-C_{\sigma(i)}\|.
\@
This bound gives root-$n_1$ negligible imbalance whenever the average
matching distance is $o_P(n_1^{-1/2})$, as required by the
approximation conditions in
\citet[Theorem~1]{AI2012} and \citet[Assumption~3]{AS2022}. However, this requirement is difficult to fulfill, which leads to requiring a large control pool that grows much faster than $n_1$.
The key insight of our analysis is that the first-order expansion may neglect the cancellation between positive
and negative pairwise differences.
To see how this cancellation may happen, by Taylor expansion,
\[
B_{n,\sigma}(g)
=
\frac1{n_1}\sum_{i=1}^{n_1}
\nabla g(T_i)^\top\{T_i-C_{\sigma(i)}\}
+R_{n,\sigma}(g),
\qquad
|R_{n,\sigma}(g)|
\le
\frac{dL}{2n_1}\sum_{i=1}^{n_1}
\|T_i-C_{\sigma(i)}\|^2.
\]
The remainder is controlled by the quadratic matching objective.
The leading sum, however, keeps the signs of the differences, and bounding each term by the magnitude in~\eqref{eq:first_order} could be loose. The main step in our analysis
is to control this signed sum for the globally optimal assignment
under the stated sampling assumptions.
For a fixed function $g$ and a fixed matrix $A_0$, our proof (Appendix~\ref{app:sec_proof_uniform_balance}) studies
the optimal value of the tilted objective
$Q_n(\sigma;A_0)+tB_{n,\sigma}(g)$ for $t\in \RR$. Its derivative at $t=0$ equals the imbalance of
the original optimizer $\sigma_n^\star$ almost surely.
We use integration by parts, as well as the Efron--Stein inequality to obtain bounds on the optimal values and translate them to bounds on the derivative. The E--S inequality is applicable because we prove a stability property of the optimal value under one-observation replacement (Lemmas~\ref{lem:stab_opt_value} and~\ref{lem:value-bounds}).
We then use a standard chaining technique to extend to uniform balance over $\cF_L$.
This distinction also connects to \citet[Proposition~1]{GR2023}, who uses the conditional matching bias to characterize the validity of paired randomization test.
Our result establishes sufficient conditions for uniformly negligible
smooth-function imbalance (bias) under the quadratic objective.
\subsection{The dimension boundary: counterexample for \texorpdfstring{$d\geq 4$}{d >= 4}}
\label{subsec:dim}
The uniform balance result holds for $d\leq 3$. A natural question is whether this is the sharp dimension boundary.
We now show that with $d\geq 4$, even for the exact optimizer $\sigma_n^\star$ for a fixed matching objective with $\hat A=A_0=I_d$ and $\xi_n=0$ can fail to yield $\sqrt{n}$-negligible imbalance. The construction of the example and the proof of Theorem~\ref{thm:dimension-obs} are in Appendix~\ref{app:sec_proof_boundary}.
\begin{theorem}\label{thm:dimension-obs}
For every fixed integer $d\ge4$, there exist fixed $C^\infty$ functions $f,e$ satisfying Assumption~\ref{assump:obs-sample} and a fixed $C^\infty$ function $g$ on a neighborhood of $\Omega$, with $0\le g\le1$ and $\sup_\Omega\|\nabla g\|\le1$, such that the exact minimizer $\sigma_n^\star$ under $A_0=I_d$ has the following properties. If $d=4$, there are constants $b_*,p_*>0$ such that
\begin{equation}\label{sharp:eq:iid-prob}
\liminf_{n\to\infty}\PP\{\sqrt{n_1}\,B_{n,\sigma_n^\star}(g)\ge b_*\}\ge p_*.
\end{equation}
If $d>4$, then $\sqrt n B_{n,\sigma_n^\star}(g)\to+\infty$ in probability. The imbalance is set to be zero on $\mathcal I_n^c$.
\end{theorem}
In the construction, the expected imbalance is of order $n^{-2/d}$, while its variance is of order at most $n^{-1}$. The former coincides with the sampling-error scale at $d=4$ and decays more slowly when $d>4$. For matching with replacement, \cite{AI2006} pointed out the bias scale of $n^{-2/d}$, which coincides with ours in the matching without replacement setting.
Finally, the single-$g$ bound clarifies that the failure of uniform balance is due to inexact matching, instead of the intrinsic sampling uncertainty due to taking the uniform over the function class for $d=4$ (which happens to be $\sqrt{n}$).
\section{Wald-type and bootstrap inference from one matched design}
\label{sec:wald_inference}
A useful implication of the uniform balance in Theorem~\ref{thm:observational} is that one matched design provides negligible bias in a wide range of functions, so that the impact of inexact matching on various downstream inference is negligible.
In this section, we show how the same matched design allows convenient Wald-type inference and Bootstrap approximation for various targets, including the average treatment effect on the treated (ATT) (Section~\ref{subsec:wald_att}), distribution effects (Section~\ref{subsec:wald_cdf}), and quantile treatment effects (Section~\ref{subsec:wald_quantile}).
Throughout this section, for notational convenience, we denote the treated outcomes as $\{Y_{1,i}\}_{i=1}^{n_1}$ and the control outcomes as $\{Y_{0,j}\}_{j=1}^{n_0}$. Then, a matched design $\hat\sigma_n$ matches each treated sample $i$ to the control sample $\hat\sigma_n(i)$.
We also consider a bootstrap approximation to the Wald-type inference.
\paragraph{Bootstrapping the matched pairs.} We consider the bootstrap inference when resampling the matched pairs. Given the observed data, the original matching $\hat\sigma_n$, and bootstrap size $B\in \NN^+$, for each $b=1,\dots,B$, we draw $I_1^{*(b)},\ldots,I_{n_1}^{*(b)}$ independently and uniformly from $\{1,\ldots,{n_1}\}$ with replacement.
The $b$-th bootstrap matched pairs are
$(Y_{1,I_i^{*(b)}}, Y_{0, \hat\sigma_n(I_i^{*(b)})} )$, $i=1,\dots, n_1$.
Note that we do not re-compute the matched design.
\vspace{1em}
Throughout, we assume the conditions in Theorem~\ref{thm:observational}, which we state for easier reference. The only additional requirement is that $\hat{A}$ must also use the outcome-blind information as the matching.
\begin{assumption}\label{assump:match}
Suppose $d\leq 3$.
Let $\cG_n = \sigma((X_i,Z_i)_{i=1}^n)$ be the $\sigma$-field used in constructing the matching.
Let $\hat{A}$ be a data-dependent symmetric matrix that is $\cG_n$-measurable such that $\norm{\hat{A}-A_0}_{\textnormal{op}} = O_P(1/\sqrt{n})$ for some fixed symmetric positive-definite matrix $A_0$.
On the event $\cI_n=\{2\leq n_1\leq n_0\}$, let
$\hat\sigma_n\in\mathcal A_{n_1,n_0}$ be a $\cG_n$-measurable injection with optimization gap
$\xi_n = Q_n(\hat\sigma_n;\hat{A})
-\min_{\sigma\in\mathcal A_{n_1,n_0}}Q_n(\sigma;\hat{A}) = o_P(1/n)$.
\end{assumption}
Throughout, the probability and asymptotics concern the super-population sampling uncertainty, so the covariates, $\hat\sigma_n$, and outcomes are random, and the extra bootstrap sampling randomness when applicable.
\subsection{Average treatment effect on the treated}
\label{subsec:wald_att}
The population ATT is defined as $\theta_{{\textnormal{ATT}}} := \EE[Y(1)-Y(0)\given Z=1]$. Given a matched design $\hat\sigma_n$, we estimate with the difference-in-means estimator, with the Wald-type variance estimator
\@\label{eq:def_Ri_theta_att}
\hat\theta_{{\textnormal{ATT}}} = \frac{1}{n_1}\sum_{i=1}^{n_1} R_i,\quad \text{where}\quad R_i := Y_{1i} - Y_{0\hat\sigma_n(i)},\quad \hat{V}_{{\textnormal{ATT}}} := \frac{1}{n_1-1} \sum_{i=1}^{n_1} (R_i-\hat\theta_{{\textnormal{ATT}}})^2.
\@
For bootstrap inference, we consider the studentized t-statistics. For $b=1,\dots,B_n$ and some fixed $B_n\in \NN^+$,
\@\label{eq:def_bootstrap_att_tstat}
T^{*(b)} = \frac{\hat\theta_{{\textnormal{ATT}}}^{*(b)} - \hat\theta_{\textnormal{ATT}}}{\sqrt{\hat{V}_{\textnormal{ATT}}^{*(b)}/n_1}},\quad \hat\theta^{*(b)}_{\textnormal{ATT}} = \frac{1}{n_1}\sum_{i=1}^{n_1} R_{I_i^{*(b)}} ,\quad
\hat{V}_{\textnormal{ATT}}^{*(b)} = \frac{1}{n_1-1}\sum_{i=1}^{n_1} (R_{I_i^{*(b)}} - \hat\theta_{\textnormal{ATT}}^{*(b)} )^2.
\@
We make standard moment and smoothness conditions on the outcomes. Let $\mu_z(x)=\EE[Y(z)\given X=x]$.
We write the conditional treatment effect (CATE) function as $\tau(x) = \mu_1(x)-\mu_0(x)$, and the conditional variance function $\sigma_z^2(x):=\EE[(Y(z)-\mu_z(x))^2\given x=x]$. Let $\Var_z(\cdot)$ denote the variance conditional on $Z=z$.
\begin{assumption}\label{assump:outcome}
The function $\mu_0$ has a $C^2$ extension near $\Omega$, $\mu_1$ is bounded and measurable, and $\sigma_0^2$ is continuous.
In addition, $\EE[(Y(z)-\mu_z(x))^4\given X=x]\leq M$ for some a common constant $M>0$ for all $x$ and $z\in \{0,1\}$. Define $v_\tau=\Var_1\{\tau(X)\}$, $v_\epsilon=\E_1\{\sigma_1^2(X)+\sigma_0^2(X)\}$,
and suppose $V_{\textnormal{ATT}}:=v_\tau+v_\epsilon>0$.
\end{assumption}
We show that when $d\leq 3$, the matched-pair difference-in-means estimator is asymptotically normal with variance $V_{\textnormal{ATT}}$ defined above, and permits valid Wald-type and bootstrap inference.
\begin{theorem}[Population ATT inference]\label{thm:ATT}\label{thm:bootstrap}
Under Assumptions~\ref{assump:obs-sample},~\ref{assump:match}, and~\ref{assump:outcome}, as $n\to \infty$,
\begin{equation}\label{eq:ATT-main}
\sqrt {n_1}(\hat\theta_{\textnormal{ATT}}-\theta_{\textnormal{ATT}})\Rightarrow N(0,V_{\textnormal{ATT}}),
\qquad \hat V_{{\textnormal{ATT}}}\to V_{{\textnormal{ATT}}}\quad\text{in probability}.
\end{equation}
Hence the Wald-CI $\hat\theta_{\textnormal{ATT}}\pm z_{1-\alpha/2}\sqrt{\hat V_{\textnormal{ATT}}/{n_1}}$ has limiting coverage of $1-\alpha$, where $z_{1-\alpha/2}$ is the $1-\alpha/2$-th quantile of a standard normal distribution. Furthermore, the matched-pair bootstrap satisfies
\begin{equation}\label{eq:bootstrap}
\sup_{x\in \RR}\left|\PP^*\{
T^{*(1)}\leq x \}-\Phi(x)
\right|\to 0
\quad\text{in probability}.
\end{equation}
For bootstrap, let $B=B_n$ be a
deterministic sequence of positive integers satisfying
$B_n\to\infty$ as $n\to\infty$.
For the t-statistics in~\eqref{eq:def_bootstrap_att_tstat}, define $c_{q}^*$ as the $q$-th empirical quantile of $\{T^{*(b)}\}_{b=1}^B$. Then the bootstrap confidence interval $\mathrm{CI}_{n,1-\alpha}^{\mathrm{boot}} = [\hat\theta_{{\textnormal{ATT}}} - c_{1-\alpha/2}^* \cdot(\hat{V}_{\textnormal{ATT}}/n_1)^{1/2}, \hat\theta_{{\textnormal{ATT}}} - c_{\alpha/2}^* \cdot (\hat{V}_{\textnormal{ATT}}/n_1)^{1/2}]$ has limiting coverage of $1-\alpha$.
\end{theorem}
The proof of Theorem~\ref{thm:ATT} is in Appendix~\ref{app:subsec_att_proof}.
Here we outline the intuitions for the root-$n$ Wald-type inference. Let $\tau(x)=\EE[Y(1)\given X=x]-\EE[Y(0)\given X=x]$ be the conditional average treatment effect. Then
\$
\sqrt{n_1}(\hat\theta_{{\textnormal{ATT}}}-\theta_{\textnormal{ATT}}) = \frac{1}{\sqrt{n_1}} \sum_{i=1}^{n_1} \{\tau(T_i) - \theta_{\textnormal{ATT}}\}
+ \sqrt{n_1} B_{n,\hat\sigma_n}(\mu_0) + \frac{1}{\sqrt{n_1}}\sum_{i=1}^{n_1} (\varepsilon_{1,i} - \varepsilon_{0,\hat\sigma_n(i)}),
\$
where $\epsilon_{1,i}=Y_{1,i}-\mu_1(T_i)$ and $\epsilon_{0,j} = Y_{0,j}-\mu_0(C_j)$ are the independent, mean-zero residuals.
The first term is a standard centered sample mean.
The second term is the imbalance in the (yet-unknown) outcome model $\mu_0$, but the uniform balance in Theorem~\ref{thm:observational} ensures it is root-$n$ negligible as long as the outcome model belongs to that $C^2$ smooth class.
Finally, since the design is outcome-blind and the matching is an injection without re-using control samples in the matching pairs, the third term is a summation over independent and mean-zero variables conditional on the matched design.
Finally, the (design-conditional) central limit theorem for the third term can be combined with the first two to obtain the desired result~\eqref{eq:ATT-main}.
\subsection{Counterfactual distributions and distributional effects}
\label{subsec:wald_cdf}
Apart from mean effects, distributional effects are widely used to assess the variation of the outcomes. We show that if the underlying distribution functions belong to the smooth function class, then the matched design also seamlessly supports the inference on distributional effects.
Let $\EE_z$ denote the expectation under the treatment-conditional distribution for $z\in \{0,1\}$. We define the conditional distribution function and the cumulative distribution function (in the treated)
\$
H_z(y\given x) = \PP(Y(z)\leq y\given X=x),\quad F_{z,T}(y) = \EE_1[H_z(y\given X)].
\$
The distributional effect (in the treated) is given by
\$
D(y) := F_{1,T}(y) - F_{0,T}(y).
\$
Here both $F_{1,T}$ and $F_{0,T}$ are defined by the covariate distribution in the treated group.
\paragraph{Point-wise Wald inference.} With the matched design $\hat\sigma_n$, the distributional effect can be estimated by
\@\label{eq:def_emp_cdf}
\hat{D}(y):= \hat{F}_1(y) - \hat{F}_0(y),\quad
\text{where}\quad
\hat{F}_1(y) = \frac{1}{n_1}\sum_{i=1}^{n_1} \ind\{Y_{1,i}\leq y\},\quad
\hat{F}_0(y) = \frac{1}{n_1}\sum_{i=1}^{n_1}\ind\{Y_{0,\hat\sigma_n(i)} \leq y\},
\@
whose variance can be estimated by
\$
\hat{V}_D(y) = \frac{1}{n_1-1} \sum_{i=1}^{n_1} \big\{ R_i(y) - \hat{D}(y)\big\}^2,\quad
R_i(y) = \ind\{Y_{1,i}\leq y\} - \ind\{Y_{0,\hat\sigma_n(i)} \leq y\}.
\$
We will show that $(\hat{F}_1(y),\hat{F}_0(y))$ converges to a Gaussian process, so pointwise Wald-type inference is convenient and valid. However, inference on confidence bands is arguable more convenient with bootstrap, which we present next.
\paragraph{Bootstrap inference.} The bootstrap inference will use the same resampling scheme introduced at the beginning of Section~\ref{sec:wald_inference}. For $b=1,\dots,B_n$, given the resampled pair indices $I_i^{*(b)}$, we define
\@\label{eq:def_cdf_boot}
\hat{D}^{*(b)} = \hat{F}_1^{*(b)} - \hat{F}_0^{*(b)},\quad
\hat F_1^{*(b)}(y)
=\frac1{n_1}\sum_{i=1}^{n_1}\ind\{Y_{1,I_i^{*(b)}}\le y\},~
\hat F_0^{*(b)}(y)
&=\frac1{n_1}\sum_{i=1}^{n_1}
\ind\{Y_{0,\hat\sigma_n(I_i^{*(b)})}\le y\}.
\@
To perform inference on the bands, we define
\$
T_F^{*(b)} = \sqrt{n_1} \max_{z\in \{0,1\}} \sup_{y\in \RR} \big|\hat{F}_z^{*(b)}(y) - \hat{F}_z(y)\big|,\quad
T_D^{*(b)} = \sqrt{n_1} \sup_{y\in \RR} \big| \hat{D}^{*(b)}(y) - \hat{D}(y)\big|.
\$
For a fixed $\alpha\in (0,1)$ and $H\in \{F,D\}$, we define $c_H^*$ as the empirical $(1-\alpha)$-quantile of $\{T_H^{*(b)}\}_{b=1}^{B_n}$.
Define the joint CDF confidence bands $\text{CI}_{z,n,1-\alpha}^{\text{joint}}$, $z\in \{0,1\}$ and the confidence band for the contrast $\text{CI}_{n,1-\alpha}^{\text{boot}}(y)$ as
\$
\text{CI}_{z,n,1-\alpha}^{\text{joint}} &= [L_z(y),U_z(y)] , \quad L_z(y) = \hat{F}_z(y)-c_F^*/\sqrt{n_1} ,\quad U_z(y) = \hat{F}_z(y)+c_F^*/\sqrt{n_1} , \\
\text{CI}_{n,1-\alpha}^{\text{contrast}}(y) &= [\hat{D}(y) - c_D^*/\sqrt{n_1}, \hat{D}(y) + c_D^*/\sqrt{n_1}].
\$
By default, $\text{CI}_{z,n,1-\alpha}^{\text{joint}}$ is truncated to $[0,1]$ and $\text{CI}_{n,1-\alpha}^{\text{contrast}}(y)$ is truncated to $[-1,1]$.
To state our results, we require the conditional CDF of the \textit{control} outcome to be smooth.
\begin{assumption}\label{assump:cdf}
The functions $H_0,H_1$ are jointly measurable, and there exists a fixed $L<\infty$ such that
$\{H_0(y\given\cdot):y\in\RR\}\subseteq\mathcal F_L$.
\end{assumption}
For a random variable $W\in \RR$ with CDF $G_W(\cdot)$, we say it has a \emph{regular $p$-quantile} if there exists some $c\in \RR$ such that $G_W(\cdot)$ is continuous at $c$, $G_W(c)=p$, and $G_W(c-\varepsilon)<p<G_W(c+\varepsilon)$ for every $\varepsilon>0$.
We present Wald inference for $D(y)$ in a pointwise fashion since the corresponding asymptotic variance can be explicitly estimated.
For confidence bands, we state the inferential results for bootstrap inference.
To describe the joint limiting distribution, we let $(\BB_1,\BB_0)$ be centered Gaussian process with covariance kernel
\@\label{dist:eq:kernel}
K_{zz}(s,t)&=F_{z,T}(s\wedge t)-F_{z,T}(s)F_{z,T}(t),\notag\\
K_{10}(s,t)&=\operatorname{Cov}_1\{H_1(s\given X),H_0(t\given X)\},
\qquad K_{01}(s,t)=K_{10}(t,s).
\@
Theorem~\ref{thm:wald_cdf} shows the convergence of the point estimates $(\hat{F}_1,\hat{F}_0)$, as well as the validity of Wald CI and bootstrap confidence bands. The proof is in Appendix~\ref{app:subsec_thm_wald_cdf}. The key idea is that the uniform balance in Theorem~\ref{thm:observational} applied to the conditional CDFs makes the matching bias negligible in the CDF inference.
\begin{theorem}
\label{thm:wald_cdf}
Under Assumptions~\ref{assump:obs-sample},~\ref{assump:match}, and~\ref{assump:cdf}, as $n\to \infty$,
\begin{equation}\label{dist:eq:CDFCLT}
\sqrt{n_1}(\hat F_1-F_{1,T},\hat F_0-F_{0,T})
\Rightarrow(\mathbb B_1,\mathbb B_0)
\quad\text{in }\ell^\infty(\mathbb R)^2.
\end{equation}
For each $y\in \RR$, the Wald interval $\textnormal{CI}_{n,1-\alpha}^{\textnormal{contrast,Wald}}(y) = \hat{D}(y) \pm z_{1-\alpha/2} \sqrt{\hat{V}_D(y)/n_1}$ has a limiting coverage of $1-\alpha$ for $D(y)$.
For bootstrap, suppose $B_n\to \infty$. Define
$W_F=\max_z\sup_y|\mathbb B_z(y)|$ and
$W_D=\sup_y|\mathbb B_1(y)-\mathbb B_0(y)|$.
If $W_F$ has a regular $(1-\alpha)$-quantile, then
\begin{equation}\label{dq:eq:joint-cdf-coverage}
\PP\big\{ F_{z,T}(y) \in \textnormal{CI}_{z,n,1-\alpha}^{\textnormal{joint}}
\text{ for every }y\in\RR,\ z\in\{0,1\}\big\}
\to 1-\alpha.
\end{equation}
If $W_D$ has a regular $(1-\alpha)$-quantile, then
\begin{equation}\label{dq:eq:direct-cdf-coverage}
\mathbb P\{D(y)\in\mathrm{CI}^{\textnormal{contrast}}_{n,1-\alpha}(y)
\text{ for every }y\in\mathbb R\}
\to1-\alpha.
\end{equation}
The coverage probabilities include both data and bootstrap sampling randomness.
\end{theorem}
\subsection{Quantile treatment effects}
\label{subsec:wald_quantile}
The same matched design also support inference on quantile treatment effects at regular quantiles.
For $u\in (0,1)$, the quantile treatment effect is defined as
\$
\Delta_Q(u) =Q_1(u)-Q_0(u),\quad \text{where}\quad Q_z(u) = \inf\{y\colon F_{z,T}(y)\geq u\}.
\$
An estimator is $\hat{\Delta}_Q(u) = \hat{Q}_1(u)-\hat{Q}_0(u)$, where $\hat{Q}_z(u)$ is the $u$-th quantile for the empirical CDFs $\hat{F}_1(\cdot)$ and $\hat{F}_0(\cdot)$ given in~\eqref{eq:def_emp_cdf}.
As the asymptotic variance in Wald-type inference requires involves the density function and requires additional estimation, we shall primarily focus on the convenient bootstrap inference.
Following the approach of~\cite{chernozhukov2020generic}, one can invert the CDF inference in Section~\ref{subsec:wald_cdf} for the quantile treatment effect. Namely, one can invert the upper and lower confidence bands for $F_{z,T}(y)$ obeying~\eqref{dq:eq:joint-cdf-coverage} to obtain a valid confidence bands for $\{\Delta_Q(u)\}_{u\in (0,1)}$. However, this approach can be conservative and may produce unbounded confidence bands near the tails.
Here, we present a direct bootstrap inference approach for regular quantiles to show the wide range of inferential targets the matched design can support. We use the same resampling scheme introduced earlier. For $b=1,\dots,B_n$, given the resampled pair indices $I_i^{*(b)}$, we define
\$
\hat{Q}_z^{*(b)}(u) = \inf\{y\colon \hat{F}_z^{*(b)}(y) \geq u\},\quad \hat\Delta_Q^{*(b)} = \hat{Q}_1^{*(b)}(u) - \hat{Q}_0^{*(b)}(u),
\quad S_Q^{*(b)}(u) = \sqrt{n_1} \{\hat\Delta_Q^{*(b)} - \hat\Delta_Q(u)\}.
\$
where the bootstrap empirical CDFs $\hat{F}_z^{*(b)}(\cdot)$ are given in~\eqref{eq:def_cdf_boot}.
The inference will be limited to regular quantiles. Consider a range $\cU = [u_-,u_+]$ with $0<u_-<u_+<1$.
\begin{assumption}\label{assump:quantile}
For each $z\in \{0,1\}$, there is a compact interval $J_z = [a_z,b_z]$ such that $F_{z,T}(a_z)<u_-$ and $F_{z,T}(b_z)>u_+$. On a neigborhood of $J_z$, $F_{z,T}$ is continuously differentiable and its derivative $f_{z,T}$ is bounded away from zero on $J_z$.
\end{assumption}
Let $c_Q^*$ be the empirical $(1-\alpha)$-quantile of $T_Q^{*(b)} = \sup_{u\in \cU}|S_Q^{*(b)}(u)|$ for $b=1,\dots,B_n$. We define the simultaneous confidence band
\$
\textnormal{CI}_{n,1-\alpha}^Q(u) = \big[\hat\Delta_Q(u) - c_Q^*/\sqrt{n_1}, \hat\Delta_Q(u) + c_Q^*/\sqrt{n_1}\big].
\$
For a fixed $u\in (0,1)$, let $c_{Q}^*(u)$ be the empirical $(1-\alpha)$-th quantile of $\{|S_Q^{*(b)}(u)|\}_{b=1}^{B_n}$.
We then define the pointwise bootstrap interval is
\$
\textnormal{CI}_{n,1-\alpha}^{Q,\textnormal{pt}}(u) = \big[\hat\Delta_Q(u) - c_{Q}^*(u)/\sqrt{n_1}, \hat\Delta_Q(u) + c_{Q}^*(u)/\sqrt{n_1}\big].
\$
The proof of Theorem~\ref{thm:quantile} is in Appendix~\ref{app:subsec_proof_quantile}.
\begin{theorem}\label{thm:quantile}
Under Assumptions~\ref{assump:obs-sample},~\ref{assump:match},~\ref{assump:cdf} and~\ref{assump:quantile}, as $n\to \infty$,
\begin{equation}\label{dist:eq:Qlimit}
\sqrt{n_1}(\hat\Delta_Q-\Delta_Q)
\Rightarrow\mathcal Z_Q
\quad\text{in }\ell^\infty(\mathcal U),\qquad
\mathcal Z_Q(u)=
-\frac{\mathbb B_1\{Q_1(u)\}}{f_{1,T}\{Q_1(u)\}}
+\frac{\mathbb B_0\{Q_0(u)\}}{f_{0,T}\{Q_0(u)\}},
\end{equation}
where the joint process $(\BB_1,\BB_0)$ is defined in Section~\ref{subsec:wald_cdf}.
Suppose $B_n\to\infty$. If
$W_Q=\sup_{u\in\mathcal U}|\mathcal Z_Q(u)|$
has a regular $(1-\alpha)$-quantile, then
\begin{equation}\label{dq:eq:quantile-coverage}
\PP\big\{\Delta_Q(u)\in\mathrm{CI}^{Q}_{n,1-\alpha}(u) \text{ for every }u\in\cU\big\} \to 1-\alpha.
\end{equation}
For every fixed $u\in\mathcal U$ with
$\operatorname{Var}\{\mathcal Z_Q(u)\}>0$,
\begin{equation}\label{dq:eq:quantile-pointwise-coverage}
\PP\big\{\Delta_Q(u)\in\mathrm{CI}^{Q,\mathrm{pt}}_{n,1-\alpha}(u)\big\} \to 1-\alpha.
\end{equation}
The above coverage probabilities include both the data and bootstrap sampling randomness.
\end{theorem}
\section{Asymptotic validity of randomization inference}
\label{sec:rand_inference}
Paired randomization tests are commonly used after matching in observational studies. They compare the observed statistic with a reference distribution obtained by independently exchanging the treated and control roles within each completed pair. Ideally, if the matching finds pairs with exactly the same covariate values, i.e., $T_i=C_{\hat\sigma_n(i)}$, then conditional on the covariates, the treatment is as-if randomized in the matched pair. However, with continuously-valued covariates, the matching is often inexact and leaves bias, in the sense that the treatment conditional on the matched pair is no long as-if uniform within the pair.
Recognizing that matching is inexact, \cite{pimentel2024covariate} proposed to modify the permutation probabilities using estimated
within-pair propensity-score discrepancies.
Here, our focus is instead when the randomization inference is asymptotically valid in the super-population sense, similar to~\cite{GR2023}, without
outcome adjustment or propensity-score modification.
\subsection{Paired randomization test of the weak and sharp null}
We consider two common hypotheses
\$
H_{0,\text{sharp}}\colon Y(1)=Y(0)~\text{almost surely},\qquad
H_{0,\text{weak}} \colon \theta_{\textnormal{ATT}} = \EE[Y(1)-Y(0)\given Z=1]=0.
\$
Throughout, we shall use the difference-in-mean test statistic.
Recall the paired difference $R_i$ and the estimator $\hat\theta_{\textnormal{ATT}}$ defined in~\eqref{eq:def_Ri_theta_att}.
Fix some $B_n\in \NN^+$. For $b=1,\dots,B_n$, given the observed data, we generate independent signs $\omega_i^{(b)}\sim\text{Unif}(\{-1,1\})$, independently of the data.
We then define the randomization p-values
\begin{equation}\label{eq:rimean:statistics}
p_{{\textnormal{RI}},n}
=\frac{1+\sum_{b=1}^{B_n}
\ind\{|T_{{\textnormal{RI}},n}^{\omega^{(b)}}| \ge |T_{{\textnormal{RI}},n}|\}}{B_n+1},
\qquad
T_{{\textnormal{RI}},n}= \hat\theta_{\mathrm{ATT}},
\qquad
T_{{\textnormal{RI}},n}^{\omega^{(b)}}
=\frac1{n_1}\sum_{i=1}^{n_1}\omega_i^{(b)}R_i,
\end{equation}
On the event $\cI_n^c$, we use the convention $p_{{\textnormal{RI}},n}=1$.
Theorem~\ref{thm:rand_inf} shows the validity of the difference-in-mean-based randomization inference, whose proof is in Appendix~\ref{app:subsec_proof_rand_inf}.
\begin{theorem}\label{thm:rand_inf}
Suppose Assumptions~\ref{assump:obs-sample},~\ref{assump:match},~\ref{assump:outcome} hold, and $B_n\to \infty$.
Under $H_{0,\mathrm{weak}}$, for every fixed
$\alpha\in(0,1)$,
\[
\mathbb P\{p_{{\textnormal{RI}},n}\le\alpha \}\to\alpha.
\]
The same conclusion holds under $H_{0,\mathrm{sharp}}$.
\end{theorem}
The probability in Theorem~\ref{thm:rand_inf} includes the randomness in both data sampling and randomization.
The key idea of the randomization inference validity is that, due to the negligible imbalance proved in Section~\ref{subsec:uniform_balance}, the bias in the difference-in-mean statistic caused by inexact matching is asymptotically negligible. Note that while randomization inference for sharp null, assuming exact matching, can in principle use any test statistic, our results are limited to the difference-in-mean statistic so that the uniform balance works.
\subsection{A dimension boundary at \texorpdfstring{$d=4$}{d = 4}}
The preceding conclusion can fail when matching bias is not negligible.
We continue to use the counterexample from Theorem~\ref{thm:dimension-obs}, keeping the covariate and treatment distribution and the function $g(x)$ and setting
\begin{equation}\label{eq:rimean:binary-model}
Y(1)=Y(0)=\ind\{U\le 1/4+g(X)/2\},
\qquad U\sim\operatorname{Unif}(0,1),\quad U\perp(X,Z).
\end{equation}
We consider the exact optimizer $\hat\sigma_n = \sigma_n^\star$ for $A_0=I_4$.
The next proposition shows that at $d=4$, the randomization p-value can fail to be asymptotically valid in this example.
\begin{prop}
\label{prop:rand_boundary}
Assume $B_n\to \infty$.
For the model \eqref{eq:rimean:binary-model}, for every fixed
$\alpha\in(0,1)$, there is a constant $\epsilon_\alpha>0$ such that
\[
\liminf_{n\to\infty}\mathbb P\{p_{{\textnormal{RI}},n}\le\alpha\}
\ge\alpha+\epsilon_\alpha.
\]
\end{prop}
The proof of Proposition~\ref{prop:rand_boundary} is in Appendix~\ref{app:subsec_rand_dimension}.
The key idea is as follows: the observed statistic is asymptotically equal to a centered outcome-noise
term plus the imbalance term
$\tfrac12\sqrt{n_1}B_{n,\sigma_n^\star}(g)$.
The randomization statistics approximate the centered noise distribution, but the imbalance term is non-negligible by Theorem~\ref{thm:dimension-obs} and translates to the violation of the type-I error control.
\section{Numerical demonstrations}
\label{sec:simu}
We corroborate the theory via simulations, where we vary the dimension of the covariates and evaluate the coverage of Wald and Bootstrap confidence intervals.
Below, we introduce the matching algorithms evaluated, whereas the data-generating processes and evaluation metrics are introduced separately in each subsection.
\vspace{-0.5em}
\paragraph{Matching algorithm.} Throughout the experiments, we evaluate both global matching and greedy matching. For each sample of size $n$, let $T_i$ and $C_j$ denote the covariates of
the $n_1$ treated and $n_0$ control observations, respectively. We estimate $\hat A$ as the inverse covariance of all $n$
covariate vectors, using denominator $n-1$, and form the
cost matrix $c_{ij}=\tfrac12(T_i-C_j)^\top\hat A(T_i-C_j)$. Global matching casts the optimization as a linear assignment problem, and uses SciPy's floating-point assignment solver to minimize
$Q_n(\sigma;\hat A)=n_1^{-1}\sum_i c_{i,\sigma(i)}$ over injections, which is guaranteed to find the globally optimal solution.
The Greedy matching approach visits treated observations in an independently drawn random
order and assigns each the remaining unmatched control with the lowest matching cost; this approach does not necessarily find the optimal solution and is shown as a robustness check.
\subsection{Population ATT inference}
\label{subsec:simu_att}
We first evaluate population ATT inference established in Section~\ref{subsec:wald_att}.
We generate $n$ independent observations with
$X_i\sim\operatorname{Unif}([0,1]^d)$ and
$Z_i\mid X_i\sim\operatorname{Bernoulli}\{e(X_i)\}$, where
$e(x)=0.1+0.3x_1$. Potential outcomes are generated by
$Y_i(z)=\mu_z(X_i)+(0.5+0.5X_{i1})\varepsilon_{iz}$ for $z\in\{0,1\}$,
where the errors $\{\varepsilon_{iz}\}$ are independent standard normal,
independent of the covariates and treatment indicators. We set nonlinear outcomes
$\mu_0(x) =2x_1+\frac{2}{\sqrt d}\sum_{j=1}^d(x_j^2-1/3)+g_d(x)$,
$\mu_1(x)=\mu_0(x)+a(x_1-3/5)$, where $a\in \{0,2\}$ controls the heterogeneity of the conditional average treatment effect (CATE) function, with $g_d(x)=0$ when $d=1$ and $g_d(x) = 2(x_1-1/2)(x_2-1/2)$ when $d\geq 2$.
Since the matching design is outcome-blind, the estimation and inference require no fitting or adjustment based on outcome models. It also does not require estimating the propensity score.
We vary the total sample size $n\in\{500,2000,5000\}$ and covariate dimension $d\in\{1,2,3,4,6\}$.
The population ATT
is $\theta=0$ under both heterogeneity settings.
Given the matching $\hat\sigma_n$ by the two matching algorithms introduced at the beginning of this section, we compute the paired differences
$R_i=Y_{1i}-Y_{0,\hat\sigma_n(i)}$, ATT estimator, Wald-type and bootstrap (where we set $B_n=500$) confidence intervals at the $1-\alpha=95\%$ nominal level as introduced in Section~\ref{subsec:wald_att}.
For each matching algorithm and each inference method (Wald versus bootstrap), we evaluate the empirical coverage and CI length over $R=2000$ independent runs (including data generation, matching, and inference).
We summarize the results for the heterogeneous CATE case with $a=2$ in Table~\ref{tab:att-a2}, while the results for the constant CATE model with $a=0$ are deferred to Table~\ref{tab:att-a0} in Appendix~\ref{app:sec_simu} with qualitatively similar messages.
Consistent with our theory, with both global and greedy matching, both Wald and bootstrap CI achieve close to $95\%$ empirical coverage for $d\leq 3$, with the gap diminishing as the sample size grows to $n=5000$. This shows the particular benefit of the matching design: even with nonlinear outcome models, as long as $d\leq 3$, matching without outcome information still delivers valid inference.
In contrast, when $d\in \{4,6\}$, the coverage gap does not vanish as we increase the sample size: it appears to converge to some probability below $0.95$ for $d=4$ and keeps deteriorating for $d=6$. This suggests that $d=4$ can indeed be the dimension boundary beyond which the bias of inexact matching becomes non-negligible.
\begin{table}[p]
\centering
\caption{Population ATT inference results with $a=2$ (heterogeneous CATE). Coverage and its Monte Carlo SE (parentheses) are in percentage. Results are averaged over $R=2000$ independent runs.}
\label{tab:att-a2}
\small\setlength{\tabcolsep}{6pt}\begin{tabular}{rrlrrrr}
\toprule
&&&\multicolumn{2}{c}{Coverage (\%)}&\multicolumn{2}{c}{Mean CI width}\\
\cmidrule(lr){4-5}\cmidrule(lr){6-7}
$n$ & $d$ & Matching & Wald & Bootstrap-t & Wald & Bootstrap-t \\
\midrule
500 & 1 & Global & 94.65 (0.50) & 95.00 (0.49) & 0.4452 (0.0007) & 0.4493 (0.0009)\\
& & Greedy & 94.80 (0.50) & 94.60 (0.51) & 0.4452 (0.0007) & 0.4490 (0.0009)\\
\addlinespace[2pt]
& 2 & Global & 94.70 (0.50) & 95.25 (0.48) & 0.4483 (0.0007) & 0.4532 (0.0009)\\
& & Greedy & 94.30 (0.52) & 94.10 (0.53) & 0.4489 (0.0008) & 0.4539 (0.0009)\\
\addlinespace[2pt]
& 3 & Global & 93.90 (0.54) & 93.95 (0.53) & 0.4566 (0.0008) & 0.4611 (0.0009)\\
& & Greedy & 93.75 (0.54) & 93.55 (0.55) & 0.4575 (0.0008) & 0.4622 (0.0009)\\
\addlinespace[2pt]
& 4 & Global & 91.20 (0.63) & 91.10 (0.64) & 0.4661 (0.0008) & 0.4711 (0.0009)\\
& & Greedy & 91.55 (0.62) & 90.80 (0.65) & 0.4679 (0.0008) & 0.4731 (0.0009)\\
\addlinespace[2pt]
& 6 & Global & 85.25 (0.79) & 84.95 (0.80) & 0.4890 (0.0008) & 0.4934 (0.0010)\\
& & Greedy & 85.70 (0.78) & 85.50 (0.79) & 0.4914 (0.0008) & 0.4962 (0.0010)\\
\midrule
2000 & 1 & Global & 93.65 (0.55) & 93.35 (0.56) & 0.2224 (0.0002) & 0.2228 (0.0003)\\
& & Greedy & 93.50 (0.55) & 93.45 (0.55) & 0.2224 (0.0002) & 0.2227 (0.0003)\\
\addlinespace[2pt]
& 2 & Global & 94.65 (0.50) & 94.65 (0.50) & 0.2229 (0.0002) & 0.2234 (0.0003)\\
& & Greedy & 94.20 (0.52) & 94.00 (0.53) & 0.2231 (0.0002) & 0.2235 (0.0003)\\
\addlinespace[2pt]
& 3 & Global & 94.50 (0.51) & 94.25 (0.52) & 0.2247 (0.0002) & 0.2246 (0.0003)\\
& & Greedy & 94.30 (0.52) & 94.25 (0.52) & 0.2250 (0.0002) & 0.2251 (0.0003)\\
\addlinespace[2pt]
& 4 & Global & 91.10 (0.64) & 90.80 (0.65) & 0.2280 (0.0002) & 0.2280 (0.0003)\\
& & Greedy & 91.30 (0.63) & 90.75 (0.65) & 0.2285 (0.0002) & 0.2286 (0.0003)\\
\addlinespace[2pt]
& 6 & Global & 78.40 (0.92) & 77.75 (0.93) & 0.2364 (0.0002) & 0.2367 (0.0003)\\
& & Greedy & 79.20 (0.91) & 79.00 (0.91) & 0.2371 (0.0002) & 0.2374 (0.0003)\\
\midrule
5000 & 1 & Global & 95.00 (0.49) & 94.85 (0.49) & 0.1407 (0.0001) & 0.1407 (0.0002)\\
& & Greedy & 94.95 (0.49) & 94.80 (0.50) & 0.1407 (0.0001) & 0.1407 (0.0002)\\
\addlinespace[2pt]
& 2 & Global & 94.70 (0.50) & 94.55 (0.51) & 0.1408 (0.0001) & 0.1407 (0.0002)\\
& & Greedy & 94.75 (0.50) & 94.60 (0.51) & 0.1408 (0.0001) & 0.1407 (0.0002)\\
\addlinespace[2pt]
& 3 & Global & 94.05 (0.53) & 93.95 (0.53) & 0.1415 (0.0001) & 0.1415 (0.0002)\\
& & Greedy & 93.65 (0.55) & 93.35 (0.56) & 0.1416 (0.0001) & 0.1417 (0.0002)\\
\addlinespace[2pt]
& 4 & Global & 91.70 (0.62) & 90.90 (0.64) & 0.1430 (0.0001) & 0.1430 (0.0002)\\
& & Greedy & 91.70 (0.62) & 91.65 (0.62) & 0.1431 (0.0001) & 0.1432 (0.0002)\\
\addlinespace[2pt]
& 6 & Global & 71.50 (1.01) & 70.80 (1.02) & 0.1470 (0.0001) & 0.1469 (0.0002)\\
& & Greedy & 72.05 (1.00) & 71.60 (1.01) & 0.1474 (0.0001) & 0.1472 (0.0002)\\
\bottomrule\end{tabular}
\par\vspace{7pt}
\end{table}
\subsection{Distributional treatment effects}
We then evaluate the distributional inference based on a matched design in Section~\ref{subsec:wald_cdf}.
We generate $n$ independent observations with
$X_i\sim\operatorname{Unif}([0,1]^d)$ and
$Z_i\mid X_i\sim\operatorname{Bernoulli}\{e(X_i)\}$, where the propensity score is $e(x)=0.10+\frac{0.30}{1+\exp\{-6h_d(x)\}}$ for
$h_d(x)=\frac{1}{\sqrt d}\sum_{j=1}^d(x_j^2-1/3)$.
The potential outcomes are generated as $Y_i(z)=4h_d(X_i)+\varepsilon_{iz}$ for
$z\in\{0,1\}$, where the errors are independent standard normal and
independent of everything else. Thus
$H_1(y\given x)=H_0(y\given x)=\Phi\{y-4h_d(x)\}$ and the population
contrast is $D(y)=F_{1,T}(y)-F_{0,T}(y)=0$ for every $y$.
We vary the total sample size $n\in\{500,2000,5000\}$ and covariate dimension $d\in\{2,3,4,6\}$.
Given the matching $\hat\sigma_n$ by the two algorithms described, we form the bootstrap (with $B_n=500$) joint confidence bands and confidence band for the contrast $D(y)$ at the $1-\alpha=95\%$ nominal level as described in Section~\ref{subsec:wald_cdf}.
For each matching algorithm and each type of confidence bands, we evaluate the empirical coverage and mean band width over $R=2000$ independent runs.
Table~\ref{tab:cdf-bands} summarizes the results. As expected, the coverage is close to nominal for $d\leq 3$. For $d=4$, as we increase the sample size $n$, the coverage tends to a constant level below $0.95$. The coverage for $d=6$ keeps deteriorating as $n$ increases. This confirms our theory, and the uniform balance violation at $d\geq 4$ appears to play a role in this experiment.
\begin{table}[!htbp]
\centering
\caption{Coverage and band width in CDF inference. Joint coverage requires both $F_{0,T}$ and $F_{1,T}$ to be covered
at every threshold; contrast coverage requires $D=F_{1,T}-F_{0,T}$ to
be covered at every threshold. Coverage and its Monte Carlo SE (parentheses) are in percentage. The confidence bands are truncated to $[0,1]$ before computing the mean width. Results are averaged over $R=2000$ independent runs.}
\label{tab:cdf-bands}
\small\setlength{\tabcolsep}{6pt}\begin{tabular}{rrlrrrr}
\toprule
& & & \multicolumn{2}{c}{Coverage (\%)} & \multicolumn{2}{c}{Mean band width} \\
\cmidrule(lr){4-5}\cmidrule(lr){6-7}
$n$ & $d$ & Matching & Joint CDF & Contrast & Joint CDF & Contrast \\
\midrule
500 & 1 & Global & 94.00 (0.53) & 96.60 (0.41) & 0.2566 (0.00028) & 0.2740 (0.00036) \\
& & Greedy & 93.65 (0.55) & 96.70 (0.40) & 0.2568 (0.00028) & 0.2742 (0.00036) \\
\addlinespace[2pt]
& 2 & Global & 94.00 (0.53) & 95.35 (0.47) & 0.2546 (0.00028) & 0.2788 (0.00037) \\
& & Greedy & 93.80 (0.54) & 95.40 (0.47) & 0.2546 (0.00028) & 0.2788 (0.00037) \\
\addlinespace[2pt]
& 3 & Global & 92.35 (0.59) & 94.75 (0.50) & 0.2543 (0.00028) & 0.2816 (0.00038) \\
& & Greedy & 92.20 (0.60) & 94.50 (0.51) & 0.2544 (0.00028) & 0.2820 (0.00037) \\
\addlinespace[2pt]
& 4 & Global & 90.10 (0.67) & 92.35 (0.59) & 0.2550 (0.00028) & 0.2857 (0.00039) \\
& & Greedy & 90.25 (0.66) & 92.30 (0.60) & 0.2547 (0.00028) & 0.2855 (0.00038) \\
\addlinespace[2pt]
& 6 & Global & 77.95 (0.93) & 83.95 (0.82) & 0.2542 (0.00027) & 0.2913 (0.00036) \\
& & Greedy & 78.85 (0.91) & 83.95 (0.82) & 0.2543 (0.00026) & 0.2913 (0.00037) \\
\addlinespace[2pt]
2000 & 1 & Global & 94.90 (0.49) & 95.55 (0.46) & 0.1310 (0.00009) & 0.1382 (0.00011) \\
& & Greedy & 95.20 (0.48) & 96.00 (0.44) & 0.1311 (0.00009) & 0.1381 (0.00011) \\
\addlinespace[2pt]
& 2 & Global & 94.00 (0.53) & 96.00 (0.44) & 0.1303 (0.00009) & 0.1411 (0.00011) \\
& & Greedy & 93.80 (0.54) & 95.80 (0.45) & 0.1303 (0.00009) & 0.1411 (0.00011) \\
\addlinespace[2pt]
& 3 & Global & 93.65 (0.55) & 94.25 (0.52) & 0.1299 (0.00009) & 0.1420 (0.00012) \\
& & Greedy & 93.60 (0.55) & 94.85 (0.49) & 0.1297 (0.00009) & 0.1420 (0.00012) \\
\addlinespace[2pt]
& 4 & Global & 91.30 (0.63) & 91.95 (0.61) & 0.1299 (0.00009) & 0.1432 (0.00012) \\
& & Greedy & 91.30 (0.63) & 91.80 (0.61) & 0.1299 (0.00009) & 0.1433 (0.00012) \\
\addlinespace[2pt]
& 6 & Global & 68.20 (1.04) & 67.50 (1.05) & 0.1298 (0.00009) & 0.1457 (0.00012) \\
& & Greedy & 69.10 (1.03) & 68.55 (1.04) & 0.1298 (0.00009) & 0.1456 (0.00012) \\
\addlinespace[2pt]
5000 & 1 & Global & 94.50 (0.51) & 95.35 (0.47) & 0.0835 (0.00005) & 0.0876 (0.00005) \\
& & Greedy & 94.55 (0.51) & 95.50 (0.46) & 0.0834 (0.00005) & 0.0877 (0.00006) \\
\addlinespace[2pt]
& 2 & Global & 94.30 (0.52) & 94.40 (0.51) & 0.0829 (0.00005) & 0.0896 (0.00006) \\
& & Greedy & 94.55 (0.51) & 94.55 (0.51) & 0.0830 (0.00005) & 0.0896 (0.00006) \\
\addlinespace[2pt]
& 3 & Global & 94.10 (0.53) & 93.15 (0.56) & 0.0828 (0.00005) & 0.0901 (0.00006) \\
& & Greedy & 93.85 (0.54) & 93.40 (0.56) & 0.0828 (0.00005) & 0.0901 (0.00006) \\
\addlinespace[2pt]
& 4 & Global & 90.65 (0.65) & 89.40 (0.69) & 0.0828 (0.00005) & 0.0906 (0.00006) \\
& & Greedy & 90.70 (0.65) & 88.80 (0.71) & 0.0828 (0.00005) & 0.0906 (0.00006) \\
\addlinespace[2pt]
& 6 & Global & 57.15 (1.11) & 56.30 (1.11) & 0.0827 (0.00005) & 0.0920 (0.00006) \\
& & Greedy & 57.65 (1.10) & 57.40 (1.11) & 0.0827 (0.00005) & 0.0920 (0.00006) \\
\bottomrule
\end{tabular}
\end{table}
\subsection{Quantile treatment effects}
We further demonstrate the inference for quantile treatment effects in Section~\ref{subsec:wald_quantile}.
We use the same covariate distribution $X\sim\operatorname{Unif}([0,1]^d)$ and treatment assignment
mechanism $e(x)$ as in the CDF experiment.
The potential outcomes are generated as
$Y_i(0)=4h(X_i)+\varepsilon_{i0},$ and
$Y_i(1)=0.5+4h(X_i)+\varepsilon_{i1}$, where the errors are independent standard normal and independent of everything else. The target is the quantile treatment effect
$\Delta_Q(u)=F_{1,T}^{-1}(u)-F_{0,T}^{-1}(u)$ for
$u\in U=[0.1,0.9]$.
We consider the total sample size $n\in\{500,2000,5000\}$ and covariate dimension $d\in\{1,2,3,4,6\}$.
Given the matching $\hat\sigma_n$ from the described global and greedy procedures, we construct the bootstrap simultaneous confidence bands for $\textnormal{CI}_{n,1-\alpha}^{Q}(u)$ as stated in Section~\ref{subsec:wald_quantile} over $B_n=500$ boostrap draws.
We then report coverage over the entire curve on $U$ and the mean width of the confidence band, averaged over $R=2000$ independent runs as before.
The simultaneous coverage is reported in Table~\ref{tab:quantile-risk-aligned} column named \texttt{Raw cov.}, where we observe valid coverage for $d\leq 3$ and deteriorating coverage for $d=6$ when $n$ increases.
However, the coverage appears conservative: it remains above $97\%$ for $d\leq 4$ when $n$ is as large as $5000$.
To diagnose this conservativeness, we construct an oracle version (column name \texttt{Oracle cov.}) by removing the matching-induced bias in the quantile treatment effects, and observe the over-coverage for all dimensions; this confirms that the inference without matching bias is itself conservative at this sample size.
In addition, we measure the matching-induced bias by averaging $D_Q=\sup_{u\in U}|\delta_n(u)|$ and $\sqrt{n_1}D_Q$, where $Q_{T,n}$ and $Q_{M,n}$
are the quantile functions of the unit-variance Gaussian mixtures
with means $4h(X_i)$ over treated and matched-control covariates,
respectively. These quantities are reported in the column \texttt{Mean design bias}.
The results confirm our theory: the scaled bias $\sqrt{n_1} D_Q$ appears to diminish for $d\leq 3$ (although a bit slow for $d=3$), converge for $d=4$, while keeps increasing for $d=6$.
Together, these results demonstrate the dimension boundary for the imbalance, yet in this example, the effect on simultaneous coverage is diluted especially for $d=4$ due to the intrinsic conservativeness of the inference procedure.
Finally, we report the coverage and width of point-wise confidence intervals for the $u$-th quantile treatment effects introduced in Section~\ref{subsec:wald_quantile}, for $u\in \{0.25, 0.5, 0.75\}$ in Table~\ref{tab:pointwise-qte} in Appendix~\ref{app:sec_simu}. The point-wise confidence intervals converge to the nominal level faster than simultaneous ones for $d\leq 3$, and illustrates the systematic undercoverage for $d=4$ and deteriorating coverage for $d=6$.
\begin{table}[!htbp]
\centering
\small\setlength{\tabcolsep}{4pt}\caption{Simultaneous coverage for quantile treatment effects. Coverage is for the entire $U=[0.1,0.9]$, in percent, with Monte Carlo SEs in parentheses.
Oracle coverage is a diagnostic which removes the matching-induced bias from the construction of confidence bands.}
\label{tab:quantile-risk-aligned}
\begin{tabular}{rrlccccc}
\toprule
&&&\multicolumn{3}{c}{Direct simultaneous band}&\multicolumn{2}{c}{Mean design bias}\\
\cmidrule(lr){4-6}\cmidrule(lr){7-8}
$n$&$d$&Matching&Raw cov. (\%) & Oracle cov. (\%) &Width&$D_Q$&$\sqrt{n_1}D_Q$\\
\midrule
500 & 1 & Global & 99.90 (0.07) & 99.85 (0.09) & 1.8362 (0.0040) & 0.00862 (0.00033) & 0.0956 (0.0037) \\
& & Greedy & 99.75 (0.11) & 99.70 (0.12) & 1.8412 (0.0041) & 0.00964 (0.00033) & 0.1069 (0.0037) \\
\addlinespace[2pt]
& 2 & Global & 99.30 (0.19) & 99.60 (0.14) & 1.8309 (0.0045) & 0.04471 (0.00085) & 0.4979 (0.0095) \\
& & Greedy & 99.30 (0.19) & 99.35 (0.18) & 1.8333 (0.0045) & 0.04540 (0.00083) & 0.5054 (0.0093) \\
\addlinespace[2pt]
& 3 & Global & 99.55 (0.15) & 99.70 (0.12) & 1.8313 (0.0046) & 0.09649 (0.00120) & 1.0747 (0.0136) \\
& & Greedy & 99.45 (0.17) & 99.75 (0.11) & 1.8349 (0.0046) & 0.09558 (0.00120) & 1.0648 (0.0135) \\
\addlinespace[2pt]
& 4 & Global & 99.15 (0.21) & 99.40 (0.17) & 1.8423 (0.0045) & 0.15228 (0.00150) & 1.6937 (0.0169) \\
& & Greedy & 99.15 (0.21) & 99.65 (0.13) & 1.8435 (0.0046) & 0.14997 (0.00149) & 1.6681 (0.0168) \\
\addlinespace[2pt]
& 6 & Global & 97.65 (0.34) & 99.95 (0.05) & 1.8600 (0.0043) & 0.24923 (0.00174) & 2.7746 (0.0196) \\
& & Greedy & 97.35 (0.36) & 99.80 (0.10) & 1.8610 (0.0043) & 0.24762 (0.00174) & 2.7565 (0.0196) \\
\midrule
2000 & 1 & Global & 99.10 (0.21) & 99.10 (0.21) & 0.8899 (0.0012) & 0.00096 (0.00004) & 0.0212 (0.0009) \\
& & Greedy & 99.35 (0.18) & 99.35 (0.18) & 0.8900 (0.0012) & 0.00119 (0.00004) & 0.0262 (0.0009) \\
\addlinespace[2pt]
& 2 & Global & 99.15 (0.21) & 99.20 (0.20) & 0.8877 (0.0014) & 0.01150 (0.00019) & 0.2552 (0.0042) \\
& & Greedy & 99.00 (0.22) & 98.90 (0.23) & 0.8871 (0.0014) & 0.01149 (0.00020) & 0.2550 (0.0044) \\
\addlinespace[2pt]
& 3 & Global & 98.85 (0.24) & 98.90 (0.23) & 0.8799 (0.0014) & 0.03833 (0.00039) & 0.8526 (0.0086) \\
& & Greedy & 98.50 (0.27) & 99.00 (0.22) & 0.8804 (0.0014) & 0.03718 (0.00039) & 0.8271 (0.0087) \\
\addlinespace[2pt]
& 4 & Global & 97.80 (0.33) & 99.30 (0.19) & 0.8833 (0.0014) & 0.07535 (0.00055) & 1.6772 (0.0122) \\
& & Greedy & 98.25 (0.29) & 99.20 (0.20) & 0.8841 (0.0014) & 0.07304 (0.00055) & 1.6259 (0.0122) \\
\addlinespace[2pt]
& 6 & Global & 92.45 (0.59) & 99.40 (0.17) & 0.8862 (0.0014) & 0.15845 (0.00078) & 3.5290 (0.0175) \\
& & Greedy & 91.80 (0.61) & 99.45 (0.17) & 0.8879 (0.0014) & 0.15515 (0.00077) & 3.4557 (0.0174) \\
\midrule
5000 & 1 & Global & 99.00 (0.22) & 99.00 (0.22) & 0.5485 (0.0006) & 0.00023 (0.00001) & 0.0079 (0.0003) \\
& & Greedy & 98.80 (0.24) & 98.85 (0.24) & 0.5486 (0.0006) & 0.00029 (0.00001) & 0.0101 (0.0003) \\
\addlinespace[2pt]
& 2 & Global & 98.60 (0.26) & 98.55 (0.27) & 0.5480 (0.0006) & 0.00474 (0.00007) & 0.1662 (0.0026) \\
& & Greedy & 98.85 (0.24) & 98.80 (0.23) & 0.5481 (0.0006) & 0.00465 (0.00007) & 0.1632 (0.0025) \\
\addlinespace[2pt]
& 3 & Global & 98.30 (0.29) & 98.65 (0.26) & 0.5450 (0.0006) & 0.02051 (0.00019) & 0.7208 (0.0065) \\
& & Greedy & 98.55 (0.27) & 98.85 (0.24) & 0.5453 (0.0006) & 0.01976 (0.00019) & 0.6942 (0.0066) \\
\addlinespace[2pt]
& 4 & Global & 97.25 (0.37) & 98.75 (0.25) & 0.5443 (0.0006) & 0.04776 (0.00029) & 1.6792 (0.0101) \\
& & Greedy & 97.80 (0.33) & 98.95 (0.23) & 0.5446 (0.0006) & 0.04590 (0.00029) & 1.6140 (0.0101) \\
\addlinespace[2pt]
& 6 & Global & 83.95 (0.82) & 98.85 (0.24) & 0.5471 (0.0006) & 0.11701 (0.00044) & 4.1194 (0.0155) \\
& & Greedy & 85.10 (0.80) & 98.90 (0.23) & 0.5477 (0.0006) & 0.11387 (0.00043) & 4.0086 (0.0153) \\
\bottomrule
\end{tabular}
\end{table}
\subsection{Randomization inference}
Finally, we demonstrate the super-population validity of randomization inference for $d\leq 3$ and the dimension boundary for $d\geq 4$.
We use the same covariate distribution, propensity score, outcome regression, and matching procedures as in the ATT experiment in Section~\ref{subsec:simu_att}.
We consider two null settings. Under the weak null, we keep the heterogeneous CATE with $a=2$, so that $\theta_{\textnormal{ATT}}=0$.
Under the sharp null, we instead set
$Y_i(1)=Y_i(0)=\mu_0(X_i)+s(X_i)\eta_i$,
where the $\eta_i$'s are independent standard normal variables and independent of everything else, leading to the same observed distribution as the $a=0$ case. We vary the total sample size $n\in \{500, 2000, 5000\}$ and the covariate dimension $d\in \{1,2,3,4,6\}$.
Given the matching $\hat\sigma_n$ by the two algorithms, we form the randomization p-value in Section~\ref{sec:rand_inference} with $B_n=4999$. We then evaluate the type-I error at nominal levels $\alpha\in \{0.025, 0.05, 0.1\}$.
The results for the weak null are in Table~\ref{tab:ri-weak-heterogeneous-alphas}, while we defer those for the sharp null to Table~\ref{tab:ri-sharp-shared-alphas} in Appendix~\ref{app:sec_simu}.
Across both null settings and various nominal levels, both matching methods yield close-to-nominal rejection rates for $d\leq 3$. The rejection rate for $d=4$ appears to converge to some value above the nominal level, and that for $d=6$ keeps increasing as the total sample size $n$ increases. This is consistent with our theory on asymptotic, super-population validity of randomization inference based on the difference-in-mean statistic.
\begin{table}
\centering
\small
\setlength{\tabcolsep}{4.5pt}
\caption{Rejection rates in randomization inference under the weak null. The rejection rates and its Monte Carlo SE (parentheses) are in percentage, averaged over $R=2000$ independent runs.}
\label{tab:ri-weak-heterogeneous-alphas}
\begin{tabular}{@{}rr*{6}{r}@{}}
\toprule
& & \multicolumn{3}{c}{Global} & \multicolumn{3}{c}{Greedy} \\
\cmidrule(lr){3-5}\cmidrule(lr){6-8}
$n$ & $d$ & $\alpha=0.025$ & $\alpha=0.05$ & $\alpha=0.10$ & $\alpha=0.025$ & $\alpha=0.05$ & $\alpha=0.10$ \\
\midrule
500 & 1 & 2.65 (0.36) & 5.10 (0.49) & 9.00 (0.64) & 2.70 (0.36) & 4.90 (0.48) & 9.30 (0.65) \\
500 & 2 & 2.60 (0.36) & 5.10 (0.49) & 10.70 (0.69) & 2.65 (0.36) & 5.70 (0.52) & 10.20 (0.68) \\
500 & 3 & 3.15 (0.39) & 5.90 (0.53) & 10.95 (0.70) & 3.10 (0.39) & 5.90 (0.53) & 10.75 (0.69) \\
500 & 4 & 5.00 (0.49) & 8.65 (0.63) & 15.05 (0.80) & 4.60 (0.47) & 8.25 (0.62) & 14.70 (0.79) \\
500 & 6 & 9.15 (0.64) & 14.45 (0.79) & 23.05 (0.94) & 9.00 (0.64) & 13.95 (0.77) & 22.35 (0.93) \\
\addlinespace
2000 & 1 & 3.60 (0.42) & 6.30 (0.54) & 11.40 (0.71) & 3.75 (0.42) & 6.80 (0.56) & 11.80 (0.72) \\
2000 & 2 & 3.05 (0.38) & 5.25 (0.50) & 10.80 (0.69) & 3.30 (0.40) & 5.65 (0.52) & 11.10 (0.70) \\
2000 & 3 & 2.50 (0.35) & 5.50 (0.51) & 10.55 (0.69) & 2.80 (0.37) & 5.60 (0.51) & 10.45 (0.68) \\
2000 & 4 & 4.60 (0.47) & 8.80 (0.63) & 15.25 (0.80) & 3.80 (0.43) & 8.35 (0.62) & 14.85 (0.80) \\
2000 & 6 & 14.15 (0.78) & 21.50 (0.92) & 32.65 (1.05) & 13.00 (0.75) & 20.60 (0.90) & 31.40 (1.04) \\
\addlinespace
5000 & 1 & 2.75 (0.37) & 4.75 (0.48) & 9.95 (0.67) & 2.50 (0.35) & 5.15 (0.49) & 9.85 (0.67) \\
5000 & 2 & 2.70 (0.36) & 5.10 (0.49) & 10.50 (0.69) & 2.60 (0.36) & 5.45 (0.51) & 10.85 (0.70) \\
5000 & 3 & 3.15 (0.39) & 5.85 (0.52) & 11.55 (0.71) & 3.45 (0.41) & 6.40 (0.55) & 11.75 (0.72) \\
5000 & 4 & 4.60 (0.47) & 8.35 (0.62) & 15.15 (0.80) & 5.00 (0.49) & 8.05 (0.61) & 15.25 (0.80) \\
5000 & 6 & 20.25 (0.90) & 28.65 (1.01) & 40.75 (1.10) & 19.60 (0.89) & 27.75 (1.00) & 39.35 (1.09) \\
\bottomrule
\end{tabular}
\end{table}
\section{Related literature}\label{sec:literature}
This work adds to the large-sample theory of matching, especially one-to-one matching without replacement. We discuss related literature in several strands.
\vspace{-0.5em}
\paragraph{Matching with and without replacement.}
Much of the large-sample theory for matching concerns nearest-neighbor (NN) matching with replacement. \citet{AI2006} establish asymptotic and bias results for estimators with a fixed number of matches, and \citet{AI2011} develop regression-based bias correction. In these procedures, control samples can be reused for different treated units. More recent work has revisited NN matching from
the perspective of increasing numbers of matches and semiparametric efficiency~\citep{lin2023estimation,he2024propensity,cattaneo2025rosenbaum,lin2026consistency}. In particular, \citet{lin2023estimation} shows that the matching nearest neighbors essentially estimates a density ratio, and with regression adjustment, NN matching can be doubly robust and semiparametrically efficient.
In contrast, we study one-to-one matching without replacement, a classical design
\citep{rosenbaum1989optimal} with applied precedents for the squared Mahalanobis
objective~\citep{bind2019bridging,zhang2023matching}.
Assigning a control to one treated observation changes the controls available to the others. Consequently, the local nearest-neighbor calculations do not directly apply in our setting, which necessitates distinct analyses.
\vspace{-0.5em}
\paragraph{Error cancellation, bias, and balance in matching without replacement}
This work expands upon the classical bias-reduction theory for matching. Under ellipsoidal symmetry and related mixture assumptions, \citet{RT1992} and \citet{RS2006} establish equal-percent bias-reduction properties for affinely invariant matching. These results describe how matching improves covariate comparability but do not study the rate for empirical imbalance. Similarly, the balancing-score identity of \citet{RR1983} provides a population-level perspective without controlling the finite-sample discrepancy.
For matching without replacement, even consistency depends on the availability of controls in the regions occupied by treated observations. \citet{Savje2022} shows that optimal matching on the propensity score can be inconsistent when some regions have too few controls. We rule out this complication by assuming the propensity score is strictly bounded below $1/2$, providing a clean basis for analyzing the matching bias in large samples.
\citet{GR2023} study the finite-sample matching bias using ordinary Mahalanobis distance, similarly assuming the propensity score is strictly bounded below $1/2$. They establish asymptotic validity of an unadjusted paired randomization test with one continuous covariate, conjecture validity with up to three, and give a failure example for $d=4$. Our results establish a quadratic counterpart, yielding dimensionality results echoing their conjecture.
Second-order bias and the relevance of dimension four also have nearest-neighbor precedents. The recent work of \citet{VTY2025} derive higher-order bias bounds and expansions under geometric and smoothness conditions in the NN matching setting, yet matching without replacement requires distinct analysis techniques. Our contribution is to show that
the global use-once constraint does not create a first-order bias, and to control
the resulting signed imbalance uniformly.
\vspace{-0.5em}
\paragraph{Inference from matched samples.}
Inference after matching without replacement already has a substantial foundation. \citet{AI2012} derive a martingale representation, establish an asymptotic normal distribution for a without-replacement matching estimator, and give a consistent paired-difference variance estimator. \citet{AS2022} develop post-matching inference that allows regression misspecification and justify resampling complete matched sets (similar to our bootstrap procedures). These works isolate the approximation error from inexact matching, often assuming sufficiently small absolute distance, defined as
$A_n=\frac1m\sum_{i=1}^m\|T_i-C_{\sigma(i)}\|$.
A sufficient condition for the results in \citet[Theorem~1]{AI2012} and \citet[Assumption~3]{AS2022} is $\sqrt m\,A_n = o_P(1)$. The primitive sufficient conditions in \citet[Proposition~1]{AI2012} include $m^r/N=O(1)$ for some $r>d$, requiring a control pool that grows much faster than the treated sample. By contrast, when $m\asymp N\asymp n$ and the covariates have bounded full-dimensional densities, their argument gives an $N^{-1/d}$ lower bound for average absolute distance with probability tending to one. Hence the absolute-distance condition cannot hold in dimensions two and above.
By contrast, our theory controls the imbalance $B_{n,\sigma_n^\star}(g)$ directly instead of establishing the distance bound, assuming quadratic objective and smoothness conditions.
\vspace{-0.5em}
\paragraph{Related balancing and distributional methods.}
Explicit balancing offers another route to small aggregate error. \citet{Kallus2020} formulate matching in terms of function-class imbalance and variance, while \citet{WZ2023} establish root-$n$ and efficiency results for matching subject to aggregate balance constraints. Their asymptotic analysis allows replacement and directly imposes aggregate balance constraints. Our procedure concerns minimizing the matching objective only, without any other constraints. The uniform balance is a consequence of minimizing geometric cost.
Our inference procedures also connect to a literature for counterfactual distribution and quantile effects, including regression-based~\citep{CFM2013}, inverse-propensity weighting-based~\citep{DH2014}, and matching-with-replacement based~\citep{YZ2023} procedures. We provide a justification for paired bootstrap after inexact matching without replacement.
\section*{Discussion}
In this paper, we study when inexact one-to-one matching without replacement, on its
own, yields balance sufficient for valid inference. For the quadratic objective
with $d \le 3$ continuous covariates and abundant local control samples, the empirical
imbalance between the treated and matched control samples is $o_P(n^{-1/2})$
uniformly over smooth functions, even though individual pairs may remain at distance
of order $n^{-1/d}$.
As a consequence, a single outcome-blind matched design supports Wald-type and bootstrap inference for the ATT, distributional and quantile treatment effects, and paired randomization tests, without the need to tailor the design to the target or to fit outcome models. We also show that the dimension restriction is sharp: for $d = 4$ the imbalance of a fixed smooth function is non-negligible at the root-$n$ scale,
and for $d > 4$ it dominates the root-$n$ scale, even under the same favorable conditions.
The practical implication of this restriction is that we may rely on matching alone when few continuous covariates are matched inexactly.
When there are discrete features, the matching could be performed exactly in these discrete features while being inexact in the continuous ones.
Otherwise, our analysis suggests that one may need to reduce the dimension of the
matching variables or combine matching with other adjustments such as outcome regression.
Finally, our analysis relies on the quadratic structure of the objective, and whether the same results hold for the ordinary Mahalanobis distance, as conjectured by \citet{GR2023}, remains open as future work.
\section*{Acknowledgments}
Generative AI tools were used for language editing, code assistance, and checking mathematical derivations and numerical implementations. All theoretical arguments, numerical results, and final manuscript content were reviewed and verified by the author, who takes full responsibility for the work.
\bibliographystyle{apalike}
\bibliography{references}
\newpage