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.
84,731 characters
Tuning-Free Tests for Treatment Effect Heterogeneity via Expectile Processes
\begin{bibunit}
\begin{spacing}{1.2}
\title{Tuning-Free Tests for Treatment Effect Heterogeneity\\
via Expectile Processes}
\author{Haokun Lu\thanks{Department of Business Statistics and Econometrics,
Guanghua School of Management, Peking University, Beijing, 100871, China. Email: \texttt{[email removed]}.}\\Peking University
\and Xiaojun Song\thanks{Department of Business Statistics and Econometrics, Guanghua School of Management, Peking University, Beijing, 100871, China. Email: \texttt{[email removed]}. This work was supported by the National Natural Science Foundation of China [Grant Numbers 72373007 and 72333001]. The author also gratefully acknowledges the research support from the Center for Statistical Science of Peking University, China, and the Key Laboratory of Mathematical Economics and Quantitative Finance (Peking University) of the Ministry of Education, China.}\\Peking University}
\maketitle
\vspace{-0.5cm}
\begin{abstract}
This paper develops tuning-free tests for treatment-effect heterogeneity. By characterizing the null hypothesis as constancy of the expectile treatment-effect curve, we construct test statistics as continuous functionals of the mean-centered sample expectile process. We establish a uniform Bahadur representation, derive the asymptotic properties of the test statistics under the null and fixed alternatives, and show nontrivial power against local alternatives approaching the null at the $n^{-1/2}$ rate. For feasible inference, we develop a fast multiplier bootstrap that consistently approximates the nonpivotal null distribution while capturing the first-order effect of estimating the unknown common shift. The procedure is tuning-free, requiring no density or score-function estimation, bandwidth selection, or repeated refitting. The framework further incorporates a pointwise-studentized refinement and a covariate adjustment that reduces asymptotic variation without altering the population discrepancy. Monte Carlo comparisons with existing methods show generally reliable size control, competitive power, and substantial computational savings. An application to a large-scale water-conservation experiment illustrates the proposed tests.
\par\vspace{0.5em}
\noindent\textbf{Keywords}: expectile processes;
treatment-effect heterogeneity;
multiplier bootstrap; studentization; covariate adjustment.
\noindent\textbf{JEL Classification Number:} C12; C14; C15.
\end{abstract}
\end{spacing}
\thispagestyle{empty}
\newpage
\normalsize
\section{Introduction}
\label{sec:intro}
Over the past several decades, the study of treatment-effect heterogeneity has become a central concern in causal inference and program evaluation \citep[see, e.g.,][]{heckman1997making,manski2004statistical,crump2008nonparametric,athey2016recursive,wager2018estimation}. An important motivation is that average treatment effects alone cannot fully characterize how an intervention changes the distribution of outcomes. Indeed, interventions with similar average effects may generate markedly different changes in dispersion and tail behavior. Such distributional differences can reveal economically meaningful heterogeneity in treatment responses, with potentially important implications for inequality and welfare \citep{abadie2002instrumental,bitler2006whatmean}. In this paper, we investigate treatment-effect heterogeneity from this distributional perspective. Let $F_d$ denote the distribution function of $Y(d)$, $d\in\{0,1\}$. We consider the null hypothesis
\begin{equation}
\label{eq:H0}
H_0:\quad
\exists\,\delta\in\mathbb{R}
\ \text{ such that }
F_1(y)=F_0(y-\delta)
\quad\text{for all }y\in\mathbb{R}.
\end{equation}
The null hypothesis in \eqref{eq:H0} allows the marginal outcome distributions to differ only by a common location shift, which necessarily equals the average treatment effect. If $\delta$ were known, testing the null would reduce to a classical two-sample problem after aligning the distributions, for which Kolmogorov--Smirnov, Cram\'er--von Mises, and permutation tests are readily available \citep{lehmann2005testing}. In practice, however, estimating $\delta$ alters the first-order behavior of the empirical CDF and quantile processes commonly used in this literature. Consequently, conventional critical values or naive permutation procedures that ignore this estimation effect can lead to invalid inference, even asymptotically, a consequence of the classical Durbin problem \citep{durbin1973weak,koenker2002inference}. The literature has developed two broad approaches to valid inference. The first restores asymptotic pivotality by removing the effect of nuisance-parameter estimation, notably through the martingale transformations of \citet{khmaladze1981martingale,khmaladze1993goodness}, which \citet{koenker2002inference} apply to quantile-regression-process inference, including tests of location-shift restrictions. The second uses resampling to obtain valid critical values, as in the recentered subsampling approach of \citet{chernozhukov2005inference} and the confidence-set randomization procedure of \citet{ding2016randomization}, which maximizes Fisher randomization $p$-values over candidate shifts. More recently, \citet{chung2021permutation,chung2025quantile} combine these two approaches, drawing on earlier work on permutation tests based on suitably pivotal statistics \citep{neuhaus1993conditional,janssen1997studentized,chung2013exact}. They apply Khmaladze transformations to CDF and quantile processes, respectively, and establish asymptotically valid permutation inference using the resulting pivotal statistics. These approaches provide well-established solutions to the inferential problem, but involve different practical considerations. Khmaladze-based procedures require nonparametric estimation of density and score-related functions, with bandwidth choices that can affect finite-sample performance, as documented by \citet{koenker2002inference}. Recentered subsampling requires selecting a subsample size and repeated estimation, which can be computationally demanding, while confidence-set randomization involves repeated testing over candidate shifts. These considerations leave room for a more direct and computationally efficient approach to valid inference, without the need for auxiliary tuning parameters.
In this paper, we develop such a procedure using expectiles, introduced by \citet{newey1987asymmetric}. Just as quantiles generalize the median through asymmetric absolute loss, expectiles generalize the mean through asymmetric squared loss, assigning different weights to positive and negative deviations. Like quantiles, the full expectile curve characterizes the underlying distribution and is equivariant to location shifts. Consequently, the null hypothesis in \eqref{eq:H0} is equivalent to constancy of the difference between the treatment and control expectile curves, which we call the expectile treatment-effect (ETE) curve. Since the expectile at $\tau=1/2$ coincides with the mean, the common shift under $H_0$ equals the ETE at this index. We therefore center the ETE curve at its mean coordinate to remove the unknown shift from the population restriction. A further advantage for inference comes from the expectile estimating equation. Whereas the derivative of the quantile estimating equation involves the unknown density at the corresponding quantile, the expectile counterpart depends only on the distribution function evaluated at the expectile \citep{newey1987asymmetric,holzmann2016expectile}. Under continuity at the relevant expectiles, this derivative is well defined and uniformly bounded away from zero on any fixed interior index interval. Consequently, all quantities entering the first-order representation can be estimated using sample expectiles and simple sample averages of weighted residuals and indicator functions, without auxiliary density or score-function estimation, nonparametric smoothing, or bandwidth selection. This makes tuning-free inference feasible.
Yet the density-free first-order structure of expectiles does not, by itself, resolve the Durbin problem. Although centering the ETE curve at its mean coordinate removes the unknown common shift from the population discrepancy, estimating that shift still affects the sample process at first order. We establish a uniform Bahadur representation for the centered sample ETE process that explicitly accounts for this estimation effect and show that the process converges weakly under the null to a generally nonpivotal Gaussian limit. Rather than transforming the process to achieve asymptotic pivotality, we retain the estimation effect and use the uniform representation to directly approximate the nonpivotal null distribution through a fast multiplier bootstrap. The influence-function contributions are computed once from the original sample, after which each bootstrap replication requires only random reweighting, without resampling observations or refitting expectile curves. We establish the asymptotic validity of the bootstrap and the consistency of the resulting Kolmogorov--Smirnov and Cram\'er--von Mises tests against fixed departures from constancy. We also characterize the tests' behavior under local alternatives approaching the null at the $n^{-1/2}$ rate, showing how the resulting deterministic drift in the limiting Gaussian process yields nontrivial local power.
We further develop two refinements to improve the practical usefulness of the proposed tests. First, we introduce pointwise studentization to account for differences in sampling variability across expectile indices, which may improve power when sampling precision varies substantially along the expectile curve, as can occur with skewed or heavy-tailed outcomes. Second, since empirical researchers often have access to rich baseline information, we introduce a covariate adjustment that exploits pretreatment covariates to reduce sampling variation without altering the population discrepancy. Under the null, the adjustment reduces the asymptotic covariance kernel in the positive-semidefinite order, without requiring a correctly specified outcome regression model. Both refinements preserve the tuning-free nature of the procedure and can be accommodated within the same fast multiplier bootstrap framework.
The remainder of the paper is organized as follows. Section~\ref{sec:framework} introduces the testing problem and the centered expectile process. Section~\ref{sec:asymptotic} establishes the asymptotic null distribution and investigates the power properties of the proposed tests against fixed and local alternatives. Section~\ref{sec:inference} develops the multiplier bootstrap procedure and pointwise-studentized tests, followed by the covariate adjustment in Section~\ref{sec:covariate-adjustment}. Section~\ref{sec:simu} presents Monte Carlo experiments evaluating size, power, and computational performance relative to existing methods across a broad range of data-generating processes. Section~\ref{sec:empirical} applies the proposed tests to the large-scale water-conservation experiment of \citet{ferraro2013using} to examine treatment-effect heterogeneity under different interventions. Section~\ref{sec:con} concludes. Proofs and additional simulation results are collected in the Online Appendix. Finally, an R package implementing the proposed tests is publicly available at \url{https://github.com/HaokunLU722/embt}.
\textsc{Notation}. For an index set $\mathcal A$, let $\ell^\infty(\mathcal A)$ denote the Banach space of bounded real-valued functions on $\mathcal A$, equipped with the supremum norm $\|f\|_\infty=\sup_{a\in\mathcal A}|f(a)|$. We write $X_n\rightsquigarrow X$ for weak convergence of random elements in $\ell^\infty(\mathcal A)$ in the Hoffmann--J{\o}rgensen sense; see, for example, Definition 1.3.3 of \citet{vdvw2023}. We write $\xrightarrow{p}$ for convergence in probability. The norm $\|\cdot\|$ is Euclidean for finite-dimensional vectors, and stochastic-order statements for random functions are understood with respect to the supremum norm.
\section{Test Formulation}
\label{sec:framework}
\subsection{Setup}
\label{subsec:framework-null}
Consider a randomized experiment with binary treatment indicator $D_i\in\{0,1\}$ and potential outcomes $Y_i(1)$ and $Y_i(0)$. The observed outcome is $Y_i=D_iY_i(1)+(1-D_i)Y_i(0)$. Let $F_d$ denote the marginal distribution function of $Y_i(d)$, $d\in\{0,1\}$. Under random assignment independent of the potential outcomes, these distributions are identified from the observed outcomes in the respective treatment arms. Our question is whether treatment changes the outcome distribution in ways that cannot be explained by a common shift, as formulated by $H_0$ in \eqref{eq:H0}. We write $H_1$ for its complement.
To relate this question to treatment-effect heterogeneity at the individual level, consider the stronger hypothesis that treatment has the same effect on every unit:
\begin{equation}
\label{eq:H0S}
H_0^S:\qquad Y_i(1)-Y_i(0)=\delta\quad\text{for all }i\quad\text{for some }\delta\in\mathbb R.
\end{equation}
Since only one potential outcome is observed for each unit, individual treatment effects are not identified from the observed data, reflecting the fundamental problem of causal inference \citep{holland1986statistics}. For a prespecified $\delta=\delta_0$, \eqref{eq:H0S} is sharp in Fisher's sense; when $\delta$ is unknown, randomization inference must account for the unknown common effect, as in \citet{ding2016randomization}. The stronger null $H_0^S$ implies $H_0$, but the converse need not hold, as individual treatment effects may vary even when the marginal outcome distributions differ only by a common shift. Our tests thus examine treatment-effect heterogeneity through its implications for the outcome distributions.
For the asymptotic analysis, we adopt the following independent two-sample superpopulation framework:
\begin{equation}\label{eq:samples}
Y_1^{(1)},\ldots,Y_{n_1}^{(1)} \stackrel{\mathrm{iid}}{\sim}F_1,
\qquad
Y_1^{(0)},\ldots,Y_{n_0}^{(0)} \stackrel{\mathrm{iid}}{\sim}F_0,
\end{equation}
where the two samples are mutually independent.
\subsection{Expectile Formulation of the Null}
\label{subsec:expectiles}
We now reformulate the null hypothesis \eqref{eq:H0} in terms of expectiles. To begin, we briefly review the expectile functional and its interpretation. For $\tau\in(0,1)$ and a random variable $Y$ with finite second moment, the $\tau$-expectile is defined as the unique minimizer of the asymmetric least-squares criterion
\begin{equation}\label{eq:expectile-loss}
e(\tau)
=
\operatorname*{arg\,min}_{e\in\mathbb{R}}
\mathbb{E}\!\left[
\left|\tau-\mathbf{1}\{Y\le e\}\right|(Y-e)^2
\right].
\end{equation}
The index $\tau$ controls the relative weight placed on positive and negative residuals: when $\tau>1/2$, observations above $e$ receive greater weight, whereas when $\tau<1/2$, observations below $e$ receive greater weight. At $\tau=1/2$, the criterion reduces, up to a constant factor, to ordinary squared loss, so that $e(1/2)=\mathbb{E} Y$ \citep{newey1987asymmetric}. Equivalently, $e(\tau)$ is characterized by the first-order condition
\begin{equation}\label{eq:psi-def}
\mathbb{E}\psi_\tau\{Y,e(\tau)\}=0,
\qquad
\psi_\tau(y,e)
=
2\lvert\tau-\mathbf{1}\{y\le e\}\rvert(y-e).
\end{equation}
Expectiles are translation equivariant: for any $b\in\mathbb{R}$, $e_{Y+b}(\tau)=e_Y(\tau)+b$.
Let $e_d(\tau)$ denote the $\tau$-expectile of $Y(d)$, $d\in\{0,1\}$, and define the \emph{expectile treatment-effect} (ETE) curve
\begin{equation}\label{eq:ETE-def}
\Delta_E(\tau)
=
e_1(\tau)-e_0(\tau),
\qquad
\tau\in(0,1),
\end{equation}
following the terminology of \citet{stahlschmidt2014expectile}. The following result gives the expectile characterization of the location-shift null.
\begin{proposition}
\label{prop:full-curve}
Suppose that $F_0$ and $F_1$ have finite first moments. Then $H_0$ holds if and only if the ETE curve is constant on $(0,1)$; equivalently, there exists $\delta\in\mathbb{R}$ such that $\Delta_E(\tau)=\delta$ for every $\tau\in(0,1)$.
\end{proposition}
Let $\mu_d=\mathbb{E}\{Y(d)\}$. Since $e_d(1/2)=\mu_d$, the constant ETE under $H_0$ equals $\Delta_E(1/2)=\mu_1-\mu_0$. We therefore center the ETE curve at its $1/2$ coordinate and define
\begin{equation}\label{eq:r-def}
r(\tau)
=
\Delta_E(\tau)-\Delta_E(1/2)
=
e_1(\tau)-e_0(\tau)-(\mu_1-\mu_0).
\end{equation}
Under $H_0$, $r(\tau)=0$ for every $\tau\in(0,1)$. This motivates a process-based test of the null. We work with a fixed closed interior interval $\mathcal T=[\tau_L,\tau_U]\subset(0,1)$ containing $1/2$, and base inference on global functionals of a suitably normalized sample analogue of $r(\tau)$ over $\mathcal T$. The next subsection defines the sample process and the corresponding test statistics.
\subsection{The Expectile Process and Test Statistics}
\label{subsec:sample-process}
For $d\in\{0,1\}$, define the sample $\tau$-expectile $\widehat e_d(\tau)$ as
\begin{equation}\label{eq:sample-expectile}
\widehat e_d(\tau)
=
\operatorname*{arg\,min}_{e\in\mathbb{R}}
\frac{1}{n_d}\sum_{i=1}^{n_d}
\left|\tau-\mathbf{1}\{Y_i^{(d)}\le e\}\right|
\left(Y_i^{(d)}-e\right)^2.
\end{equation}
Equivalently, $\widehat e_d(\tau)$ is the unique solution to ${n_d}^{-1}\sum_{i=1}^{n_d} \psi_\tau\{Y_i^{(d)},\widehat e_d(\tau)\}=0.$ The sample ETE curve is then $\widehat\Delta_E(\tau) = \widehat e_1(\tau)-\widehat e_0(\tau)$. Let $\bar Y_d={n_d}^{-1}\sum_{i=1}^{n_d}Y_i^{(d)}$. Since $\widehat e_d(1/2)=\bar Y_d$, the sample analogue of the centered ETE curve is
\[
\widehat r(\tau)
=
\widehat\Delta_E(\tau)-\widehat\Delta_E(1/2)
=
\{\widehat e_1(\tau)-\widehat e_0(\tau)\}
-(\bar Y_1-\bar Y_0).
\]
Write $n=n_1+n_0$, $\lambda_n=n_1/n$, and $N=n_1n_0/n$. Using the usual two-sample normalization, define the centered ETE process
\begin{equation}\label{eq:S-def}
S_n(\tau)
=
\sqrt{N}\,\widehat r(\tau),
\qquad \tau\in\mathcal{T}.
\end{equation}
By construction, $S_n(1/2)\equiv0$ in every finite sample.
To measure departures of the centered ETE curve from zero over $\mathcal{T}$, we use the Kolmogorov--Smirnov and Cram\'er--von Mises functionals
\begin{equation}\label{eq:test-stats}
KS_n
=
\sup_{\tau\in\mathcal{T}}|S_n(\tau)|,
\qquad
CvM_n
=
\int_{\mathcal{T}}S_n^2(\tau)\,d\tau.
\end{equation}
Both tests reject $H_0$ for sufficiently large values of their respective statistics. The next section develops the asymptotic theory underlying these statistics.
\section{Asymptotic Theory}
\label{sec:asymptotic}
For arm $d\in\{0,1\}$, let $\mathbb{E}_d$ denote expectation under $F_d$, let $\mathbb{E}_{n,d}f=n_d^{-1}\sum_{i=1}^{n_d}f(Y_i^{(d)})$, and define $\mathbb G_{n,d}f=\sqrt{n_d}(\mathbb{E}_{n,d}-\mathbb{E}_d)f$. We begin with two regularity conditions governing the sampling scheme and the marginal outcome distributions.
\begin{assumption}\label{ass:sampling}
The samples in \eqref{eq:samples} are independent, observations are i.i.d. within each arm, $n_1\to\infty$, $n_0\to\infty$, and $\lambda_n\to\lambda\in(0,1)$.
\end{assumption}
\begin{assumption}\label{ass:regularity}
For each $d\in\{0,1\}$, $\mathbb{E}_d Y^2<\infty$. The CDF $F_d$ is continuous on an open neighborhood of $\mathcal E_d(\mathcal{T})=\{e_d(\tau):\tau\in\mathcal{T}\}$.
\end{assumption}
Assumption~\ref{ass:sampling} is standard for two-sample asymptotics. Assumption~\ref{ass:regularity} reflects the main regularity trade-off of expectile-based inference. Expectiles are moment-based functionals: a finite first moment is required for their definition, and the second-moment condition here supports root-$n$ Gaussian inference. This is the price paid relative to quantile methods, whose underlying distributional functionals remain well defined without such moment restrictions and are consequently more robust to heavy tails; see \citet{newey1987asymmetric}. Beyond this moment requirement, however, we impose only local continuity of the CDF on a neighborhood of the relevant expectile range, which is sufficient for the uniform first-order linearization used below; see \citet{holzmann2016expectile}. No density is required to exist, let alone to be estimated. This feature, already emphasized by \citet{newey1987asymmetric}, makes expectile-based inference particularly convenient for empirical implementation: the first-order covariance structure can be directly estimated from sample moments. This contrasts with existing tests of treatment-effect heterogeneity. The CDF-based procedure of \citet{chung2021permutation} relies on smooth absolute-continuity conditions and feasible estimation of density and score objects for its martingale transformation. The quantile-based procedure of \citet{chung2025quantile} likewise requires positive densities at the relevant quantiles and uniformly consistent estimation of density and score functions; see also \citet{koenker2002inference} for the underlying quantile-process theory. Thus, the expectile formulation trades stronger moment requirements for a tuning-free inference, i.e., it requires no smoothing-based nonparametric nuisance estimation, such as density, sparsity, or score estimation, and hence no associated tuning-parameter choice.
The absence of density terms is apparent from the Jacobian of the expectile score. For each arm, recalling the score $\psi_\tau$ defined in \eqref{eq:psi-def}, we have
\begin{equation}
\label{eq:score_jacobian}
D_d(\tau)
=
-\left.
\frac{\partial}{\partial e}
\mathbb{E}_d\psi_\tau(Y,e)
\right|_{e=e_d(\tau)}
=
2\left[\tau+(1-2\tau)F_d\{e_d(\tau)\}\right]=2\mathbb{E}_d\lvert\tau-\mathbf{1}\{Y\le e_d(\tau)\}\rvert.
\end{equation}
The influence function of the $\tau$-expectile is then
\begin{equation}\label{eq:phi-def}
\phi_{d,\tau}(y)
=
\frac{\psi_\tau\{y,e_d(\tau)\}}{D_d(\tau)},
\qquad
\mathbb{E}_d\phi_{d,\tau}=0.
\end{equation}
Since $\mathcal{T}=[\tau_L,\tau_U]$ is a fixed interior interval, $\inf_{\tau\in\mathcal{T}}D_d(\tau) \ge 2\min\{\tau_L,1-\tau_U\} >0$, so the denominator in \eqref{eq:phi-def} is uniformly bounded away from zero and $\phi_{d,\tau}$ is well defined uniformly over $\mathcal{T}$. With these definitions in hand, the following lemma provides the uniform Bahadur representation of the sample expectile curve that underlies the asymptotic analysis of the centered ETE process.
\begin{lemma}\label{lem:bahadur}
Under Assumption~\ref{ass:regularity}, for each $d\in\{0,1\}$,
\begin{equation}\label{eq:bahadur}
\widehat e_d(\tau)-e_d(\tau)
=
\mathbb{E}_{n,d}\phi_{d,\tau}
+
R_{d,n_d}(\tau),
\qquad
\sup_{\tau\in\mathcal{T}}|R_{d,n_d}(\tau)|
=
o_p(n_d^{-1/2}).
\end{equation}
Moreover, $R_{d,n_d}(1/2)=0$, and $\{\phi_{d,\tau}:\tau\in\mathcal{T}\}$ is Donsker under $F_d$ with an envelope bounded by $C(1+|Y|)$ for some constant $C$.
\end{lemma}
The $o_p(n_d^{-1/2})$ remainder in Lemma~\ref{lem:bahadur} is sufficient for the first-order asymptotics of $S_n(\tau)$ developed below. Moreover, if $F_d$ additionally admits a density bounded on a neighborhood of $\mathcal E_d(\mathcal{T})$, the proof of Lemma~\ref{lem:bahadur} further shows that the remainder sharpens to $O_p(n_d^{-1})$ uniformly over $\mathcal{T}$. We next combine these arm-specific expansions to obtain the null limit of the centered ETE process.
\subsection{Null Asymptotics}
\label{subsec:asymptotic-limit}
Since $\bar Y_d-\mu_d=\mathbb{E}_{n,d}(Y-\mu_d)$, the influence function of the mean-centered expectile $e_d(\tau)-\mu_d$ is
\begin{equation}\label{eq:xi-def}
\xi_{d,\tau}(y)
=
\phi_{d,\tau}(y)-(y-\mu_d).
\end{equation}
In particular, $\mathbb{E}_d\xi_{d,\tau}=0$ and $\xi_{d,1/2}\equiv0$.
\begin{theorem}
\label{thm:centered-limit}
Under Assumptions~\ref{ass:sampling}--\ref{ass:regularity}, uniformly in $\tau\in\mathcal{T}$,
\begin{equation}\label{eq:S-linear}
S_n(\tau)
=
\sqrt{N}r(\tau)
+
\sqrt{1-\lambda_n}\,\mathbb G_{n,1}\xi_{1,\tau}
-
\sqrt{\lambda_n}\,\mathbb G_{n,0}\xi_{0,\tau}
+o_p(1).
\end{equation}
Under $H_0$,
\begin{equation}\label{eq:S-null-limit}
S_n\rightsquigarrow\mathbb G_{\xi}
\qquad\text{in }\ell^\infty(\mathcal{T}),
\end{equation}
where $\mathbb G_{\xi}$ is a centered Gaussian process with covariance kernel $\Gamma(\tau,\tau') = \mathbb{E}_d\{\xi_{d,\tau}\xi_{d,\tau'}\}$, $d\in\{0,1\}$.
\end{theorem}
Under $H_0$, translation equivariance makes the two arm-specific influence-function laws identical, so the covariance contributions weighted by $1-\lambda$ and $\lambda$ sum to the common kernel $\Gamma$. Hence the null limit is independent of the limiting treatment fraction. By the continuous mapping theorem, we obtain the asymptotic null distributions of the proposed test statistics.
\begin{corollary}\label{cor:functionals}
Under the conditions of Theorem~\ref{thm:centered-limit} and $H_0$,
\begin{equation*}
KS_n
\rightsquigarrow
\sup_{\tau\in\mathcal{T}}|\mathbb G_{\xi}(\tau)|,\quad CvM_n
\rightsquigarrow
\int_\mathcal{T}\mathbb G_{\xi}^2(\tau)\,d\tau.
\end{equation*}
\end{corollary}
\begin{remark}[The Durbin problem]
The mean centering removes the unknown location shift from the population restriction, but not the first-order effect of estimating that shift. Under $H_0$, the common shift $\delta=\mu_1-\mu_0$ is estimated by $\bar Y_1-\bar Y_0$, giving rise to the terms $-(y-\mu_d)$ in \eqref{eq:xi-def}. If $\delta$ were known and the sample ETE curve were instead centered at this known value, these terms would disappear. This estimation effect should be distinguished from nonpivotality: even with a known shift, the limiting Gaussian process would generally have a covariance kernel that depends on the underlying outcome distribution. Rather than removing the estimation effect through a transformation, Section~\ref{sec:inference} develops a multiplier bootstrap that directly approximates the nonpivotal null distribution while retaining the full contribution from estimating $\delta$.
\end{remark}
\subsection{Asymptotic Power}
\label{subsec:asymptotic-power}
We now investigate the power properties of the centered-process tests, first under fixed alternatives and then under local departures from the null. Under fixed alternatives, the uniform representation in Theorem~\ref{thm:centered-limit} implies that the rescaled centered process consistently recovers the population discrepancy $r$. The resulting limits for the two global statistics are given below.
\begin{proposition}
\label{prop:fixed}
Suppose Assumptions~\ref{ass:sampling}--\ref{ass:regularity} hold. Then
\[
\frac{KS_n}{\sqrt N}
\xrightarrow p
\sup_{\tau\in\mathcal{T}}|r(\tau)|,
\qquad
\frac{CvM_n}{N}
\xrightarrow p
\int_\mathcal{T} r^2(\tau)\,d\tau.
\]
In particular, if $r\not\equiv0$ on $\mathcal{T}$, both statistics diverge in probability.
\end{proposition}
Thus, the tests are consistent against any fixed alternative in $H_1$ for which $r\not\equiv0$ on $\mathcal{T}$. We next turn to Pitman-type local alternatives that approach the null at the parametric rate $N^{-1/2}$. Let $\varepsilon_0$ and $\varepsilon_1$ be independent copies of a mean-zero random variable $\varepsilon\sim F$. Under the location-shift null, a convenient representation of the two marginal outcomes is
\[
Y^{(0)}=\mu+\varepsilon_0,
\qquad
Y^{(1)}=\mu+\delta+\varepsilon_1.
\]
To generate local departures from the null while holding the shift $\delta$ fixed, we perturb the treated innovation. Let $b:\mathbb{R}\to\mathbb{R}$ be a nondecreasing and globally Lipschitz perturbation direction and, for $s$ near zero, define the local perturbation map $T_s(x)=x+s b(x)$. We impose $\mathbb{E} b(\varepsilon)=0$, so that $\mathbb{E} T_s(\varepsilon)=0$ for every $s$. This normalization isolates the shape component of the perturbation: any constant component of $b$ would generate only an additional location shift and can therefore be absorbed into $\delta$. We further assume $\mathrm{Var}\{b(\varepsilon)\}>0$ to exclude the trivial perturbation.
For fixed $\kappa\ne0$, set $s_n=\kappa/\sqrt N$ and consider
\begin{equation}\label{eq:local-outcome-path}
H_{1n}(b,\kappa):
\qquad
Y_n^{(0)}\stackrel d=\mu+\varepsilon_0,
\qquad
Y_n^{(1)}\stackrel d=
\mu+\delta+T_{s_n}(\varepsilon_1).
\end{equation}
For all sufficiently large $n$, $T_{s_n}$ is strictly increasing. Hence \eqref{eq:local-outcome-path} departs from the null only through an $N^{-1/2}$ perturbation of the shape of the treated marginal distribution.
To characterize local power, we first determine the first-order response of the expectile curve to this perturbation. Let $e_s(\tau)$ denote the $\tau$-expectile of $T_s(\varepsilon)$, so that $e_0(\tau)=e(\tau)$, and let $D(\tau)$ denote the baseline score Jacobian defined in \eqref{eq:score_jacobian}. Define
\begin{equation}\label{eq:hb-def}
h_b(\tau)
:=
\left.
\frac{\partial e_s(\tau)}{\partial s}
\right|_{s=0}
=
\frac{
2\mathbb{E}\!\left[
\lvert\tau-\mathbf{1}\{\varepsilon\le e(\tau)\}\rvertb(\varepsilon)
\right]
}{
D(\tau)
}.
\end{equation}
Thus $h_b(\tau)$ is the directional derivative of the expectile curve along the perturbation path. The second equality follows by differentiating the population expectile first-order condition with respect to $s$. Equivalently, uniformly over $\mathcal{T}$, $e_s(\tau)=e(\tau)+s h_b(\tau)+o(|s|)$. Under $H_{1n}(b,\kappa)$, let $r_n(\tau):=e_{s_n}(\tau)-e(\tau)$ denote the centered population ETE curve; the equality follows because $\mathbb{E} T_{s_n}(\varepsilon)=0$. The following theorem translates this first-order response into the local drift of the centered ETE process.
\begin{theorem}
\label{thm:local}
Suppose Assumptions~\ref{ass:sampling} and \ref{ass:regularity} hold. Under $H_{1n}(b,\kappa)$,
\begin{equation}\label{eq:path-drift}
\sup_{\tau\in\mathcal{T}}
\left|\sqrt N\,r_n(\tau)-\kappa h_b(\tau)\right|
\to0,
\end{equation}
and
\begin{equation}\label{eq:S-local}
S_n\rightsquigarrow\mathbb G_\xi+\kappa h_b
\qquad\text{in }\ell^\infty(\mathcal{T}),
\end{equation}
where $\mathbb G_\xi$ is the null Gaussian process in Theorem~\ref{thm:centered-limit}. Moreover, $h_b(\tau)<0$ for $\tau<1/2$ and $h_b(\tau)>0$ for $\tau>1/2$.
\end{theorem}
Thus, under the local alternatives, the centered ETE process converges to the null Gaussian process shifted by the deterministic drift $\kappa h_b$. By the continuous mapping theorem,
\[
KS_n
\rightsquigarrow
\sup_{\tau\in\mathcal{T}}
\lvert\mathbb G_\xi(\tau)+\kappa h_b(\tau)\rvert,
\qquad
CvM_n
\rightsquigarrow
\int_\mathcal{T}
\{\mathbb G_\xi(\tau)+\kappa h_b(\tau)\}^2\,d\tau.
\]
Hence, under local alternatives approaching the null at the parametric rate $N^{-1/2}$, both statistics have nondegenerate shifted limits, yielding nontrivial local power against $H_{1n}(b,\kappa)$. The following scale perturbation provides a simple illustration of the resulting drift.
\begin{example}[Local scale perturbation]\label{ex:local-scale-path}
Take $b(x)=x$, so that $T_{s_n}(x)=(1+\kappa/\sqrt N)x$. For all sufficiently large $n$, positive scale and translation equivariance imply that the control and treatment expectiles are $\mu+e(\tau)$ and $\mu+\delta+(1+\kappa/\sqrt N)e(\tau)$, respectively. Since the perturbation leaves the location shift equal to $\delta$, $\sqrt N\,r_n(\tau)=\kappa e(\tau)$, and hence $h_b(\tau)=e(\tau)$. The local drift is therefore proportional to the baseline expectile curve. For $\kappa>0$, the scale perturbation pushes lower expectiles downward and upper expectiles upward while leaving the mean coordinate fixed; for $\kappa<0$, these movements are reversed.
\end{example}
\section{Tuning-Free Inference}
\label{sec:inference}
As discussed above, the null limit in Theorem~\ref{thm:centered-limit} is generally nonpivotal. Nevertheless, the uniform Bahadur representation provides a direct route to feasible inference, as the relevant influence-function components can be estimated from simple sample analogues. To obtain feasible critical values, in this section, we develop a fast, tuning-free multiplier bootstrap that approximates the unknown null distribution through random reweighting of these estimated components, and introduce a pointwise-studentized refinement within the same framework.
\subsection{Multiplier Bootstrap}
\label{subsec:multiplier}
We first construct feasible counterparts of the influence functions in Section~\ref{sec:asymptotic}. For arm $d\in\{0,1\}$, let $\widehat F_d$ denote the empirical CDF and define
\begin{align}
\widehat D_d(\tau)
&=
2\left[
\tau+(1-2\tau)\widehat F_d\{\widehat e_d(\tau)\}
\right],
\label{eq:Dhat}\\
\widehat\xi_{d,i,\tau}
&=
\frac{
\psi_\tau\{Y_i^{(d)},\widehat e_d(\tau)\}
}{
\widehat D_d(\tau)
}
-
\{Y_i^{(d)}-\bar Y_d\}.
\label{eq:xihat}
\end{align}
Thus, $\widehat\xi_{d,i,\tau}$ is the direct plug-in analogue of $\xi_{d,\tau}$ in \eqref{eq:xi-def}. In particular, $\mathbb{E}_{n,d}\widehat\xi_{d,\tau}=0$ for every $\tau\in\mathcal{T}$, and $\widehat\xi_{d,i,1/2}=0$.
For each arm, let $\{V_i^{(d)}\}$ be i.i.d. multipliers, independent across arms and of the data, with mean zero, unit variance, and bounded support. A convenient choice is the two-point distribution of \citet{mammen1993bootstrap}, under which $V=1-a$ with probability $a/\sqrt5$ and $V=a$ with probability $1-a/\sqrt5$, where $a=(1+\sqrt5)/2$. Using the feasible influence values $\widehat\xi_{d,i,\tau}$ in \eqref{eq:xihat}, define the multiplier process
\begin{equation}\label{eq:Sstar-def}
S_n^*(\tau)
=
\sqrt{1-\lambda_n}\frac1{\sqrt{n_1}}
\sum_{i=1}^{n_1}
V_i^{(1)}\widehat\xi_{1,i,\tau}
-
\sqrt{\lambda_n}\frac1{\sqrt{n_0}}
\sum_{j=1}^{n_0}
V_j^{(0)}\widehat\xi_{0,j,\tau},
\qquad \tau\in\mathcal{T}.
\end{equation}
For each draw, construct the bootstrap analogues of the two test statistics,
\begin{equation}\label{eq:bootstrap-stats}
KS_n^*
=
\sup_{\tau\in\mathcal{T}}|S_n^*(\tau)|,
\qquad
CvM_n^*
=
\int_\mathcal{T} \{S_n^*(\tau)\}^2\,d\tau.
\end{equation}
Repeating the multiplier draw $B$ times yields $\{KS_{n,b}^*,CvM_{n,b}^*\}_{b=1}^B$. Let $\widehat c_{KS,1-\alpha}^*$ and $\widehat c_{CvM,1-\alpha}^*$ denote the corresponding empirical $(1-\alpha)$-quantiles. The KS and CvM tests reject $H_0$ when, respectively, $KS_n>\widehat c_{KS,1-\alpha}^*$ and $CvM_n>\widehat c_{CvM,1-\alpha}^*$.
Since all quantities entering $\widehat\xi_{d,i,\tau}$ are computed once from the original sample and held fixed across multiplier draws, the bootstrap calibration reduces to repeatedly reweighting the precomputed influence values in \eqref{eq:Sstar-def}. In particular, no expectile or other nuisance quantity needs to be re-estimated across draws. Write $\rightsquigarrow^*$ in probability for weak convergence under the multiplier law conditional on the original two-sample data $\{Y_i^{(1)}\}_{i=1}^{n_1}\cup\{Y_j^{(0)}\}_{j=1}^{n_0}$. The following theorem establishes the validity of the bootstrap approximation underlying the preceding procedure.
\begin{theorem}
\label{thm:multiplier-fclt}
Suppose Assumptions~\ref{ass:sampling} and \ref{ass:regularity} hold. Under $H_0$ or the sequence of local alternatives $H_{1n}(b,\kappa)$,
\begin{equation*}
S_n^*
\rightsquigarrow^*
\mathbb G_\xi
\qquad
\text{in probability in }\ell^\infty(\mathcal{T}),
\end{equation*}
where $\mathbb G_\xi$ is the centered Gaussian process in Theorem~\ref{thm:centered-limit}.
Under the fixed alternative $H_1$,
\begin{equation*}
S_n^*
\rightsquigarrow^*
\mathbb G_\xi^{A}
\qquad
\text{in probability in }\ell^\infty(\mathcal{T}),
\end{equation*}
where $\mathbb G_\xi^{A}$ is a centered Gaussian process with covariance kernel
\begin{equation*}
\Gamma_A(\tau,\tau')
=
(1-\lambda)\mathbb{E}_1\{\xi_{1,\tau}\xi_{1,\tau'}\}
+
\lambda\mathbb{E}_0\{\xi_{0,\tau}\xi_{0,\tau'}\}.
\end{equation*}
\end{theorem}
Under $H_0$, Theorem~\ref{thm:multiplier-fclt} shows that the multiplier process converges conditionally to the same Gaussian process $\mathbb G_\xi$ as the observed process in Theorem~\ref{thm:centered-limit}. The resulting multiplier critical values therefore yield asymptotically valid KS and CvM tests. Along the local alternatives $H_{1n}(b,\kappa)$, the multiplier process continues to reproduce the centered Gaussian component, while Theorem~\ref{thm:local} shows that the observed process acquires the deterministic drift $\kappa h_b$. Hence, multiplier calibration preserves the nontrivial local-power limits characterized in Section~\ref{sec:asymptotic}.
Under any fixed alternative in $H_1$, the second part of Theorem~\ref{thm:multiplier-fclt} implies that the multiplier statistics, and hence their critical values, remain stochastically bounded, although their limiting covariance kernel generally differs from that under $H_0$. By contrast, Proposition~\ref{prop:fixed} shows that $KS_n$ and $CvM_n$ diverge whenever $r\not\equiv0$ on $\mathcal{T}$. Consequently, the multiplier-bootstrap tests are consistent against any fixed alternative in $H_1$ for which $r\not\equiv0$ on $\mathcal{T}$.
\subsection{A Studentized Refinement}
\label{subsec:studentization}
The asymptotic theory in Section~\ref{sec:asymptotic}, together with the multiplier bootstrap developed in the preceding subsection, already provides feasible inference for the centered-process KS and CvM tests. We now consider an alternative way of weighting information across expectile coordinates. Since the pointwise variability of $S_n(\tau)$ may vary over $\mathcal{T}$, global criteria can place substantial weight on relatively noisy coordinates. A natural approach is therefore to standardize each coordinate by its pointwise standard error, yielding a studentized process. This inverse-variance principle is familiar in process-based testing; for example, \citet{chernozhukov2005inference} recommend inverse pointwise covariance weighting for quantile-process tests and refer to it as Anderson--Darling-type weighting.\footnote{In the classical one-sample goodness-of-fit problem, the Anderson--Darling criterion weights the squared Brownian-bridge process by $[t(1-t)]^{-1}$, the inverse of its pointwise variance. Our studentized CvM statistic follows the same principle, while the studentized KS statistic applies the corresponding standard-error normalization inside the supremum.} The following stylized example illustrates how such reweighting can affect power.
\begin{example}[Variance-dominating noise]\label{ex:variance-dominance}
Let $Z_1,Z_2$ be independent standard normal variables and suppose that the signal $\Delta\ne0$ enters only the first coordinate. Consider
\[
T_\rho=\max\{|Z_1+\Delta|,\rho|Z_2|\},
\qquad
T_S=\max\{|Z_1+\Delta|,|Z_2|\},
\]
where each statistic is calibrated at level $\alpha$ using its own null law ($\Delta=0$). The second coordinate contains no signal, but its variance in $T_\rho$ increases with $\rho^2$. As $\rho\to\infty$, the rejection probability of $T_\rho$ converges to $\alpha$, whereas that of $T_S$ remains strictly above $\alpha$. The same conclusion holds for the quadratic criteria $(Z_1+\Delta)^2+\rho^2Z_2^2$ and $(Z_1+\Delta)^2+Z_2^2$.
\end{example}
Intuitively, the example replaces the continuum of expectile coordinates by two points. In the full process, the same consideration arises when sampling variability differs across $\tau$: standardizing pointwise variability changes how the KS and CvM criteria aggregate signal and noise over $\mathcal{T}$. In the expectile setting, such standardization is particularly convenient because the relevant pointwise variances can be estimated directly from sample moments and empirical CDF values, without nonparametric estimation. There is, however, one complication. By construction, $S_n(1/2)\equiv0$, and its pointwise variance also vanishes at the mean coordinate. Direct studentization of the centered process therefore produces an apparent $0/0$ singularity at $\tau=1/2$. We next show that this singularity can be resolved by an exact factorization of the mean-centered expectile.
Write $c_\tau=2\tau-1$ and define $m_d(\tau)=\mathbb{E}_d|Y-e_d(\tau)|$. The expectile first-order condition
gives $e_d(\tau)-\mu_d=c_\tau m_d(\tau)$. Thus, the mean-centered expectile vanishes at $\tau=1/2$ through the known factor $c_\tau$. The sample expectile likewise satisfies $\widehat e_d(\tau)-\bar Y_d=c_\tau\widehat m_d(\tau)$, where $\widehat m_d(\tau)=\mathbb{E}_{n,d}|Y-\widehat e_d(\tau)|$. Let $q(\tau)=m_1(\tau)-m_0(\tau)$ and $A_n(\tau)=\sqrt N\{\widehat m_1(\tau)-\widehat m_0(\tau)\}$. The same factorization extends to the first-order representation. The following proposition identifies the corresponding influence function.
\begin{proposition}
\label{prop:factorization}
Let $a_d(\tau)=1-2F_d\{e_d(\tau)\}$ and define
\[
\chi_{d,\tau}(y)
=
\frac{
|y-e_d(\tau)|-m_d(\tau)-a_d(\tau)(y-\mu_d)
}{
D_d(\tau)
}.
\]
Under Assumption~\ref{ass:regularity}, uniformly over $\tau\in\mathcal{T}$,
\begin{equation}\label{eq:m-linearization}
\widehat m_d(\tau)-m_d(\tau)
=
\mathbb{E}_{n,d}\chi_{d,\tau}
+
o_p(n_d^{-1/2}).
\end{equation}
For its feasible counterpart, let $\widehat a_d(\tau) =1-2\widehat F_d\{\widehat e_d(\tau)\}$ and define
\[
\widehat\chi_{d,i,\tau}
=
\frac{
|Y_i^{(d)}-\widehat e_d(\tau)|
-\widehat m_d(\tau)
-\widehat a_d(\tau)(Y_i^{(d)}-\bar Y_d)
}{
\widehat D_d(\tau)
}.
\]
Then $\xi_{d,\tau}=c_\tau\chi_{d,\tau}$ and $\widehat\xi_{d,i,\tau} =c_\tau\widehat\chi_{d,i,\tau}$, with $\mathbb{E}_{n,d}\widehat\chi_{d,\tau}=0$. Consequently, $r(\tau)=c_\tau q(\tau)$ and $S_n(\tau)=c_\tau A_n(\tau)$.
\end{proposition}
Proposition~\ref{prop:factorization} shows that the apparent singularity at the mean coordinate is generated entirely by the known factor $c_\tau$. We therefore studentize the regular process $A_n$ rather than divide $S_n$ directly by its vanishing pointwise standard error. Define the asymptotic pointwise variance of $A_n(\tau)$ and its sample analogue by
\begin{equation}\label{eq:student-variance}
v_A^2(\tau)
=
(1-\lambda)\mathbb{E}_1\chi_{1,\tau}^2
+
\lambda\mathbb{E}_0\chi_{0,\tau}^2,
\qquad
\widehat v_{A,n}^2(\tau)
=
(1-\lambda_n)\mathbb{E}_{n,1}\widehat\chi_{1,\tau}^2
+
\lambda_n\mathbb{E}_{n,0}\widehat\chi_{0,\tau}^2.
\end{equation}
The studentized process and the corresponding global statistics are
\begin{equation}\label{eq:U-def}
U_n(\tau)
=
\frac{A_n(\tau)}{\widehat v_{A,n}(\tau)},
\qquad
KS_n^{\mathrm{stud}}
=
\sup_{\tau\in\mathcal{T}}|U_n(\tau)|,
\qquad
CvM_n^{\mathrm{stud}}
=
\int_\mathcal{T} U_n^2(\tau)\,d\tau.
\end{equation}
For $\tau\neq1/2$, $|U_n(\tau)|$ and $U_n^2(\tau)$ coincide with the usual absolute and squared pointwise-standardized versions of $S_n(\tau)$; the factorization extends the same normalization through the mean coordinate.
Under Assumptions~\ref{ass:sampling} and \ref{ass:regularity}, $\widehat v_{A,n}^2$ is uniformly consistent for $v_A^2$, which is finite and continuous on $\mathcal{T}$. For studentization, we additionally require $\inf_{\tau\in\mathcal{T}}v_A^2(\tau)>0$. Under $H_0$, translation equivariance makes the two arm-specific $\chi$ processes identical in distribution, so $v_A^2(\tau)=\mathbb{E}_d\chi_{d,\tau}^2$ for $d\in\{0,1\}$, and the allocation fraction drops out.
\begin{corollary}
\label{cor:studentized}
Suppose Assumptions~\ref{ass:sampling} and \ref{ass:regularity} hold and $\inf_{\tau\in\mathcal{T}}v_A^2(\tau)>0$. Under $H_0$,
\begin{equation}\label{eq:U-null-limit}
U_n
\rightsquigarrow
\mathbb U
:=
\frac{\mathbb G_\chi}{v_A}
\qquad
\text{in }\ell^\infty(\mathcal{T}),
\end{equation}
where $\mathbb G_\chi$ is a centered Gaussian process with covariance kernel $K(\tau,\tau') = \mathbb{E}_d\{\chi_{d,\tau}\chi_{d,\tau'}\}$, $d\in\{0,1\}$. Consequently, $KS_n^{\mathrm{stud}} \rightsquigarrow \sup_{\tau\in\mathcal{T}}|\mathbb U(\tau)|$ and $CvM_n^{\mathrm{stud}} \rightsquigarrow \int_\mathcal{T} \mathbb U^2(\tau)\,d\tau$.
\end{corollary}
For feasible inference, we use the same multipliers as in Section~\ref{subsec:multiplier} and define the studentized multiplier process
\begin{equation}\label{eq:student-multiplier}
U_n^*(\tau)=\frac{1}{\widehat v_{A,n}(\tau)}\left\{\sqrt{1-\lambda_n}\frac{1}{\sqrt{n_1}}\sum_{i=1}^{n_1}V_i^{(1)}\widehat\chi_{1,i,\tau}-\sqrt{\lambda_n}\frac{1}{\sqrt{n_0}}\sum_{j=1}^{n_0}V_j^{(0)}\widehat\chi_{0,j,\tau}\right\},\qquad \tau\in\mathcal{T}.
\end{equation}
Under $H_0$, the same conditional multiplier argument, together with uniform consistency of $\widehat v_{A,n}$, yields $U_n^*\rightsquigarrow^*\mathbb U$ in probability. Applying the KS and CvM functionals to $U_n^*$ therefore provides asymptotically valid bootstrap critical values for the studentized tests. To further see how this reweighting affects power, we return to the local alternatives in Section~\ref{subsec:asymptotic-power}. Recall that $e_s(\tau)$ denotes the $\tau$-expectile of $T_s(\varepsilon)$. Define the corresponding absolute-deviation curve by $m_s(\tau)=\mathbb{E}|T_s(\varepsilon)-e_s(\tau)|$, so that $m_0(\tau)=m(\tau)$. Its first-order response to the perturbation is
\begin{equation}\label{eq:gb-def}
g_b(\tau)
:=
\left.
\frac{\partial m_s(\tau)}{\partial s}
\right|_{s=0}
=
-\frac{
2\mathbb{E}[b(\varepsilon)\mathbf{1}\{\varepsilon\le e(\tau)\}]
}{
D(\tau)
}.
\end{equation}
Along the perturbation path, $e_s(\tau)-\mathbb{E} T_s(\varepsilon)=c_\tau m_s(\tau)$. Differentiating at $s=0$ and using $\mathbb{E} b(\varepsilon)=0$ gives $h_b(\tau)=c_\tau g_b(\tau)$. Consequently, under the local alternatives $H_{1n}(b,\kappa)$,
\[
U_n
\rightsquigarrow
\frac{\mathbb G_\chi+\kappa g_b}{v_A}
\qquad
\text{in }\ell^\infty(\mathcal{T}).
\]
Thus studentization leaves the $N^{-1/2}$ local detection rate unchanged but reweights the deterministic drift across expectile coordinates: the centered process has drift $\kappa c_\tau g_b(\tau)$, whereas the studentized process has drift $\kappa g_b(\tau)/v_A(\tau)$. For the local scale perturbation in Example~\ref{ex:local-scale-path}, $b(x)=x$ and $T_s(x)=(1+s)x$. Positive scale equivariance implies $m_s(\tau)=(1+s)m(\tau)$ for $s$ near zero, and hence $g_b(\tau)=m(\tau)$. The centered and studentized local drifts therefore reduce to $\kappa c_\tau m(\tau)$ and $\kappa m(\tau)/v_A(\tau)$, respectively.
These calculations highlight the potential power gains from studentization. By placing expectile coordinates on a common pointwise noise scale, studentization can improve sensitivity when sampling variability is uneven across $\mathcal{T}$, preventing relatively noisy coordinates from receiving disproportionate weight in the global statistic. Because global power also depends on how the signal aligns with the variance and dependence structure across $\mathcal{T}$, however, studentization need not improve power for every alternative, and no uniform dominance over the raw tests can hold across all data-generating processes. Overall, this section provides tuning-free multiplier inference for the centered-process tests together with a studentized refinement. Section~\ref{sec:simu} examines the finite-sample size and power of both versions across a range of designs.
\section{Covariate Adjustment}
\label{sec:covariate-adjustment}
In modern randomized experiments, researchers often have access to a rich set of pretreatment covariates in addition to treatment assignment and outcomes. Although random assignment balances these characteristics on average, chance imbalances generally remain in realized samples. When pretreatment covariates are predictive of outcomes, such information can be used to absorb otherwise extraneous sampling variation and improve the precision of treatment-effect inference \citep{zhang2008improving,li2020rerandomization}. This argument extends naturally to tests of treatment-effect heterogeneity: removing outcome variation explained by pretreatment covariates can make treatment-effect variation more readily apparent \citep{ding2016randomization}. In this section, we develop a covariate adjustment for the centered expectile process using a fixed, finite-dimensional, prespecified covariate basis.
\begin{assumption}
\label{ass:covariate}
Let $Z=z(X)\in\mathbb R^p$, where $p$ is fixed and $z(\cdot)$ is prespecified. Suppose $\mathbb{E}\|Z\|^2<\infty$ and $\Sigma_Z:=\mathbb{E}[(Z-\mathbb{E} Z)(Z-\mathbb{E} Z)^\top]$ is positive definite.
\end{assumption}
A natural choice of $Z$ is $X$ itself, although $Z$ may also include indicators for categorical variables, low-order transformations, or prespecified interactions. Write $\widetilde Z=Z-\mathbb{E} Z$. Random assignment ensures balance of pretreatment covariates in the population, but their sample means need not coincide in a realized experiment. Writing $\bar Z_d=\mathbb{E}_{n,d}Z$, we summarize this chance imbalance by
\begin{equation}\label{eq:ca-imbalance}
B_n
=
\sqrt N(\bar Z_1-\bar Z_0).
\end{equation}
The idea of covariate adjustment is to use the realized covariate imbalance to account for sampling variation in $S_n(\tau)$ without changing the population discrepancy $r(\tau)$. Under random assignment, $B_n$ is centered and $B_n=O_p(1)$. Hence, for any deterministic $a(\tau)\in\mathbb R^p$, the adjusted process $S_n(\tau)-a(\tau)^\top B_n$ has the same deterministic component $\sqrt N\,r(\tau)$ as $S_n(\tau)$, while its first-order sampling variation may differ. We therefore seek the choice of $a(\tau)$ that minimizes the pointwise asymptotic variance.
This choice follows from the usual linear projection argument. With a slight abuse of notation, let $\mathbb{E}_d$ below denote expectation under the joint distribution of $(Y(d),Z)$. Using the influence representation of $S_n(\tau)$ in \eqref{eq:S-linear}, define the limiting covariance between the covariate imbalance and the first-order sampling error of the centered expectile process by
\[
C_\xi(\tau)
:=
\lim_{n\to\infty}
\mathrm{Cov}\!\left\{
B_n,\,
S_n(\tau)-\sqrt N\,r(\tau)
\right\}
=
(1-\lambda)\mathbb{E}_1\{\widetilde Z\,\xi_{1,\tau}\}
+
\lambda\mathbb{E}_0\{\widetilde Z\,\xi_{0,\tau}\}.
\]
Since $\mathrm{Var}(B_n)\to\Sigma_Z$, the usual least-squares projection formula implies that the variance-minimizing coefficient among adjustments of the form $a(\tau)^\top B_n$ is $\gamma(\tau):=\Sigma_Z^{-1}C_\xi(\tau)$. Thus, $\gamma(\tau)^\top B_n$ is the asymptotic linear projection of the first-order sampling error on the realized covariate imbalance, so that subtracting it removes the corresponding linearly predictable component.
The population projection has a direct sample analogue. Let $\widetilde Z_{d,i}=Z_{d,i}-\bar Z_d$ and define
\begin{align}
\widehat\Sigma_Z
&=
(1-\lambda_n)\mathbb{E}_{n,1}
(\widetilde Z_1\widetilde Z_1^\top)
+
\lambda_n\mathbb{E}_{n,0}
(\widetilde Z_0\widetilde Z_0^\top),
\nonumber\\
\widehat C_\xi(\tau)
&=
(1-\lambda_n)\mathbb{E}_{n,1}
\{\widetilde Z_1\widehat\xi_{1,\tau}\}
+
\lambda_n\mathbb{E}_{n,0}
\{\widetilde Z_0\widehat\xi_{0,\tau}\},
\qquad
\widehat\gamma(\tau)
=
\widehat\Sigma_Z^{-1}\widehat C_\xi(\tau).
\label{eq:ca-feasible}
\end{align}
The covariate-adjusted process is
\begin{equation}\label{eq:ca-process}
S_n^{\mathrm{adj}}(\tau)
=
S_n(\tau)-\widehat\gamma(\tau)^\top B_n,
\qquad \tau\in\mathcal{T},
\end{equation}
with corresponding feasible influence values $\widehat\xi_{d,i,\tau}^{\mathrm{adj}} = \widehat\xi_{d,i,\tau} -\widehat\gamma(\tau)^\top\widetilde Z_{d,i}$. The adjusted KS and CvM statistics are constructed from $S_n^{\mathrm{adj}}$ in the same way as their unadjusted counterparts. The following theorem shows that the adjustment preserves the population discrepancy while reducing the covariance of the expectile process.
\begin{theorem}
\label{thm:covariate-adjustment}
Suppose Assumptions~\ref{ass:sampling}--\ref{ass:covariate} hold. Then, uniformly over $\tau\in\mathcal{T}$,
\begin{equation}\label{eq:ca-linearization}
S_n^{\mathrm{adj}}(\tau)
=
\sqrt N\,r(\tau)
+
\sqrt{1-\lambda_n}\,
\mathbb G_{n,1}\xi_{1,\tau}^{\mathrm{adj}}
-
\sqrt{\lambda_n}\,
\mathbb G_{n,0}\xi_{0,\tau}^{\mathrm{adj}}
+
o_p(1),
\end{equation}
where $\xi_{d,\tau}^{\mathrm{adj}}(y,z) = \xi_{d,\tau}(y) - \gamma(\tau)^\top(z-\mathbb{E} Z)$. Consequently, under $H_0$,
\[
S_n^{\mathrm{adj}}
\rightsquigarrow
\mathbb G_{\mathrm{adj}}
\qquad\text{in }\ell^\infty(\mathcal{T}),
\]
where $\mathbb G_{\mathrm{adj}}$ is a centered Gaussian process with covariance kernel
\begin{equation}\label{eq:ca-kernel}
\Gamma_{\mathrm{adj}}(\tau,s)
=
\Gamma(\tau,s)
-
C_\xi(\tau)^\top\Sigma_Z^{-1}C_\xi(s),
\end{equation}
and $\Gamma$ is the covariance kernel in Theorem~\ref{thm:centered-limit}. In particular, $\Gamma_{\mathrm{adj}}\preceq\Gamma$ in the positive-semidefinite kernel order.
\end{theorem}
The positive-semidefinite ordering in Theorem~\ref{thm:covariate-adjustment} implies that covariate adjustment leaves the population discrepancy $r$ unchanged while reducing asymptotic variation in every finite linear contrast of the process. At any fixed expectile index, this yields a direct local-power ordering for the corresponding Wald-type test: under $N^{-1/2}$ local alternatives, adjustment preserves the local signal while reducing its asymptotic variance, thereby increasing the noncentrality parameter and hence local asymptotic power whenever the covariates are informative. For the global KS and CvM statistics considered here, however, the comparison is more subtle. Covariance ordering implies greater concentration of centered Gaussian limits over common symmetric convex sets, including the supremum- and quadratic-norm acceptance regions underlying KS- and CvM-type statistics \citep{anderson1955integral}. Because the adjusted and unadjusted tests are each calibrated using their own covariance-dependent null distribution, this ordering alone does not imply a uniform power ranking. Related work on CvM-type omnibus tests shows that local power depends critically on the interaction between the alternative direction and the covariance structure of the limiting process; see \citet{neuhaus1976asymptotic}. Nonetheless, covariate adjustment preserves the same deterministic signal while removing covariate-predictable sampling variation, suggesting power gains in many empirically relevant settings. The simulation study below examines these gains. The multiplier bootstrap and studentized refinement from Section~\ref{sec:inference} carry over directly to the adjusted process, and the parallel formulas are omitted for brevity.
\section{Monte Carlo Experiments}
\label{sec:simu}
In this section, we conduct Monte Carlo experiments to assess the finite-sample size and power of the proposed expectile multiplier tests relative to existing confidence-set randomization, recentered subsampling, and martingale-transformed permutation procedures. We further investigate the effects of pointwise studentization and covariate adjustment and compare computational costs across methods. Following the location-scale design of \citet{koenker2002inference}, also adopted in subsequent studies \citep{chernozhukov2005inference,chung2021permutation,chung2025quantile}, we generate potential outcomes according to
\begin{equation}\label{eq:simulation-dgp}
Y_i(0)=\varepsilon_i,
\qquad
Y_i(1)=1+(1+\sigma_\delta)\varepsilon_i.
\end{equation}
The parameter $\sigma_\delta$ controls the degree of heterogeneity, with $\sigma_\delta=0$ corresponding to the null and $\sigma_\delta\in\{0.2,0.5\}$ representing two levels of scale heterogeneity. We draw the innovations $\varepsilon_i$ independently from four distributions: $N(0,1)$, $t_5$, $\mathrm{LN}(0,1)$, and the product-normal distribution $Z_{1i}Z_{2i}$, where $Z_{1i}$ and $Z_{2i}$ are independent standard normal variables. The $t_5$ and lognormal distributions allow us to examine heavy tails and skewness, respectively. In particular, the product-normal distribution has a continuous CDF but a density that diverges to infinity as $x\to0$. This makes it a useful stress test for our procedure, which requires only finite second moments and CDF continuity near the relevant expectiles, without the density regularity conditions often imposed by existing methods.
To examine distributional departures beyond the scale changes in \eqref{eq:simulation-dgp}, we further consider normal-to-$t_5$ alternatives, both with and without variance matching. Draw $U_i\sim U(0,1)$ and set
\begin{equation}\label{eq:simulation-shape-dgp}
Y_i(0)=\Phi^{-1}(U_i),
\qquad
Y_i(1)=1+aF_{t_5}^{-1}(U_i),
\qquad
a\in\left\{\sqrt{\frac{3}{5}},1\right\}.
\end{equation}
The choice $a=\sqrt{3/5}$ matches the marginal variances, isolating differences in distributional shape, whereas $a=1$ allows both shape and scale to differ. In both cases, the marginal means differ by one.
Throughout the Monte Carlo experiments, we consider balanced treatment allocation ($n_1/n=1/2$) with total sample sizes $n\in\{100,200,400,800\}$. We also examine unbalanced allocation ($n_1/n=1/3$) with $n\in\{150,300,600\}$, for which the results are reported in Section \ref{app:additional-monte-carlo} of the Online Appendix. Each design uses 1,000 Monte Carlo replications, and all tests are conducted at the $5\%$ nominal level. We compare the proposed expectile multiplier tests with existing procedures using the following implementations:
\begin{itemize}[leftmargin=1.95cm,labelsep=0.4em,itemsep=0.25em,
topsep=0.4em,parsep=0pt]
\item[\textbf{eMBT}]
The proposed expectile multiplier-bootstrap tests, comprising KS and CvM statistics and their pointwise-studentized counterparts (S--KS and S--CvM). All four variants use 81 equally spaced expectile indices over $[0.10,0.90]$ and 999 Mammen multiplier draws.
\item[\textbf{FRT-CI}]
The Fisher randomization test with a confidence interval proposed by \citet{ding2016randomization}. The procedure maximizes the randomization $p$-value over a $(1-\gamma)$ confidence interval for the unknown treatment effect and adds $\gamma$ to obtain the adjusted $p$-value. We set $\gamma=0.001$ and use 999 permutations.
\item[\textbf{RS}]
The subsampling test of \citet{chernozhukov2005inference}, based on the recentered quantile treatment-effect process. We use 250 subsamples of size $b=\lfloor 20+n^{1/4}\rfloor$ and the quantile grid $\{0.05,0.10,\ldots,0.95\}$, with pointwise standard-error weighting and the authors' finite-sample critical-value adjustment.
\item[\textbf{D-mtPT}]
The CDF-based martingale-transformed permutation test of \citet{chung2021permutation}, implemented using \texttt{PT.Khmaladze.fit} with the default \texttt{akj.control} settings and 999 permutations.
\item[\textbf{Q-mtPT}]
The quantile-based martingale-transformed permutation test of \citet{chung2025quantile}, implemented using \texttt{PTQTE.Khmaladze.fit} with 999 permutations over the quantile grid $\{0.10,0.15,\ldots,0.90\}$.\footnote{Both procedures are implemented using \texttt{RATest} version 0.1.12.}
\end{itemize}
Table~\ref{tab:balanced-size} reports empirical rejection frequencies under the null. The four eMBT variants exhibit stable size control across all four innovation distributions, including the product-normal design with an unbounded density. Although some mild overrejection occurs at $n=100$, particularly under lognormality, all eMBT rejection frequencies lie between 0.038 and 0.065 for $n\ge 200$, with comparable performance for the centered and studentized tests. This stability is particularly encouraging given that the multiplier procedure requires neither density nor score estimation. Among the comparison procedures, FRT-CI tends to be conservative, consistent with its confidence-set construction, while RS exhibits more variation across distributions. D-mtPT and Q-mtPT are reasonably well sized under normal and $t_5$ innovations, with somewhat elevated rejection frequencies under product-normality and more pronounced overrejection under lognormality, especially in smaller samples.
Table~\ref{tab:balanced-power} reports empirical rejection frequencies under the alternatives. Overall, the eMBT tests compare favorably with existing procedures, frequently achieving the highest empirical power. For the scale alternatives in Panels A and B, the proposed tests generally outperform the comparison procedures under normal and $t_5$ innovations, particularly at moderate sample sizes. Under lognormality, eMBT remains competitive, especially with studentization, although D-mtPT and Q-mtPT often exhibit higher rejection frequencies. These comparisons should be interpreted in light of their substantial size distortions under the null (Table~\ref{tab:balanced-size}). The advantage of eMBT is particularly pronounced under product-normal innovations, where FRT-CI exhibits little power against the scale alternatives and D-mtPT has relatively low rejection frequencies, while the proposed tests consistently perform best. Panels C and D extend the comparison to normal-to-$t_5$ alternatives with and without variance matching. Q-mtPT achieves the highest empirical power in the variance-matched design, whereas eMBT is the most powerful across all sample sizes when the alternative also involves a scale difference. These results suggest that expectile-based tests are sensitive to both scale and shape departures, with particularly strong power against scale changes, while quantile-based procedures may be more sensitive to certain variance-preserving shape alternatives.
\begin{table}[htbp]
\centering
\begin{threeparttable}
\caption{Empirical size under balanced assignment}
\label{tab:balanced-size}
\small
\setlength{\tabcolsep}{3.5pt}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}llrrrrrrrr@{}}
\toprule
&&\multicolumn{4}{c}{eMBT}&
\multicolumn{4}{c}{Comparison tests}\\
\cmidrule(lr){3-6}\cmidrule(lr){7-10}
Innovation & $n$ & KS & CvM & S--KS & S--CvM &
FRT-CI & RS & D-mtPT & Q-mtPT\\
\midrule
\multirow{4}{*}{Normal}
&100&0.069&0.068&0.068&0.063&0.015&0.056&0.055&0.042\\
&200&0.048&0.051&0.053&0.051&0.035&0.075&0.039&0.044\\
&400&0.045&0.052&0.052&0.052&0.030&0.088&0.051&0.049\\
&800&0.038&0.043&0.044&0.044&0.040&0.074&0.045&0.041\\
\addlinespace[2pt]
\multirow{4}{*}{$t_5$}
&100&0.048&0.053&0.063&0.056&0.015&0.045&0.043&0.045\\
&200&0.043&0.057&0.058&0.061&0.029&0.041&0.036&0.055\\
&400&0.047&0.055&0.051&0.056&0.026&0.047&0.048&0.039\\
&800&0.045&0.052&0.042&0.050&0.025&0.031&0.050&0.045\\
\addlinespace[2pt]
\multirow{4}{*}{Lognormal}
&100&0.088&0.066&0.088&0.064&0.034&0.027&0.191&0.200\\
&200&0.054&0.053&0.057&0.045&0.038&0.023&0.161&0.136\\
&400&0.061&0.048&0.051&0.047&0.041&0.029&0.122&0.107\\
&800&0.047&0.047&0.050&0.049&0.044&0.039&0.094&0.063\\
\addlinespace[2pt]
\multirow{4}{*}{Product-normal}
&100&0.055&0.052&0.073&0.048&0.032&0.046&0.075&0.064\\
&200&0.059&0.053&0.065&0.053&0.030&0.036&0.073&0.073\\
&400&0.059&0.054&0.057&0.057&0.049&0.039&0.072&0.069\\
&800&0.046&0.048&0.053&0.051&0.048&0.031&0.063&0.073\\
\bottomrule
\end{tabular*}
\begin{tablenotes}[flushleft]
\footnotesize
\item Notes: Entries report rejection frequencies at the $5\%$
nominal level based on 1,000 Monte Carlo replications of
\eqref{eq:simulation-dgp} with $\sigma_\delta=0$ and $n_1=n_0=n/2$.
\end{tablenotes}
\end{threeparttable}
\end{table}
\begin{table}[htbp]
\centering
\begin{threeparttable}
\caption{Empirical power under balanced assignment}
\label{tab:balanced-power}
\small
\setlength{\tabcolsep}{3.5pt}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}llrrrrrrrr@{}}
\toprule
&&\multicolumn{4}{c}{eMBT}&
\multicolumn{4}{c}{Comparison tests}\\
\cmidrule(lr){3-6}\cmidrule(lr){7-10}
Innovation & $n$ & KS & CvM & S--KS & S--CvM &
FRT-CI & RS & D-mtPT & Q-mtPT\\
\midrule
\multicolumn{10}{l}{\textit{Panel A: $\sigma_\delta=0.2$}}\\[2pt]
\multirow{4}{*}{Normal}
&100&0.220&0.239&0.235&0.232&0.049&0.144&0.082&0.125\\
&200&0.359&0.401&0.372&0.393&0.112&0.295&0.159&0.214\\
&400&0.649&0.692&0.662&0.678&0.226&0.511&0.334&0.400\\
&800&0.932&0.957&0.944&0.953&0.535&0.803&0.627&0.703\\
\addlinespace[2pt]
\multirow{4}{*}{$t_5$}
&100&0.130&0.144&0.158&0.158&0.028&0.084&0.085&0.094\\
&200&0.201&0.261&0.264&0.284&0.070&0.147&0.133&0.195\\
&400&0.364&0.462&0.467&0.494&0.133&0.231&0.217&0.358\\
&800&0.616&0.751&0.756&0.794&0.317&0.465&0.445&0.603\\
\addlinespace[2pt]
\multirow{4}{*}{Lognormal}
&100&0.099&0.083&0.103&0.085&0.054&0.050&0.294&0.270\\
&200&0.109&0.113&0.141&0.136&0.118&0.067&0.387&0.349\\
&400&0.128&0.136&0.178&0.161&0.180&0.092&0.452&0.360\\
&800&0.144&0.195&0.281&0.253&0.326&0.185&0.551&0.494\\
\addlinespace[2pt]
\multirow{4}{*}{Product-normal}
&100&0.082&0.096&0.126&0.099&0.040&0.063&0.082&0.091\\
&200&0.134&0.175&0.174&0.179&0.032&0.079&0.070&0.120\\
&400&0.224&0.272&0.275&0.297&0.032&0.120&0.061&0.198\\
&800&0.402&0.510&0.501&0.547&0.038&0.240&0.097&0.311\\
\midrule
\multicolumn{10}{l}{\textit{Panel B: $\sigma_\delta=0.5$}}\\[2pt]
\multirow{4}{*}{Normal}
&100&0.725&0.791&0.756&0.787&0.219&0.535&0.202&0.470\\
&200&0.959&0.983&0.968&0.977&0.577&0.864&0.495&0.801\\
&400&1.000&1.000&0.999&1.000&0.935&0.996&0.815&0.980\\
&800&1.000&1.000&1.000&1.000&1.000&1.000&0.995&1.000\\
\addlinespace[2pt]
\multirow{4}{*}{$t_5$}
&100&0.425&0.544&0.575&0.587&0.128&0.315&0.202&0.375\\
&200&0.722&0.843&0.856&0.879&0.350&0.552&0.410&0.688\\
&400&0.947&0.988&0.989&0.993&0.739&0.882&0.734&0.953\\
&800&0.998&1.000&1.000&1.000&0.985&0.998&0.975&1.000\\
\addlinespace[2pt]
\multirow{4}{*}{Lognormal}
&100&0.178&0.195&0.244&0.230&0.173&0.144&0.560&0.568\\
&200&0.253&0.307&0.407&0.373&0.390&0.246&0.698&0.712\\
&400&0.348&0.448&0.654&0.568&0.681&0.500&0.881&0.899\\
&800&0.575&0.720&0.923&0.856&0.951&0.880&0.991&0.990\\
\addlinespace[2pt]
\multirow{4}{*}{Product-normal}
&100&0.285&0.334&0.386&0.362&0.032&0.197&0.127&0.244\\
&200&0.501&0.618&0.628&0.651&0.037&0.283&0.122&0.425\\
&400&0.783&0.880&0.875&0.904&0.059&0.529&0.158&0.636\\
&800&0.972&0.991&0.992&0.996&0.057&0.851&0.256&0.904\\
\midrule
\multicolumn{10}{l}{\textit{Panel C: variance-matched normal-to-$t_5$ alternative}}\\[2pt]
\multirow{4}{*}{Normal/$t_5$}
&100&0.075&0.091&0.102&0.096&0.038&0.077&0.020&0.107\\
&200&0.072&0.092&0.109&0.108&0.073&0.138&0.030&0.139\\
&400&0.082&0.128&0.143&0.155&0.183&0.231&0.064&0.290\\
&800&0.089&0.174&0.239&0.250&0.370&0.408&0.127&0.490\\
\midrule
\multicolumn{10}{l}{\textit{Panel D: unstandardized normal-to-$t_5$ alternative}}\\[2pt]
\multirow{4}{*}{Normal/$t_5$}
&100&0.166&0.168&0.174&0.157&0.027&0.094&0.056&0.079\\
&200&0.313&0.333&0.307&0.314&0.039&0.152&0.087&0.095\\
&400&0.609&0.636&0.581&0.568&0.065&0.222&0.128&0.150\\
&800&0.904&0.897&0.877&0.855&0.096&0.374&0.257&0.271\\
\bottomrule
\end{tabular*}
\begin{tablenotes}[flushleft]
\footnotesize
\item Notes: Entries report rejection frequencies at the $5\%$
nominal level based on 1,000 Monte Carlo replications under
balanced assignment, using \eqref{eq:simulation-dgp} (Panels A--B)
and \eqref{eq:simulation-shape-dgp} (Panels C--D).
\end{tablenotes}
\end{threeparttable}
\end{table}
Table~\ref{tab:balanced-power} also shows how pointwise studentization affects power. While its impact is modest under normal innovations, studentization often yields substantial gains under $t_5$, product-normal, and lognormal innovations, as well as under the variance-matched normal-to-$t_5$ alternative. Although some power losses occur, particularly in Panel D, they are generally limited relative to the gains observed elsewhere, and the overall results favor studentization. To further illustrate these differences, Figure~\ref{fig:studentization-diagnostic} plots the pointwise standard deviation $v_A(\tau)$ of the regular process $A_n(\tau)$ and the standardized scale signal $s(\tau)=m(\tau)/v_A(\tau)$, where $m(\tau)=\mathbb{E}|\varepsilon-e_\varepsilon(\tau)|$, each normalized by its value at $\tau=1/2$. Under normality, both profiles are relatively flat. Under $t_5$ and product-normal innovations, pointwise variability increases toward the tails while the standardized signal weakens. Under lognormality, the profiles exhibit a pronounced asymmetry: the pointwise standard deviation rises sharply toward the upper expectiles, whereas the standardized scale signal declines over the same range. These patterns indicate that coordinates with greater sampling variability need not carry correspondingly stronger standardized signals. Thus, by placing expectile coordinates on a common noise scale, studentization reduces the relative influence of noisy regions and can enhance the power of global tests against scale departures.
\begin{figure}[t]
\centering
\includegraphics[width=0.9\textwidth]{figures/figure_studentization_diagnostic.pdf}
\caption{Population diagnostics for pointwise studentization.
Panel (a) shows the pointwise standard deviation $v_A(\tau)$
of the regular expectile process, and Panel (b) shows the
standardized scale signal $s(\tau)=m(\tau)/v_A(\tau)$.
Both profiles are normalized by their values at $\tau=1/2$.}
\label{fig:studentization-diagnostic}
\end{figure}
We next examine the finite-sample benefits of covariate adjustment using a design with a binary pretreatment covariate. Specifically, let $X_i\sim\operatorname{Bernoulli}(1/2)$ and $\varepsilon_i\sim N(0,1)$ be independent, and set
\begin{equation}\label{eq:simulation-covariate-dgp}
\sigma(X_i)=
\begin{cases}
1/\sqrt{5},&X_i=0,\\
3/\sqrt{5},&X_i=1,
\end{cases}
\qquad
Y_i(0)=\sigma(X_i)\varepsilon_i,
\qquad
Y_i(1)=1+(1+\gamma)\sigma(X_i)\varepsilon_i.
\end{equation}
The conditional scale is three times larger for $X_i=1$ than for $X_i=0$, while $\mathbb{E}\{\sigma^2(X_i)\}=1$. When $\gamma=0$, the treatment effect is constant at one, so the null holds. Positive values of $\gamma$ introduce a treatment-induced scale change while preserving the average treatment effect of one. Using $z(X)=X$ as the adjustment basis, we compare each eMBT statistic with its covariate-adjusted counterpart under the same simulation settings.
\begin{table}[htbp]
\centering
\begin{threeparttable}
\caption{Comparison of eMBT tests with and without covariate adjustment}
\label{tab:covariate-adjustment}
\small
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lrrrrrrrr@{}}
\toprule
&\multicolumn{4}{c}{Unadjusted eMBT}&
\multicolumn{4}{c}{Covariate-adjusted eMBT}\\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}
$n$ & KS & CvM & S--KS & S--CvM & KS & CvM & S--KS & S--CvM\\
\midrule
\multicolumn{9}{l}{\textit{Panel A: $\gamma=0$ (size)}}\\[2pt]
100&0.071&0.068&0.091&0.068&0.066&0.068&0.086&0.072\\
200&0.050&0.051&0.068&0.049&0.053&0.052&0.065&0.054\\
400&0.056&0.053&0.056&0.048&0.053&0.057&0.060&0.061\\
800&0.047&0.045&0.043&0.042&0.057&0.058&0.050&0.055\\
\midrule
\multicolumn{9}{l}{\textit{Panel B: $\gamma=0.2$}}\\[2pt]
100&0.153&0.163&0.180&0.164&0.173&0.194&0.214&0.202\\
200&0.189&0.231&0.227&0.243&0.219&0.293&0.293&0.315\\
400&0.337&0.405&0.397&0.429&0.380&0.488&0.487&0.534\\
800&0.603&0.687&0.671&0.711&0.689&0.784&0.799&0.834\\
\midrule
\multicolumn{9}{l}{\textit{Panel C: $\gamma=0.5$}}\\[2pt]
100&0.398&0.488&0.496&0.509&0.447&0.548&0.563&0.613\\
200&0.679&0.778&0.771&0.808&0.760&0.875&0.878&0.907\\
400&0.945&0.972&0.974&0.986&0.966&0.995&0.994&0.996\\
800&0.999&1.000&1.000&1.000&1.000&1.000&1.000&1.000\\
\bottomrule
\end{tabular*}
\begin{tablenotes}[flushleft]
\footnotesize
\item Notes: Entries are rejection frequencies at the $5\%$ nominal level
based on 1,000 Monte Carlo replications of
\eqref{eq:simulation-covariate-dgp}, with $n_1=n_0=n/2$. The
adjustment basis is $z(X)=X$. S denotes pointwise studentization.
\end{tablenotes}
\end{threeparttable}
\end{table}
The results are summarized in Table~\ref{tab:covariate-adjustment}. Covariate adjustment has little effect on empirical size but yields clear improvements in power under both alternatives. These gains are evident across all four statistics, with and without pointwise studentization. For example, at $n=800$ and $\gamma=0.2$, the empirical power of CvM increases from 0.687 to 0.784, while that of S--CvM rises from 0.711 to 0.834. Recall that, in this design, the baseline covariate affects the conditional scales of the potential outcomes but not their conditional means. The power gains thus reflect its predictive content for the expectile influence functions rather than the outcome means, consistent with the variance reduction established in Theorem~\ref{thm:covariate-adjustment}.
\begin{table}[htbp]
\centering
\begin{threeparttable}
\caption{Mean single-core computation time (seconds)}
\label{tab:computation-time}
\small
\setlength{\tabcolsep}{9pt}
\begin{tabular}{lrrrr}
\toprule
Method & $n=100$ & $n=200$ & $n=400$ & $n=800$ \\
\midrule
eMBT & 0.010 & 0.017 & 0.031 & 0.060 \\
FRT-CI & 12.178 & 13.873 & 17.098 & 24.206 \\
RS & 0.229 & 0.232 & 0.228 & 0.235 \\
D-mtPT & 0.295 & 0.377 & 0.620 & 1.356 \\
Q-mtPT & 5.008 & 6.740 & 10.263 & 19.164 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item Notes: Entries are mean elapsed seconds over 100 independently generated
datasets for each balanced normal-null design. The eMBT row reports a joint
call computing all four reported eMBT variants. Computations use one core of
an Apple M4 Pro chip with 24 GB of memory.
\end{tablenotes}
\end{threeparttable}
\end{table}
Finally, we compare the computational costs of the proposed tests with those of the competing procedures. The results are reported in Table~\ref{tab:computation-time}. Computing all four eMBT variants jointly takes only 0.010 seconds at $n=100$ and 0.060 seconds at $n=800$, which reflects the computational simplicity of our procedure: it requires no nonparametric nuisance estimation, and the multiplier bootstrap reweights precomputed influence-function quantities without repeated re-estimation. In contrast, the competing procedures involve more computationally intensive steps, including subsampling, repeated permutation calculations, and nonparametric nuisance estimation. At $n=800$, their computation times range from 0.235 seconds for RS to 24.206 seconds for FRT-CI. Taken together, the simulations show that eMBT achieves generally reliable size control and competitive power across distributions with substantially different tail behavior, skewness, and density characteristics. Covariate adjustment can further enhance power by exploiting baseline information without altering the population discrepancy being tested. Together with its tuning-free construction and computationally efficient multiplier bootstrap, these properties make eMBT a practical and flexible tool for empirical studies of treatment-effect heterogeneity.
\FloatBarrier
\section{Empirical Application}
\label{sec:empirical}
In this section, we revisit the large-scale randomized field experiment of \citet{ferraro2013using},\footnote{The data and original replication files are publicly available through Harvard Dataverse at \url{https://doi.org/10.7910/DVN1/22633}.} conducted in partnership with the Cobb County Water System in metropolitan Atlanta to examine how nonpecuniary incentives affect residential water consumption. More than 100,000 households were randomly assigned within meter-reading routes to a control group or one of three interventions. The first, \textit{Technical Advice}, provided practical water-saving tips. The second, \textit{Weak Social Norm}, supplemented these tips with an appeal to households' responsibility to conserve water. The third, \textit{Strong Social Norm}, further provided personalized information on each household's previous water consumption relative to the county median and its percentile rank among county residents.
The experiment has since provided evidence on both the effectiveness of these interventions and heterogeneity in household responses. \citet{ferraro2013using} find that technical advice alone has little impact, whereas prosocial appeals and personalized social comparisons produce progressively larger average water savings, with the latter particularly effective among households with higher baseline consumption. Building on these findings, \citet{ferraro2013heterogeneous} examine subgroup and quantile treatment effects, observing pronounced variation under the strong social-norm intervention but weaker evidence for the other messages. More recently, \citet{oka2026regression} show that the strong intervention shifts probability mass from higher to more moderate levels of water consumption. Despite these findings, a key question about treatment-effect heterogeneity remains formally untested: that is, whether each intervention changes the marginal distribution of water consumption only through an unknown common shift or also alters its shape.
To investigate this question, we apply the proposed eMBT tests to each intervention relative to the control group. Following \citet{ferraro2013using}, the outcome is total household water consumption during the summer of 2007, measured in thousands of gallons and constructed from July--October billing records.\footnote{We follow the sample exclusions in \citet{ferraro2013using} for discontinued water service, unusable billing records, and catastrophic leaks. The resulting sample reproduces the group sizes and mean outcomes reported in their Table~1. Due to billing lags, July--October records correspond to June--September consumption.} The resulting sample contains 71,643 control households, 11,675 households assigned to Technical Advice, 11,675 assigned to Weak Social Norm, and 11,676 assigned to Strong Social Norm. For each comparison, we implement the KS and CvM tests and their pointwise-studentized versions, both with and without covariate adjustment. For adjustment, we use the two pretreatment water-consumption measures included in Model A of \citet{ferraro2013using}, covering the previous watering season and spring 2007. All $p$-values are computed using 9,999 Mammen multiplier bootstrap draws.
\begin{table}[t]
\centering
\begin{threeparttable}
\caption{Tests for treatment-effect heterogeneity in household water consumption}
\label{tab:water-application}
\small
\setlength{\tabcolsep}{4pt}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lcccccccc@{}}
\toprule
& \multicolumn{4}{c}{Unadjusted} & \multicolumn{4}{c}{Covariate adjusted} \\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}
Intervention & KS & CvM & S--KS & S--CvM
& KS & CvM & S--KS & S--CvM \\
\midrule
Technical Advice
& 0.3442 & 0.4216 & 0.4882 & 0.6040
& 0.4147 & 0.5060 & 0.6174 & 0.7705 \\
Weak Social Norm
& 0.4853 & 0.2361 & 0.0630 & 0.1084
& 0.4986 & 0.1906 & \textbf{0.0131} & \textbf{0.0453} \\
Strong Social Norm
& \textbf{0.0002} & \textbf{0.0001} & \textbf{0.0001} & \textbf{0.0001}
& \textbf{0.0001} & \textbf{0.0001} & \textbf{0.0001} & \textbf{0.0001} \\
\bottomrule
\end{tabular*}
\begin{tablenotes}[flushleft]
\footnotesize
\item Notes: Entries are $p$-values based on 9,999 Mammen multiplier bootstrap draws. S denotes pointwise studentization; covariate adjustment uses the two pretreatment water-consumption measures described in the text.
Bold entries indicate $p\leq 0.05$.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table~\ref{tab:water-application} reports the results of testing for treatment-effect heterogeneity under each intervention. For Technical Advice, none of the tests provides significant evidence of heterogeneity. By contrast, all eight tests reject for Strong Social Norm, with $p$-values no greater than 0.0002, providing strong evidence of heterogeneous responses to personalized social comparisons. The evidence for Weak Social Norm is less uniform. The unadjusted KS test yields a $p$-value of 0.4853, which falls to 0.0630 with pointwise studentization and further to 0.0131 when covariate adjustment is also applied. The covariate-adjusted S--CvM test also rejects at the 5\% level, whereas neither unstudentized test rejects, even after adjustment. We view these findings as providing suggestive evidence of treatment-effect heterogeneity under Weak Social Norm, although less robust than that for Strong Social Norm. As shown in Section~\ref{sec:simu}, studentization and covariate adjustment often improve power by accounting for pointwise variability and exploiting baseline information, respectively, which may help explain the different rejection outcomes. Notably, in our analysis, covariate adjustment yields median pointwise variance reductions of approximately 43--47\% across the three comparisons. Similar precision gains from regression adjustment are reported by \citet{oka2026regression} in their analysis of the same experiment.
To better understand these differences, Figure~\ref{fig:water-application} plots the mean-anchored ETE curves for the three interventions, with and without covariate adjustment. For Strong Social Norm, the unadjusted curve declines sharply, from approximately 1.0 thousand gallons at $\tau=0.1$ to $-2.1$ at $\tau=0.9$. Since the curves measure deviations from the average treatment effect, this decline indicates smaller water-use reductions at lower expectiles and larger reductions at higher ones. Equivalently, after aligning the two distributions at their means, lower expectiles move upward and upper expectiles move downward, with a larger relative movement at the upper end. The resulting $0.1$--$0.9$ inter-expectile spread is approximately 3.1 thousand gallons narrower under treatment. Weak Social Norm exhibits a similar but less pronounced pattern, whereas Technical Advice shows no comparable contraction. These patterns remain largely unchanged after covariate adjustment. This evidence is consistent with \citet[Table~3]{ferraro2013using}, who report substantially larger water savings among households with higher baseline consumption (2.35 versus 0.57 thousand gallons), and complements the distributional changes documented by \citet{oka2026regression}. One economic interpretation, consistent with the mechanism evidence in \citet{ferraro2013heterogeneous}, is that personalized social comparisons make conservation norms more salient to high-use households, which may also have greater scope to reduce discretionary consumption, particularly outdoor water use. The contrast across interventions suggests that personalized social comparisons can both lower average water demand and narrow the gap between higher and lower consumption levels.
\begin{figure}[t]
\centering
\includegraphics[width=0.9\textwidth]{figures/figure_water_use_centered_ete_curves.pdf}
\caption{Mean-anchored expectile treatment-effect curves. Each curve plots $\widehat\Delta_E(\tau)-\widehat\Delta_E(1/2)$, measured in thousands of gallons, for an intervention relative
to the control group. The left and right panels show unadjusted and covariate-adjusted estimates, respectively.}
\label{fig:water-application}
\end{figure}
\FloatBarrier
\section{Conclusion}
\label{sec:con}
This paper develops tuning-free tests for treatment-effect heterogeneity via the sample expectile process. By centering the difference between the treatment and control expectile curves at the mean, we test whether treatment changes the marginal outcome distribution beyond a common shift. We establish the asymptotic properties of the resulting process, accounting for estimation of the unknown common effect, and construct a fast multiplier bootstrap that requires no density or score estimation, bandwidth selection, or repeated refitting. We also develop pointwise-studentized tests and covariate adjustment that reduces asymptotic variation without altering the population discrepancy. Simulations show generally reliable size control, competitive power, and substantial computational savings. In a large-scale water-conservation experiment, the tests provide strong evidence of heterogeneous responses to personalized social comparisons, with additional evidence for the weaker social-norm intervention emerging after studentization and covariate adjustment.
Several avenues for further research are worth exploring. Extending the proposed tests to observational studies would broaden their applicability beyond randomized experiments, while adapting the framework to panel data and difference-in-differences settings would allow researchers to study distributional policy effects over time. Another promising direction is to develop expectile-based tests for stochastic dominance and related partial orderings of marginal potential-outcome distributions, moving beyond the detection of heterogeneity to assessing the direction of distributional changes. We leave these interesting extensions for future research.
\begingroup \setstretch{1}
\setlength{\bibsep}{0pt}
\putbib
\endgroup
\end{bibunit}
\clearpage
\begin{bibunit}