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.
75,691 characters
Testing Conditional Stochastic Dominance via Copula Derivatives
\maketitle
\begin{abstract}
Comparing two populations at the same physical covariate value requires more than conditional means or isolated target-point decisions: researchers may need evidence about an entire conditional-distribution ordering over a continuum, even when covariate margins differ. This paper makes that common-value comparison estimable under an explicit structure--flexibility tradeoff and turns the resulting surface into simultaneous evidence for first-order stochastic dominance. Population-specific margins map the common covariate value into each group, while a fitted copula-derivative representation links conditional distributions across the region. Uniform inference propagates uncertainty from both the margins and dependence model through a one-sided statistic with unknown binding locations. Under correct specification within a finite copula class, smoothness and trimming conditions, and a uniquely best candidate family, the procedure admits uniform control and consistent calibration. Simulations show increasing rejection as alternatives become more distinguishable, alongside model-selection sensitivity and small-sample size distortion. In a descriptive PSID application, the high--low parental-education comparison satisfies the two-direction criterion after multiplicity adjustment, whereas adjacent education-group comparisons remain inconclusive. The framework therefore supports region-wide distributional comparison while making its structural and inferential boundaries explicit.
\end{abstract}
\section{Introduction}
Comparing populations at a common covariate value is often a distributional problem rather than a question about conditional means. Income opportunities, treatment responses, risks, and performance can differ in dispersion, tails, or crossing patterns even when their averages are similar. The relevant object is therefore the full conditional distribution, $H_g(y\mid x)=\mathbb{P}(Y_g\leq y\mid X_g=x)$. Researchers may also need to know whether an ordering holds throughout an empirically relevant conditioning region, not only at selected values. This distinction is especially consequential when populations have different covariate margins: a common percentile and a common physical value of $x$ then define different comparisons. This paper studies first-order dominance at common values of a continuous covariate over a trimmed continuum of conditioning and outcome values, while allowing each population to retain its own covariate margin.
The paper belongs to a broad inferential lineage rather than a narrow estimator class. Unconditional stochastic-dominance testing established how distributional inequalities can be evaluated without reducing welfare or risk comparisons to a few moments \citep{McFadden1989,Anderson1996,DavidsonDuclos2000,BarrettDonald2003,LintonMaasoumiWhang2005,LintonSongWhang2010}. Conditional work has since developed several complementary branches: global tests and conditional-moment formulations \citep{DelgadoEscanciano2013,AndrewsShi2013,AndrewsShi2017}, conditional treatment and distributional-effect tests \citep{Abadie2002,ChangLeeWhang2015}, inference at local or target covariate values \citep{DonaldHsuBarrett2012,ShenZhang2016,GoldmanKaplan2018,QuYoon2019,BugniCanayKim2025}, and structured or dynamic dominance procedures \citep{GonzaloOlmo2014,LintonSeoWhang2023}. These branches establish a substantial intellectual neighborhood for comparisons of conditional distributions, but they address different conditioning domains, statistical objects, and inferential targets.
The principal benchmarks sharpen the research question at hand. Delgado and Escanciano test global conditional-dominance restrictions through integrated joint-distribution differences \citep{DelgadoEscanciano2013}, while Chang, Lee, and Whang directly test conditional stochastic dominance between treatment and control distributions over covariate values \citep{ChangLeeWhang2015}. Linton, Seo, and Whang study first- and second-order dominance with high- or growing-dimensional dynamic information under a location–scale structure \citep{LintonSeoWhang2023}, and Gonzalo and Olmo develop conditional-dominance tests in dynamic settings \citep{GonzaloOlmo2014}. These papers establish that neither continuum-wide CSD nor structured conditional dominance is new. A distinct finite/local branch asks for distributional inference at selected conditioning values: Bugni, Canay, and Kim consider one or finitely many prespecified targets \citep{BugniCanayKim2025}, with related cutoff and local-distribution procedures developed by Shen and Zhang and Goldman and Kaplan \citep{ShenZhang2016,GoldmanKaplan2018}. The present paper addresses a different inferential target: one simultaneous two-population ordering over a trimmed continuum of common physical $(x,y)$ values under separate margins and an explicit dependence class. This difference matters because local decisions cannot authorize a region-wide conclusion, while the principal global benchmarks use different statistical representations and maintained structures.
Moving from selected targets to a surface changes the inferential problem. The inequality may bind at unknown locations on a two-dimensional continuum, while estimated marginal mappings and dependence parameters affect every point of that surface. The paper addresses this difficulty through a copula derivative. The established representation links a conditional distribution to the derivative of its joint copula \citep{JanssenSwanepoelVeraverbeke2017}: population-specific margins locate the same physical $x$ within each population, while the copula supplies coherence across conditioning values. Combined with semiparametric copula estimation \citep{GenestGhoudiRivest1995,Tsukahara2005}, this yields one fitted conditional-distribution surface for each population. The copula derivative is thus enabling literature, not the novelty protagonist. It is useful here because its separation of marginal location and dependence matches the added common-value surface problem and makes the maintained structure transparent.
The paper makes two research-level contributions. First, it makes the full common-value, two-population conditional distribution surface estimable under an explicit structure–flexibility trade-off. This capability permits coherent comparison over the continuum rather than a collection of unrelated local fits. Empirical margins and a fitted dependence representation are the enabling components, not the contribution headline. Second, the paper turns that fitted surface into simultaneous evidence about a uniform ordering. Building on contact-set methods for stochastic dominance \citep{LintonSongWhang2010} and directional inference for nonsmooth functionals \citep{FangSantos2019}, the analysis propagates uncertainty from estimated margins and dependence through a supremum whose least-favorable locations are unknown. This contribution is inseparable from the research question: without uniform calibration, the surface would remain a visualization or a collection of unadjusted pointwise contrasts rather than a defensible statement of dominance.
The theory validates the complete procedure under a correctly specified finite copula class, smoothness and trimming conditions, and a uniquely best candidate family. At a high level, it establishes uniform control of the estimated conditional-distribution contrast and consistent calibration of the one-sided dominance statistic. The simulations show why both the guarantee and its scope matter. Rejection probabilities generally rise as alternatives become easier to distinguish, while small-sample calibration remains sensitive to model selection and design. The numerical evidence therefore supports the procedure as a workable, structured option, but does not equate asymptotic validity with automatic robustness to misspecification or selection ambiguity.
An application to intergenerational income mobility illustrates the substantive payoff. The analysis compares adult family-income-rank distributions across parental-education groups at common childhood-income ranks, with childhood-income margins differing across groups. The high–low education comparison satisfies the paper's two-direction criterion after adjustment for multiplicity, whereas the adjacent high–middle and middle–low comparisons remain inconclusive. The result identifies bounded endpoint distributional separation without forcing a complete education gradient or a causal interpretation. The remainder of the paper develops the structured comparison, establishes its inferential properties, evaluates its finite-sample behavior, and reports the mobility application.
\section{Preliminaries}
\label{sec:preliminaries}
Conditional stochastic dominance is used to compare outcome distributions across populations after accounting for covariates. Existing conditional dominance procedures often focus on objects of the form $\mathbb{P}(Y\le y\mid X\le x)$ or rely on nonparametric conditional distribution estimators. These approaches are useful, but they do not directly address the pointwise conditional distribution $\mathbb{P}(Y\le y\mid X=x)$ when the conditioning variable is continuously distributed. In that case, the event $\{X=x\}$ has probability zero, and the conditional distribution is not obtained by dividing the joint distribution function by the marginal distribution function.
This paper studies pointwise conditional stochastic dominance through a copula derivative representation. If the joint distribution of $(X,Y)$ admits the copula representation $F(x,y)=C\{F_X(x),F_Y(y)\}$ and the relevant derivatives exist, then
\[
\mathbb{P}(Y\le y\mid X=x)=\frac{\partial_x F(x,y)}{f_X(x)}
=\partial_u C\{F_X(x),F_Y(y)\}.
\]
Thus, the comparison of pointwise conditional distributions can be transformed into a comparison of copula partial derivatives evaluated at the marginal probability indices. This transformation avoids direct estimation of conditional densities or bandwidth-dependent conditional distribution functions.
The proposed test combines this identification result with a semiparametric copula estimator. The marginal distributions are estimated nonparametrically by empirical distribution functions. The dependence structure is estimated by a finite collection of parametric copula families using canonical maximum pseudo-likelihood, and the candidate models are averaged by smooth AIC weights. Under a unique best-fitting copula condition, the model-averaged estimator has the same first-order behavior as the oracle estimator based on the selected copula family. The resulting test is a one-sided Kolmogorov--Smirnov statistic for the supremum of the difference between the estimated conditional distributions. Critical values are obtained by simulating the limiting Gaussian process using estimated influence functions.
\subsection{Pointwise Conditional Distribution}
\label{subsec:pointwise_conditional_distribution}
Let \(Z_{gi}=(X_{gi},Y_{gi})\), \(i=1,\ldots,n_g\), denote an i.i.d. sample from population \(g\in\{1,2\}\). The two samples are independent, and
\[
\frac{n_1}{n_1+n_2}\to \lambda\in(0,1).
\]
Let \(F_g(x,y)\) be the joint distribution function of \((X_g,Y_g)\), with marginal distribution functions \(F_{gX}\) and \(F_{gY}\). Throughout, \(X_g\) is continuously distributed with density \(f_{gX}\). Inference is conducted on a compact trimmed set \(\Lambda_\varepsilon\) such that, for some fixed \(\varepsilon\in(0,1/2)\),
\[
F_{gX}(x),\,F_{gY}(y)\in[\varepsilon,1-\varepsilon],
\qquad
(x,y)\in\Lambda_\varepsilon,\quad g=1,2.
\]
This trimming excludes boundary regions where copula derivatives may be unstable.
By Sklar's theorem, there exists a copula \(C_g\) such that
\[
F_g(x,y)
=
C_g\!\left(F_{gX}(x),F_{gY}(y)\right).
\]
The object of interest is the pointwise conditional distribution of \(Y_g\) given \(X_g=x\),
\[
H_g(y\mid x)
=
\mathbb{P}(Y_g\le y\mid X_g=x).
\]
When \(X_g\) is continuous, \(H_g(y\mid x)\) is not equal to \(F_g(x,y)/F_{gX}(x)\). Instead, if \(F_g(x,y)\) is differentiable with respect to \(x\) and \(f_{gX}(x)>0\), then
\[
H_g(y\mid x)
=
\frac{\partial F_g(x,y)/\partial x}{f_{gX}(x)}.
\]
Using the copula representation and the chain rule,
\[
\frac{\partial F_g(x,y)}{\partial x}
=
\partial_u C_g\!\left(F_{gX}(x),F_{gY}(y)\right)
f_{gX}(x),
\]
where \(\partial_u C_g\) denotes the partial derivative of \(C_g\) with respect to its first argument. Therefore,
\[
H_g(y\mid x)
=
\partial_u C_g\!\left(F_{gX}(x),F_{gY}(y)\right).
\]
This identity is the key observation underlying our approach: the pointwise conditional distribution can be represented by the first partial derivative of the copula evaluated at the marginal probability levels of \((x,y)\).
\subsection{Hypotheses and Test Statistics}
\label{subsec:hypotheses}
We say that population 1 dominates population 2 in the pointwise conditional stochastic dominance sense over \(\Lambda_\varepsilon\) if
\[
H_1(y\mid x)\le H_2(y\mid x),
\qquad
\forall (x,y)\in\Lambda_\varepsilon.
\]
Equivalently, define the conditional dominance contrast
\[
\Delta(x,y)
=
H_1(y\mid x)-H_2(y\mid x).
\]
The null hypothesis is
\[
\mathcal H_0:
\Delta(x,y)\le 0,
\qquad
\forall (x,y)\in\Lambda_\varepsilon,
\]
or equivalently,
\[
\mathcal H_0:
\sup_{(x,y)\in\Lambda_\varepsilon}\Delta(x,y)\le 0.
\]
The alternative hypothesis is
\[
\mathcal H_1:
\sup_{(x,y)\in\Lambda_\varepsilon}\Delta(x,y)>0.
\]
Using the copula derivative representation, the dominance contrast can be written as
\[
\Delta(x,y)
=
\partial_u C_1\!\left(F_{1X}(x),F_{1Y}(y)\right)
-
\partial_u C_2\!\left(F_{2X}(x),F_{2Y}(y)\right).
\]
Thus, the testing problem is transformed into a comparison of two copula derivative functions evaluated at population-specific marginal probability levels.
The null hypothesis is one-sided and can be written as
\[
\mathcal H_0:
\sup_{(x,y)\in\Lambda_\varepsilon}\Delta(x,y)\le 0.
\]
Therefore, a natural test is based on the largest positive deviation of an estimated dominance contrast from zero.
Let \(\widehat H_g(y\mid x)\) denote an estimator of \(H_g(y\mid x)\), whose construction is described in the next subsection. Define
\[
\widehat\Delta(x,y)
=
\widehat H_1(y\mid x)-\widehat H_2(y\mid x),
\qquad
(x,y)\in\Lambda_\varepsilon .
\]
Let
\[
s_N
=
\sqrt{\frac{n_1n_2}{n_1+n_2}} .
\]
We consider the one-sided Kolmogorov--Smirnov statistic
\[
T_N
=
s_N
\sup_{(x,y)\in\Lambda_\varepsilon}
\widehat\Delta(x,y).
\]
Large positive values of \(T_N\) indicate that the estimated conditional distribution of population 1 exceeds that of population 2 at some point in \(\Lambda_\varepsilon\), and hence provide evidence against pointwise conditional stochastic dominance. The asymptotic distribution of \(T_N\), the construction of the relevant contact set, and the critical value calculation are developed in the next section.
\subsection{Estimation Strategy}
\label{subsec:estimation_strategy}
We estimate \(H_g(y\mid x)\) using a semiparametric copula approach. Let
\[
\mathcal M
=
\{C_1,\ldots,C_Q\}
\]
be a finite set of candidate parametric copula families. For each population \(g\), the method proceeds in three steps.
First, estimate the marginal distributions nonparametrically by the scaled empirical distribution functions
\[
\widehat F_{gX}(x)
=
\frac{1}{n_g+1}
\sum_{i=1}^{n_g}
\mathbf 1(X_{gi}\le x),
\qquad
\widehat F_{gY}(y)
=
\frac{1}{n_g+1}
\sum_{i=1}^{n_g}
\mathbf 1(Y_{gi}\le y).
\]
The corresponding pseudo-observations are
\[
\widehat U_{gi}
=
\widehat F_{gX}(X_{gi}),
\qquad
\widehat V_{gi}
=
\widehat F_{gY}(Y_{gi}).
\]
Second, for each candidate copula family \(C_\ell(\cdot,\cdot;\theta_\ell)\), estimate the copula parameter by canonical maximum pseudo-likelihood:
\[
\widehat\theta_{g\ell}
=
\arg\max_{\theta_\ell\in\Theta_\ell}
\sum_{i=1}^{n_g}
\log c_\ell(\widehat U_{gi},\widehat V_{gi};\theta_\ell),
\]
where \(c_\ell\) is the density associated with \(C_\ell\).
Third, assign AIC weights to the candidate copulas. Let
\[
\operatorname{AIC}_{g\ell}
=
-2
\sum_{i=1}^{n_g}
\log c_\ell(\widehat U_{gi},\widehat V_{gi};\widehat\theta_{g\ell})
+
2k_\ell,
\]
where \(k_\ell\) is the dimension of \(\theta_\ell\). Define
\[
\widehat w_{g\ell}
=
\frac{
\exp\!\left[
-\frac{1}{2}
\left\{
\operatorname{AIC}_{g\ell}
-
\min_{1\le j\le Q}\operatorname{AIC}_{gj}
\right\}
\right]
}{
\sum_{m=1}^{Q}
\exp\!\left[
-\frac{1}{2}
\left\{
\operatorname{AIC}_{gm}
-
\min_{1\le j\le Q}\operatorname{AIC}_{gj}
\right\}
\right]
}.
\]
The model-averaged estimator of the pointwise conditional distribution is
\[
\widehat H_g(y\mid x)
=
\sum_{\ell=1}^{Q}
\widehat w_{g\ell}
\,
\partial_u C_\ell
\!\left(
\widehat F_{gX}(x),
\widehat F_{gY}(y);
\widehat\theta_{g\ell}
\right).
\]
Substituting this estimator into the contrast gives the feasible process
\[
\widehat\Delta(x,y)
=
\widehat H_1(y\mid x)
-
\widehat H_2(y\mid x),
\]
which is used in the statistic \(T_N\) defined above.
\section{Asymptotic Theory}
\label{sec:asymptotic_theory}
This section establishes the asymptotic properties of the semiparametric
copula derivative estimator and the resulting supremum test statistic. The main
steps are as follows. First, we derive a uniform asymptotic linear
representation for \(\widehat H_g(y\mid x)\). Second, we combine the two-sample
weak convergence result with the directional differentiability of the supremum
functional to obtain the limiting distribution of \(T_N\). Finally, we construct
multiplier critical values based on estimated influence functions.
\subsection{Regularity Conditions}
\label{subsec:regularity_conditions}
For \(g\in\{1,2\}\), write
\[
u_g(x)=F_{gX}(x),
\qquad
v_g(y)=F_{gY}(y).
\]
Let \(\ell_g^\ast\) denote the population copula model selected in the limit,
and let \(\theta_g\) be the corresponding parameter. For notational simplicity,
write
\[
C_g(u,v)
=
C_{\ell_g^\ast}(u,v;\theta_g).
\]
All copula derivatives below are evaluated at
\((u_g(x),v_g(y);\theta_g)\). Specifically, define
\[
C_{g,uu}(x,y)
=
\partial_{uu} C_{\ell_g^\ast}
\{u_g(x),v_g(y);\theta_g\},
\]
\[
C_{g,uv}(x,y)
=
\partial_{uv} C_{\ell_g^\ast}
\{u_g(x),v_g(y);\theta_g\},
\]
and
\[
C_{g,u\theta}(x,y)
=
\partial_{u\theta} C_{\ell_g^\ast}
\{u_g(x),v_g(y);\theta_g\}.
\]
\begin{assumption}[Sampling, trimming, and smoothness]
\label{ass:basic_regular}
The two samples are independent, and
\(n_1/(n_1+n_2)\to\lambda\in(0,1)\). The set
\(\Lambda_\varepsilon\) is compact and satisfies
\(F_{gX}(x),F_{gY}(y)\in[\varepsilon,1-\varepsilon]\) for all
\((x,y)\in\Lambda_\varepsilon\) and \(g=1,2\). The marginal distributions are
continuous, and the density \(f_{gX}\) is bounded away from zero on the relevant
support. For each \(g\), the derivatives
\(\partial_u C_{\ell_g^\ast}\), \(\partial_{uu} C_{\ell_g^\ast}\),
\(\partial_{uv} C_{\ell_g^\ast}\), and
\(\partial_{u\theta} C_{\ell_g^\ast}\) exist, are uniformly bounded, and are
uniformly continuous on the trimmed domain.
\end{assumption}
\begin{assumption}[Candidate copulas and first-stage estimation]
\label{ass:copula_first_stage}
The candidate set $\mathcal M=\{C_1,\ldots,C_Q\}$ is finite and fixed.
For each population $g$ and each candidate family $\ell$, let
\[
M_{g\ell}(\theta)
=
E\{\log c_\ell(U_g,V_g;\theta)\}.
\]
The corresponding normalized sample pseudo-likelihood criterion is
\[
\widehat M_{g\ell}(\theta)
=
\frac{1}{n_g}
\sum_{i=1}^{n_g}
\log c_\ell(\widehat U_{gi},\widehat V_{gi};\theta).
\]
The criterion $M_{g\ell}$ has a unique maximizer
$\theta_{g\ell}^{\ast}$, and $\widehat M_{g\ell}$ satisfies
\[
\sup_{\theta\in\Theta_\ell}
\left|
\widehat M_{g\ell}(\theta)-M_{g\ell}(\theta)
\right|
\overset{p}{\longrightarrow}0.
\]
The true copula belongs to the candidate set and is represented by a
unique family $\ell_g^\ast$. Moreover, for every
$\ell\ne\ell_g^\ast$,
\[
\kappa_{g\ell}
=
M_{g\ell_g^\ast}(\theta_g)
-
M_{g\ell}(\theta_{g\ell}^{\ast})
>0.
\]
For the limiting family $\ell_g^\ast$, $\theta_g$ is an interior point of
$\Theta_{\ell_g^\ast}$. On a neighborhood of $\theta_g$, the copula
log-density score and its derivatives with respect to the parameter and the two
copula arguments are measurable, continuously differentiable in their relevant
arguments, and dominated by square-integrable envelopes. The corresponding
local score and derivative classes are $P$-Donsker, the score map is
continuously differentiable under marginal-distribution perturbations, and the
negative expected parameter Hessian is nonsingular.
The canonical maximum pseudo-likelihood estimator in family
$\ell_g^\ast$ satisfies
\[
\sqrt{n_g}(\widehat\theta_g-\theta_g)
=
\frac{1}{\sqrt{n_g}}
\sum_{i=1}^{n_g}\psi_g(Z_{gi})
+o_p(1),
\]
where $E\{\psi_g(Z_g)\}=0$,
$E\|\psi_g(Z_g)\|^2<\infty$, and
$\operatorname{Var}\{\psi_g(Z_g)\}$ is positive definite.
\end{assumption}
\begin{remark}[Influence function of the CMPL estimator]
\label{rem:psi_g}
Let
\[
U_g=F_{gX}(X_g),
\qquad
V_g=F_{gY}(Y_g),
\]
and let $(U_g',V_g')$ be an independent copy of $(U_g,V_g)$.
For the limiting copula family $\ell_g^\ast$, define the
log copula density
\[
\ell_g(u,v;\theta)
=
\log c_{\ell_g^\ast}(u,v;\theta).
\]
The score evaluated at the true parameter is
\[
s_g(u,v)
=
\left.
\nabla_\theta
\ell_g(u,v;\theta)
\right|_{\theta=\theta_g},
\]
and the sensitivity matrix is
\[
J_g
=
-E\left[
\left.
\nabla_{\theta\theta^\top}^{2}
\ell_g(U_g,V_g;\theta)
\right|_{\theta=\theta_g}
\right].
\]
For $u_0,v_0\in(0,1)$, define the marginal-estimation correction
functions
\[
\begin{aligned}
r_{gX}(u_0)
&=
E\left[
\partial_u s_g(U_g',V_g')
\left\{
\mathbf 1(u_0\le U_g')-U_g'
\right\}
\right],\\
r_{gY}(v_0)
&=
E\left[
\partial_v s_g(U_g',V_g')
\left\{
\mathbf 1(v_0\le V_g')-V_g'
\right\}
\right].
\end{aligned}
\]
Then the influence function in
Assumption~\ref{ass:copula_first_stage} is
\[
\psi_g(Z_g)
=
J_g^{-1}
\left[
s_g(U_g,V_g)
+
r_{gX}(U_g)
+
r_{gY}(V_g)
\right].
\]
The first term is the usual copula score, while the last two terms account
for replacing the unknown marginal distributions by empirical distribution
functions. The representation requires the combined corrected score to be
well defined and square integrable and $J_g$ to be nonsingular.
\end{remark}
Under Assumption~\ref{ass:copula_first_stage}, the AIC weights concentrate on
the uniquely identified copula model. Hence the model-averaged estimator is
first-order equivalent to the oracle estimator based on
\(\ell_g^\ast\).
\begin{lemma}[Oracle equivalence of the model-averaged estimator]
\label{lem:oracle_equivalence}
Under Assumptions~\ref{ass:basic_regular} and
\ref{ass:copula_first_stage},
\[
\sup_{(x,y)\in\Lambda_\varepsilon}
\left|
\widehat H_g(y\mid x)
-
\partial_u C_{\ell_g^\ast}
\{\widehat F_{gX}(x),\widehat F_{gY}(y);\widehat\theta_g\}
\right|
=
o_p(n_g^{-1/2}),
\]
for \(g=1,2\).
\end{lemma}
\subsection{Uniform Linear Expansion}
\label{subsec:uniform_linear_expansion}
Define the influence function of \(\widehat H_g(y\mid x)\) at
\(z=(x,y)\) by
\[
\begin{aligned}
\phi_{g,z}(Z_g)
=
&\ C_{g,uu}(x,y)
\left\{\mathbf 1(X_g\le x)-F_{gX}(x)\right\}
\\
&+
C_{g,uv}(x,y)
\left\{\mathbf 1(Y_g\le y)-F_{gY}(y)\right\}
\\
&+
C_{g,u\theta}(x,y)^\top
\psi_g(Z_g).
\end{aligned}
\]
The first two terms capture the effect of estimating the marginal
distributions, while the last term captures the effect of estimating the copula
parameter.
\begin{theorem}[Uniform asymptotic linearity]
\label{thm:linear_h}
Under Assumptions~\ref{ass:basic_regular} and
\ref{ass:copula_first_stage}, for \(g=1,2\),
\[
\sup_{(x,y)\in\Lambda_\varepsilon}
\left|
\sqrt{n_g}
\{\widehat H_g(y\mid x)-H_g(y\mid x)\}
-
\frac{1}{\sqrt{n_g}}
\sum_{i=1}^{n_g}
\phi_{g,(x,y)}(Z_{gi})
\right|
=
o_p(1).
\]
Consequently,
\[
\sqrt{n_g}
\{\widehat H_g-H_g\}
\rightsquigarrow
\mathbb G_g
\quad
\text{in }
\ell^\infty(\Lambda_\varepsilon),
\]
where \(\mathbb G_g\) is a tight mean-zero Gaussian process with covariance
kernel
\[
\Sigma_g(z,z')
=
E\{\phi_{g,z}(Z_g)\phi_{g,z'}(Z_g)\}.
\]
\end{theorem}
The empirical-process conditions required for
Theorem~\ref{thm:linear_h} are not imposed as separate high-level assumptions.
They follow from the VC property of lower-orthant indicator classes, the
bounded smoothness of the copula derivative weights on the trimmed domain, and
the finite-dimensional asymptotic linear representation of
\(\widehat\theta_g\).
\subsection{Asymptotic Properties of the Statistic}
\label{subsec:two_sample_supremum_limit}
Before stating the result, define
\[
\mathbb G_\Delta(z)
=
\sqrt{1-\lambda}\,\mathbb G_1(z)
-
\sqrt{\lambda}\,\mathbb G_2(z),
\qquad
z\in\Lambda_\varepsilon,
\]
where \(\mathbb G_1\) and \(\mathbb G_2\) are independent copies of the
population-specific Gaussian limits in Theorem~\ref{thm:linear_h}. The limit process has covariance kernel
\[
\Sigma_\Delta(z,z')
=
(1-\lambda)\Sigma_1(z,z')
+
\lambda\Sigma_2(z,z').
\]
Under the continuity and Donsker conditions used in
Theorem~\ref{thm:linear_h}, the process
\(\mathbb G_\Delta\) admits a version with almost surely uniformly
continuous sample paths on the compact set
\(\Lambda_\varepsilon\). Also define
\[
\Gamma^\ast(\Delta)
=
\left\{
z\in\Lambda_\varepsilon:
\Delta(z)
=
\sup_{\tilde z\in\Lambda_\varepsilon}\Delta(\tilde z)
\right\}.
\]
Under the boundary null, \(\sup_{z\in\Lambda_\varepsilon}\Delta(z)=0\), this set
coincides with the contact set
\[
\Gamma(\Delta)
=
\{z\in\Lambda_\varepsilon:\Delta(z)=0\}.
\]
\begin{theorem}[Two-sample weak convergence]
\label{thm:two_sample_sup_limit}
Under Assumptions~\ref{ass:basic_regular} and
\ref{ass:copula_first_stage},
\[
s_N
\{\widehat\Delta-\Delta\}
\rightsquigarrow
\mathbb G_\Delta
\quad
\text{in }
\ell^\infty(\Lambda_\varepsilon).
\]
Moreover,
\[
T_N
-
s_N
\sup_{z\in\Lambda_\varepsilon}
\Delta(z)
\rightsquigarrow
\sup_{z\in\Gamma^\ast(\Delta)}
\mathbb G_\Delta(z).
\]
\end{theorem}
\begin{remark}[Local alternatives]
\label{rem:local_power}
Let \(\Delta_0\in C(\Lambda_\varepsilon)\) be a boundary-null contrast satisfying
\[
\Delta_0(z)\le0
\quad\text{for all }z\in\Lambda_\varepsilon,
\qquad
\sup_{z\in\Lambda_\varepsilon}\Delta_0(z)=0,
\]
and consider a sequence of local perturbations such that, uniformly over
\(\Lambda_\varepsilon\),
\[
\Delta_N(z)
=
\Delta_0(z)
+
s_N^{-1}h(z)
+
o(s_N^{-1}),
\]
where \(h\in C(\Lambda_\varepsilon)\). Suppose the conditions of
Theorem~\ref{thm:two_sample_sup_limit} hold uniformly along this sequence and
the centered contrast process has the same Gaussian weak limit
\(\mathbb G_\Delta\). Then
\[
T_N
\rightsquigarrow
\sup_{z\in\Gamma(\Delta_0)}
\left\{
h(z)+\mathbb G_\Delta(z)
\right\},
\]
where
\(
\Gamma(\Delta_0)
=
\{z\in\Lambda_\varepsilon:\Delta_0(z)=0\}
\)
is the contact set of the limiting null.
Let \(c_{1-\alpha}(\Delta_0)\) denote the \((1-\alpha)\)-quantile of
\(\sup_{z\in\Gamma(\Delta_0)}\mathbb G_\Delta(z)\). If its distribution is
continuous at this quantile, the local asymptotic rejection probability is
\[
\mathbb{P}\left[
\sup_{z\in\Gamma(\Delta_0)}
\left\{h(z)+\mathbb G_\Delta(z)\right\}
>c_{1-\alpha}(\Delta_0)
\right].
\]
Thus, first-order local power is determined by the restriction of \(h\) to the
contact set. Positive drift on that set can increase rejection probability,
whereas directions satisfying
\(\sup_{z\in\Gamma(\Delta_0)}h(z)\le0\) do not constitute an outward
first-order perturbation of the null at its binding points.
\end{remark}
Let \(c_{1-\alpha}\) denote the
\((1-\alpha)\)-quantile of
\( \sup_{z\in\Gamma(\Delta)}\mathbb G_\Delta(z)\)
when \(\sup_{z\in\Lambda_\varepsilon}\Delta(z)=0\), and set
\(c_{1-\alpha}=0\) when \(\sup_{z\in\Lambda_\varepsilon}\Delta(z)<0\).
When the null is on the boundary, we assume that
\(\sup_{z\in\Gamma(\Delta)}\mathbb G_\Delta(z)\) is continuous at
\(c_{1-\alpha}\).
The decision rule is to reject \(\mathcal{H}_0\) whenever
\[
T_N>c_{1-\alpha}.
\]
\begin{theorem}[Oracle size control and consistency]
\label{thm:oracle_test_validity}
Suppose Assumptions~\ref{ass:basic_regular} and
\ref{ass:copula_first_stage} hold. Then, under \(\mathcal{H}_0\),
\[
\lim_{n_1,n_2\to\infty}
\mathbb{P}(T_N>c_{1-\alpha})
\le \alpha .
\]
Moreover, under \(\mathcal{H}_1\),
\[
\mathbb{P}(T_N>c_{1-\alpha})\to1 .
\]
\end{theorem}
\begin{remark}[Feasible size control]
\label{rem:feasible_size}
Theorem~\ref{thm:oracle_test_validity} states size control for the oracle
critical value \(c_{1-\alpha}\). In practice, the critical value is estimated by
the multiplier or bootstrap procedures of Section~\ref{sec:critical_values},
whose consistency is established in
Theorems~\ref{thm:multiplier} and~\ref{thm:bootstrap_cv_consistency}.
Combining the two results, under \(\mathcal{H}_0\),
\[
\limsup_{n_1,n_2\to\infty}
\mathbb{P}\left\{T_N>\widehat c_{1-\alpha}^{(m)}\right\}
\le \alpha,
\]
with equality along the boundary null
\(\sup_{z\in\Lambda_\varepsilon}\Delta(z)=0\) when the limiting distribution of
\(\sup_{z\in\Gamma(\Delta)}\mathbb G_\Delta(z)\) is continuous at
\(c_{1-\alpha}\). The analogous statement holds for \(\widehat c_{1-\alpha}^{*}\).
\end{remark}
The preceding results concern a prespecified direction. In applications,
however, we test both directions so that a dominance conclusion requires not
only failure to reject the proposed ordering but also rejection of the reverse
ordering. For
\(j,k\in\{1,2\}\), \(j\ne k\), define
\[
\Delta_{12}(z)=\Delta(z),
\qquad
\Delta_{21}(z)=-\Delta(z),
\]
and, correspondingly,
\[
\widehat\Delta_{12}(z)=\widehat\Delta(z),
\qquad
\widehat\Delta_{21}(z)=-\widehat\Delta(z).
\]
The two directional null hypotheses are
\[
\mathcal H_{jk}:
\Delta_{jk}(z)\leq0
\quad\text{for all }z\in\Lambda_\varepsilon,
\]
with test statistics
\[
T_{jk,N}
=
s_N\sup_{z\in\Lambda_\varepsilon}
\widehat\Delta_{jk}(z).
\]
Let \(\widehat c_{jk,1-\alpha}\) denote the corresponding feasible multiplier
or bootstrap critical value. The two-direction rule declares population 1
to dominate population 2 when the event
\[
\mathcal D_{12,N}
=
\left\{
T_{12,N}\leq\widehat c_{12,1-\alpha},
\quad
T_{21,N}>\widehat c_{21,1-\alpha}
\right\},
\]
occurs, and defines \(\mathcal D_{21,N}\) analogously.
The following result follows directly from
\[
\mathcal D_{12,N}
\subseteq
\{T_{21,N}>\widehat c_{21,1-\alpha}\}
\]
and
\[
\mathbb{P}(\mathcal D_{12,N})
\geq
1
-
\mathbb{P}(T_{12,N}>\widehat c_{12,1-\alpha})
-
\mathbb{P}(T_{21,N}\leq\widehat c_{21,1-\alpha}),
\]
together with directional size control and fixed-alternative consistency. No
independence assumption between the two directional tests is required.
\begin{corollary}[Two-direction operating characteristics]
\label{cor:two_direction_rule}
Suppose that both directional tests satisfy feasible size control and
fixed-alternative consistency.
\begin{enumerate}
\item If \(\mathcal H_{21}\) holds, then
\[
\limsup_{N\to\infty}
\mathbb{P}(\mathcal D_{12,N})
\leq\alpha.
\]
\item If \(\mathcal H_{12}\) holds and \(\mathcal H_{21}\) is false under
a fixed alternative, then
\[
\liminf_{N\to\infty}
\mathbb{P}(\mathcal D_{12,N})
\geq1-\alpha.
\]
If, in addition,
\[
\sup_{z\in\Lambda_\varepsilon}\Delta_{12}(z)\leq-\eta
\]
for some \(\eta>0\), then
\[
\mathbb{P}(\mathcal D_{12,N})\longrightarrow1.
\]
\item If both \(\mathcal H_{12}\) and \(\mathcal H_{21}\) are false under
fixed alternatives, then both directional tests reject with probability
approaching one and
\[
\mathbb{P}(
\mathcal D_{12,N}\cup\mathcal D_{21,N}
)
\longrightarrow0.
\]
\end{enumerate}
The conclusions hold symmetrically after reversing the population labels.
If Holm-adjusted rejection decisions are used to define
\(\mathcal D_{jk,N}\) across a finite family of directional tests, the
probability of making any directional declaration for which the reverse null
is true is asymptotically bounded by the familywise level \(\alpha\). In
particular, under equality, the probability of declaring either direction is
asymptotically bounded by \(\alpha\).
\end{corollary}
The fixed-alternative qualification is important. Under
\(s_N^{-1}\)-local deviations from the boundary, the corresponding
classification probabilities are governed by
Remark~\ref{rem:local_power} and need not converge to either zero or one.
\begin{remark}[Discretized evaluation region]
\label{rem:grid}
Let \(\Lambda_N\subseteq\Lambda_\varepsilon\) be a
\(J_N\times J_N\) grid with mesh size
\(h_N=\sup_{z\in\Lambda_\varepsilon}d(z,\Lambda_N)\). Suppose that the relevant
maximizers are interior, \(\Delta\) is twice continuously differentiable with
bounded Hessian near them, and the contrast and resampling processes are
stochastically equicontinuous at scale \(h_N\). If \(s_Nh_N^2\to0\), then
\[
s_N\left\{
\sup_{z\in\Lambda_\varepsilon}\widehat\Delta(z)
-\max_{z\in\Lambda_N}\widehat\Delta(z)
\right\}=o_p(1),
\]
so grid evaluation does not affect the first-order limiting distribution or
resampling calibration. For a regular grid, \(h_N=O(J_N^{-1})\), so it suffices
that \(s_N/J_N^2\to0\), or equivalently \(J_N/N^{1/4}\to\infty\) when
\(s_N\asymp\sqrt N\). Under standard quantile regularity conditions, the same
mesh order holds in probability for the pooled empirical-quantile grid.
The simulations and application use \(J=25\) as a finite-sample choice. For
all reported sample sizes and pairwise comparisons, \(s_N/25^2\le0.036\).
Nevertheless, a grid held fixed as \(N\to\infty\) need not be asymptotically
equivalent to the continuum and may miss violations between grid points.
\end{remark}
\section{Choice of Critical Values}
\label{sec:critical_values}
\subsection{Multiplier Critical Values}
\label{subsec:multiplier_cv}
Throughout this section, the notation \(\rightsquigarrow_\xi\) (respectively
\(\rightsquigarrow_B\)) denotes weak convergence of a random element generated
by the multiplier variables \(\{\xi_{gi}\}\) (respectively the bootstrap weights
\(\{B_{gi}^{*}\}\)), conditional on the observed data, in probability.
The limiting distribution in Theorem~\ref{thm:two_sample_sup_limit} depends on
the unknown contact set and the unknown covariance structure of
\(\mathbb G_\Delta\). We approximate it by an analytic influence-function
(AIF) multiplier procedure that simulates the limiting Gaussian process using
estimated influence functions. We refer to this procedure as the AIF method
below.
Let \(\widehat\ell_g\) denote the AIC-selected copula family for
population \(g\), and write
\[
\widehat\theta_g
=
\widehat\theta_{g\widehat\ell_g}.
\]
For the selected family, define the score
\[
s_{g}(u,v;\theta)
=
\partial_{\theta}
\log c_{\widehat\ell_g}(u,v;\theta),
\]
and let
\[
\widehat s_{gi}
=
s_g(\widehat U_{gi},\widehat V_{gi};
\widehat\theta_g).
\]
The estimated sensitivity matrix is
\[
\widehat J_g
=
-
\frac{1}{n_g}
\sum_{j=1}^{n_g}
\partial_\theta
s_g(\widehat U_{gj},\widehat V_{gj};
\widehat\theta_g).
\]
To account for estimation of the two marginal distributions, define
the observation-specific correction terms
\[
\begin{aligned}
\widehat r_{gX,i}
&=
\frac{1}{n_g}
\sum_{j=1}^{n_g}
\partial_u
s_g(\widehat U_{gj},\widehat V_{gj};
\widehat\theta_g)
\left\{
\mathbf 1(\widehat U_{gi}\le\widehat U_{gj})
-
\widehat U_{gj}
\right\},\\
\widehat r_{gY,i}
&=
\frac{1}{n_g}
\sum_{j=1}^{n_g}
\partial_v
s_g(\widehat U_{gj},\widehat V_{gj};
\widehat\theta_g)
\left\{
\mathbf 1(\widehat V_{gi}\le\widehat V_{gj})
-
\widehat V_{gj}
\right\}.
\end{aligned}
\]
Let
\[
\widehat e_{gi}
=
\widehat s_{gi}
+
\widehat r_{gX,i}
+
\widehat r_{gY,i},
\qquad
\overline e_g
=
\frac{1}{n_g}
\sum_{j=1}^{n_g}\widehat e_{gj}.
\]
The plug-in influence-function estimate is
\[
\widehat\psi_{gi}
=
\widehat J_g^{-1}
\left(
\widehat e_{gi}-\overline e_g
\right).
\label{eq:theta_if_plugin}
\]
For the independence copula, which has no estimated dependence
parameter, the corresponding parameter influence term is absent.
The score derivatives and the sensitivity matrix are evaluated
numerically when closed-form expressions are unavailable.
Using \(\widehat\psi_{gi}\), define
\[
\begin{aligned}
\widehat\phi_{g,z}(Z_{gi})
=
&\ \widehat C_{g,uu}(z)
\left\{\mathbf 1(X_{gi}\le x)-\widehat F_{gX}(x)\right\}
\\
&+
\widehat C_{g,uv}(z)
\left\{\mathbf 1(Y_{gi}\le y)-\widehat F_{gY}(y)\right\}
\\
&+
\widehat C_{g,u\theta}(z)^\top
\widehat\psi_{gi},
\end{aligned}
\]
where the copula derivatives are evaluated at
\[
\left(
\widehat F_{gX}(x),
\widehat F_{gY}(y);
\widehat\theta_{g\widehat\ell_g}
\right).
\]
Let \(\{\xi_{gi}\}\) be independent multipliers, independent of the data, with
\(E(\xi_{gi})=0\), \(E(\xi_{gi}^2)=1\), and finite \(2+\delta\) moment for some
\(\delta>0\). Define
\[
\widehat{\mathbb G}_N^{(m)}(z)
=
\sqrt{1-\widehat\lambda}
\frac{1}{\sqrt{n_1}}
\sum_{i=1}^{n_1}
\xi_{1i}\widehat\phi_{1,z}(Z_{1i})
-
\sqrt{\widehat\lambda}
\frac{1}{\sqrt{n_2}}
\sum_{i=1}^{n_2}
\xi_{2i}\widehat\phi_{2,z}(Z_{2i}),
\]
where \(\widehat\lambda=n_1/(n_1+n_2)\).
In the simulations and empirical application, we use i.i.d. Mammen
two-point multipliers \citep{Mammen1993},
\[
\xi=
\begin{cases}
(1-\sqrt{5})/2,
& \text{with probability }(\sqrt{5}+1)/(2\sqrt{5}),\\
(1+\sqrt{5})/2,
& \text{with probability }(\sqrt{5}-1)/(2\sqrt{5}),
\end{cases}
\]
such that $E[\xi]=0$, $E[\xi^2]=1$, and $E[\xi^3]=1$. The multipliers are generated independently across observations,
populations, and resampling draws.
Let \(a_N\to\infty\) and \(a_N/s_N\to0\). Estimate the contact set by
\begin{equation}
\label{eq:contact_set}
\widehat\Gamma_N
=
\left\{
z\in\Lambda_\varepsilon:
s_N\widehat\Delta(z)>-a_N
\right\}.
\end{equation}
The multiplier statistic is
\[
T_N^{(m)}
=
\sup_{z\in\widehat\Gamma_N}
\widehat{\mathbb G}_N^{(m)}(z),
\]
with the convention that \(T_N^{(m)}=0\) if
\(\widehat\Gamma_N=\varnothing\). Let \(\widehat c_{1-\alpha}^{\,(m)}\) be the
conditional \((1-\alpha)\)-quantile of \(T_N^{(m)}\) given the data.
\begin{theorem}[Consistency of multiplier critical values]
\label{thm:multiplier}
Suppose Assumptions~\ref{ass:basic_regular} and
\ref{ass:copula_first_stage} hold. Let the multipliers be independent of the
data, have mean zero, variance one, and finite \(2+\delta\) moment for some
\(\delta>0\). Under the boundary null,
\[
T_N^{(m)}
\rightsquigarrow_\xi
\sup_{z\in\Gamma(\Delta)}
\mathbb G_\Delta(z)
\quad
\text{in probability}.
\]
If the distribution of
\(\sup_{z\in\Gamma(\Delta)}\mathbb G_\Delta(z)\) is continuous at its
\((1-\alpha)\)-quantile \(c_{1-\alpha}\), then
\[
\widehat c_{1-\alpha}^{\,(m)}
\overset{p}{\longrightarrow}
c_{1-\alpha}.
\]
\end{theorem}
\begin{remark}[Nonlinear AIF implementation]
\label{rem:nonlinear}
Theorem~\ref{thm:multiplier} establishes the first-order validity of the
influence-function multiplier procedure under the maintained
unique-best-family condition. In finite samples, however, several candidate
copula families may receive similar empirical support, so variation in the
model-selection step can remain non-negligible. To partially accommodate this
additional source of finite-sample uncertainty, our implementation augments the
first-order multiplier approximation with a nonlinear perturbation of the AIC
weights.
For each multiplier draw, we perturb both the candidate conditional CDFs by
their estimated influence functions and the candidate log-likelihoods by their
centered observation-level contributions, and then recompute the AIC weights
through the exact softmax map. In compact form,
\begin{align*}
\widehat H_g^{(b)}(z)
&=
\sum_{\ell=1}^{Q}\widehat w_{g\ell}^{(b)}
\left\{\widehat h_{g\ell}(z)
+\frac{1}{n_g}\sum_{i=1}^{n_g}
\xi_{gi}^{(b)}\widehat\phi_{g\ell,z}(Z_{gi})\right\},
\\
\widehat w_{g\ell}^{(b)}
&\propto
\exp\!\left\{-\frac{1}{2}\operatorname{AIC}_{g\ell}
+\sum_{i=1}^{n_g}\xi_{gi}^{(b)}\widehat q_{g\ell i}\right\},
\end{align*}
where \(\widehat q_{g\ell i}\) is the centered fitted log-density contribution.
Thus, the relative empirical support for the candidate families is allowed to
vary across multiplier draws, rather than being fixed at its full-sample value.
This nonlinear perturbation is a finite-sample stabilization device rather
than an additional first-order ingredient of the asymptotic theory. Under the
maintained unique-best-family condition, the AIC weight on the limiting family
converges to one, while the weights on the remaining families vanish
exponentially. The nonlinear weight perturbation is therefore asymptotically
inactive at first order and does not alter either the oracle theory or the
validity result in Theorem~\ref{thm:multiplier}.
\end{remark}
\subsection{Bootstrap Critical Values}
\label{subsec:bootstrap_cv}
The multiplier critical value in Section~\ref{subsec:multiplier_cv} is based on the
estimated influence-function representation and therefore avoids re-estimating the
copula parameters. As an alternative, we also consider a full weighted bootstrap
procedure that recomputes the marginal distributions, copula parameters, AIC
weights, and the conditional distribution estimator in each bootstrap replication.
To avoid confusion with the multiplier variables \(\xi_{gi}\) used above, we denote
the bootstrap frequency weights by \(B_{gi}^{*}\).
We impose the following condition on the bootstrap weights.
\begin{assumption}[Bootstrap weights]
\label{ass:bootstrap_weights}
For each population \(g\in\{1,2\}\), let
\(B_g^{*}=(B_{g1}^{*},\ldots,B_{gn_g}^{*})^\top\) be generated independently of the
data and independently across populations. The vector \(B_g^{*}\) is exchangeable,
\(B_{gi}^{*}\ge 0\), and \(\sum_{i=1}^{n_g}B_{gi}^{*}=n_g\). Moreover, for a generic
component \(B_{g1}^{*}\),
\[
\limsup_{n_g\to\infty}\|B_{g1}^{*}\|_{2,1}<\infty,
\qquad
\lim_{\lambda\to\infty}\limsup_{n_g\to\infty}
\sup_{t\ge \lambda}t^2\mathbb{P}(B_{g1}^{*}>t)=0,
\]
where \(\|B_{g1}^{*}\|_{2,1}=\int_0^\infty \mathbb{P}(B_{g1}^{*}\ge u)^{1/2}\,du\). Finally,
\[
\frac{1}{n_g}\sum_{i=1}^{n_g}(B_{gi}^{*}-1)^2
\overset{p}{\longrightarrow}1 .
\]
These conditions are satisfied, for example, by the multinomial bootstrap weights
\[
(B_{g1}^{*},\ldots,B_{gn_g}^{*})
\sim
\mathrm{Multinomial}\{n_g;(1/n_g,\ldots,1/n_g)\}.
\]
\end{assumption}
For each bootstrap replication, compute the weighted empirical marginal
distribution functions
\[
\widehat F_{gX}^{*}(x)
=
\frac{1}{n_g+1}\sum_{i=1}^{n_g}B_{gi}^{*}1(X_{gi}\le x),
\qquad
\widehat F_{gY}^{*}(y)
=
\frac{1}{n_g+1}\sum_{i=1}^{n_g}B_{gi}^{*}1(Y_{gi}\le y).
\]
The corresponding bootstrap pseudo-observations are
\[
\widehat U_{gi}^{*}=\widehat F_{gX}^{*}(X_{gi}),
\qquad
\widehat V_{gi}^{*}=\widehat F_{gY}^{*}(Y_{gi}).
\]
For each candidate copula family \(C_\ell(\cdot,\cdot;\theta_\ell)\), define the
normalized weighted pseudo-likelihood criterion
\[
\widehat M_{g\ell}^{*}(\theta)
=
\frac{1}{n_g}
\sum_{i=1}^{n_g}
B_{gi}^{*}
\log c_\ell(\widehat U_{gi}^{*},\widehat V_{gi}^{*};\theta)
\]
and its maximizer
\[
\widehat\theta_{g\ell}^{*}
=
\arg\max_{\theta_\ell\in\Theta_\ell}
\widehat M_{g\ell}^{*}(\theta_\ell).
\]
The bootstrap AIC criterion is
\[
AIC_{g\ell}^{*}
=
-2\sum_{i=1}^{n_g}
B_{gi}^{*}
\log c_\ell(\widehat U_{gi}^{*},\widehat V_{gi}^{*};\widehat\theta_{g\ell}^{*})
+2k_\ell,
\]
and the corresponding bootstrap AIC weights are
\[
\widehat w_{g\ell}^{*}
=
\frac{
\exp[-\{AIC_{g\ell}^{*}-\min_{1\le j\le Q}AIC_{gj}^{*}\}/2]
}{
\sum_{m=1}^Q
\exp[-\{AIC_{gm}^{*}-\min_{1\le j\le Q}AIC_{gj}^{*}\}/2]
} .
\]
The bootstrap estimator of the pointwise conditional distribution is
\[
\widehat H_g^{*}(y\mid x)
=
\sum_{\ell=1}^Q
\widehat w_{g\ell}^{*}
\partial_u C_\ell
\{\widehat F_{gX}^{*}(x),\widehat F_{gY}^{*}(y);
\widehat\theta_{g\ell}^{*}\}.
\]
\begin{assumption}[Bootstrap criterion regularity]
\label{ass:bootstrap_criterion}
For each population \(g\in\{1,2\}\), conditionally on the observed data and in
probability,
\[
\max_{1\le\ell\le Q}
\sup_{\theta\in\Theta_\ell}
\left|
\widehat M_{g\ell}^{*}(\theta)
-
\widehat M_{g\ell}(\theta)
\right|
=
o_p^{*}(1).
\]
\end{assumption}
\begin{remark}[Regularity of the bootstrap criterion]
Assumption~\ref{ass:bootstrap_criterion} is the conditional bootstrap analogue
of the uniform convergence requirement in
Assumption~\ref{ass:copula_first_stage}. It follows from the standard uniform
law for exchangeably weighted empirical processes and stochastic
equicontinuity of the smooth pseudo-likelihood criteria under empirical-margin
perturbations \citep{PraestgaardWellner1993Exchangeably,Tsukahara2005}. With a
fixed finite candidate set, compact parameter spaces away from singular
parameter values, and the usual integrable-envelope conditions for rank-based
canonical maximum pseudo-likelihood, this is a standard regularity condition
for the commonly used parametric copula families.
\end{remark}
Define the bootstrap contrast
\[
\widehat\Delta^{*}(x,y)
=
\widehat H_1^{*}(y\mid x)-\widehat H_2^{*}(y\mid x).
\]
Use the same estimated contact set
in (\ref{eq:contact_set}), the bootstrap statistic is
\[
T_N^{*}
=
\sup_{z\in\widehat\Gamma_N}
s_N\{\widehat\Delta^{*}(z)-\widehat\Delta(z)\},
\]
with the convention that \(T_N^{*}=0\) if \(\widehat\Gamma_N=\emptyset\). Let
\(\widehat c_{1-\alpha}^{\,*}\) be the conditional \((1-\alpha)\)-quantile of
\(T_N^{*}\) given the data.
\begin{theorem}[Consistency of bootstrap critical values]
\label{thm:bootstrap_cv_consistency}
Suppose Assumptions~\ref{ass:basic_regular},
\ref{ass:copula_first_stage}, \ref{ass:bootstrap_weights}, and
\ref{ass:bootstrap_criterion} hold. Under the boundary null
\(\sup_{z\in\Lambda_\varepsilon}\Delta(z)=0\),
\[
T_N^{*}
\rightsquigarrow_{B}
\sup_{z\in\Gamma(\Delta)}\mathbb{G}_\Delta(z)
\quad\text{in probability}.
\]
If the distribution of
\(\sup_{z\in\Gamma(\Delta)}\mathbb{G}_\Delta(z)\) is continuous at its
\((1-\alpha)\)-quantile \(c_{1-\alpha}\), then
\[
\widehat c_{1-\alpha}^{\,*}
\overset{p}{\longrightarrow}
c_{1-\alpha}.
\]
\end{theorem}
\subsection{Model-selection uncertainty and computational trade-offs}
\label{subsec:selection_tradeoff}
The multiplier and bootstrap procedures provide two different ways of
approximating the same first-order limiting distribution. Under
Assumption~\ref{ass:copula_first_stage}, the candidate set contains a uniquely
best copula family that is separated from the remaining candidates in the
population criterion. As a result, the AIC weight on the limiting family
converges to one, while the weights on the other families vanish
asymptotically. Model-selection uncertainty is therefore negligible at first
order under the maintained theory. This is why both
Theorem~\ref{thm:multiplier} and
Theorem~\ref{thm:bootstrap_cv_consistency} consistently estimate the same
oracle limiting distribution.
This asymptotic equivalence does not imply that model selection is innocuous
in finite samples. When several candidate copula families fit the data
similarly well, their AIC criteria can be close, so relatively small sampling
perturbations may induce substantial changes in the selected family or in the
associated model-averaging weights. Such variation is not a separate
first-order component under the unique-best-family asymptotics, but it can
nevertheless affect finite-sample calibration.
The two resampling procedures treat this source of uncertainty differently.
The full bootstrap recomputes the empirical margins, candidate-specific copula
parameters, AIC criteria, model weights, and conditional distribution
estimator within every bootstrap replication. It therefore reproduces the
entire estimation-and-selection pipeline and allows model support to vary
endogenously across resamples. By contrast, the influence-function multiplier
procedure starts from a first-order linear approximation around the fitted
sample. The nonlinear weight perturbation in
Remark~\ref{rem:nonlinear} is designed to make this approximation responsive
to finite-sample variation in relative model support, but it does not require
full re-estimation of all candidate copula models in every multiplier draw.
The distinction creates a practical trade-off. The bootstrap more directly
propagates uncertainty generated by the complete estimation and model-selection
procedure, but this comes at a substantially higher computational cost. The
influence-function multiplier procedure is much cheaper because the
candidate-model estimates and the required influence-function components are
computed once and subsequently reused across multiplier draws. Its purpose is
therefore not to replace the full bootstrap uniformly, but to provide a
computationally efficient approximation that retains the same first-order
validity under the maintained assumptions while partially accommodating
finite-sample model-selection variability. Section~\ref{sec:simulation}
examines how these differences translate into finite-sample size, power, and
computational performance.
\section{Simulation Study}
\label{sec:simulation}
This section evaluates the finite-sample calibration, power, and computational
trade-offs of the two inference procedures discussed in
Section~\ref{subsec:selection_tradeoff}. We consider three heterogeneous-\(X\)
designs. The two independent samples have equal size,
\(n_1=n_2=n\in\{50,100,200,500\}\). Every design--sample-size cell is based
on 1,000 Monte Carlo replications at the nominal 5\% level. Both procedures
use a \(25\times25\) empirical-quantile grid, \(\varepsilon=0.02\),
\(a_N=\sqrt{\log(n_1+n_2)}\), and 999
resampling draws.
We compare the AIF multiplier procedure with the
full bootstrap, using the same canonical maximum pseudo-likelihood and
AIC model-averaged point estimator in both cases. As described in
Section~\ref{sec:critical_values}, the full bootstrap recomputes the margins,
candidate copula estimates, AIC weights, and conditional CDFs in every
resampling draw, whereas the AIF procedure uses the nonlinear multiplier
implementation in Remark~\ref{rem:nonlinear} to perturb the fitted
candidate-specific conditional CDFs and their relative model weights without
fully re-estimating all candidate models in each draw.
The candidate set comprises the independence, Gaussian, Student-\(t\),
Clayton, Gumbel, Frank, and Joe copulas, together with the
survival, \(90^\circ\)-rotated, and \(270^\circ\)-rotated versions of the
Clayton and Gumbel copulas (13 families in total). For the Student-\(t\)
copula, the degrees of freedom range over 3, 4, 5, 7, 10, 15, and 30.
\subsection{Data-generating processes}
The tuning parameter \(\delta\) controls a location shift in population 2.
At \(\delta=0\), each design is a boundary-null experiment; positive values
generate \(H_1(y\mid x)-H_2(y\mid x)>0\) on part of the evaluation region.
\begin{itemize}
\item \textbf{P-DGP 1 (correctly specified Gaussian copula).}
Let \(\rho=0.5\), \(X^{(1)}\sim N(0,1)\), and
\(X^{(2)}\sim N(0.5,1)\). With independent standard normal errors,
\[
Y^{(1)}=\rho X^{(1)}+\sqrt{1-\rho^2}\,e^{(1)},\qquad
Y^{(2)}=\rho X^{(2)}+\delta+\sqrt{1-\rho^2}\,e^{(2)},
\]
where \(\delta\in\{0,0.1,\ldots,0.6\}\). The two \(X\) margins differ, but
the Gaussian copula is contained in the candidate set.
\item \textbf{P-DGP 2 (correctly specified Clayton copula with heterogeneous
margins).}
Independently for \(g=1,2\), draw \((U_g,V_g)\) from a Clayton copula with
parameter \(\theta=2\). Define
\[
q(p)=\bigl[|p-1/2|-0.4\bigr]_+^2,
\qquad
L(p)=\frac{\sqrt{3}}{\pi}\log\!\left(\frac{p}{1-p}\right),
\]
and set
\[
\begin{aligned}
X^{(1)}&=\Phi^{-1}(U_1),
&Y^{(1)}&=L(V_1)+10q(V_1),\\
X^{(2)}&=\Phi^{-1}(U_2)+10q(U_2),
&Y^{(2)}&=L(V_2)+\delta .
\end{aligned}
\]
The tail perturbation is zero on the central 80\% probability region. The
four group-specific marginal maps are all strictly increasing on \((0,1)\):
\(\Phi^{-1}\), \(L+10q\), \(\Phi^{-1}+10q\), and \(L+\delta\). Indeed,
\(10q'(p)\geq -2\), whereas
\((\Phi^{-1})'(p)\geq\sqrt{2\pi}>2\) and
\(L'(p)\geq 4\sqrt{3}/\pi>2\). Hence each \((X^{(g)},Y^{(g)})\) is
obtained from its latent Clayton pair through strictly increasing, though
group-specific, marginal transformations, so its copula remains exactly
Clayton with parameter \(\theta=2\). The perturbation creates markedly
different marginal distributions while preserving the central contact
region. We take \(\delta\in\{0,0.1,\ldots,0.6\}\).
\item \textbf{P-DGP 3 (misspecified nonlinear stress design).}
Let \(X^{(1)}\sim U(0,1)\), \(X^{(2)}\sim\operatorname{Beta}(2,2)\),
\(\sigma=0.5\), and
\[
m(x)=2(x-0.5)+0.5(x-0.5)^2.
\]
Generate
\[
Y^{(1)}=m(X^{(1)})+\sigma e^{(1)},\qquad
Y^{(2)}=m(X^{(2)})+\delta+\sigma e^{(2)},
\]
for \(\delta\in\{0,0.1,\ldots,0.5\}\). The conditional distribution is
identical across groups at \(\delta=0\), but the induced copula is generally
outside the finite candidate set. This design therefore provides robustness
evidence under misspecification and is not covered by the correct-
specification condition in Assumption~\ref{ass:copula_first_stage}.
\end{itemize}
\subsection{Simulation results}
Figures~\ref{fig:power-dgp1}--\ref{fig:power-dgp3} summarize the rejection
probabilities over \(\delta\). The value at \(\delta=0\) reports empirical
size, while positive values of \(\delta\) correspond to increasingly separated
alternatives. The AIF procedure is close to the nominal level in P-DGP 1 and
P-DGP 3 and is conservative in P-DGP 2. The bootstrap shows some
over-rejection at smaller sample sizes in P-DGP 1 and P-DGP 3 and is
conservative in P-DGP 2 at \(n=500\).
Across all three designs, power increases with both \(n\) and \(\delta\),
showing that both procedures become more discriminating as the alternative
moves farther from the boundary null and as sampling information increases.
The AIF and bootstrap procedures exhibit broadly similar power patterns
throughout the simulations. In P-DGP 1, the two power curves are closely
aligned under correct copula specification. A similarly close pattern is
observed in the misspecified P-DGP 3, which provides finite-sample robustness
evidence outside the formal validity conditions of the theory. P-DGP 2 shows
somewhat more conservative rejection behavior for the AIF procedure, but the
overall power pattern remains comparable across the two methods and the
difference narrows as \(n\) and \(\delta\) increase.
\begin{figure}[p]
\centering
\includegraphics[width=\textwidth]{figures/power_comparison_dgp1.pdf}
\caption{Power curves in P-DGP 1.}
\label{fig:power-dgp1}
\end{figure}
\begin{figure}[p]
\centering
\includegraphics[width=\textwidth]{figures/power_comparison_dgp2.pdf}
\caption{Power curves in P-DGP 2.}
\label{fig:power-dgp2}
\end{figure}
\begin{figure}[p]
\centering
\includegraphics[width=\textwidth]{figures/power_comparison_dgp3.pdf}
\caption{Power curves in the misspecified P-DGP 3.}
\label{fig:power-dgp3}
\end{figure}
Table~\ref{tab:dgp1-selection-size} examines size in P-DGP 1 and includes an
additional boundary-null experiment at \(n=1000\), using the same 1,000 Monte
Carlo replications and 999 resampling draws. The bootstrap method's
small-sample over-rejection diminishes as \(n\) increases, with the rejection
rate falling from 0.093 at \(n=50\) to 0.043 at \(n=500\) and 0.042 at
\(n=1000\). The AIF procedure remains closer to the nominal 5\% level at the
smaller sample sizes, with rejection rates between 0.045 and 0.059 across the
five sample sizes. Its rejection rate is not monotone in \(n\), rising from
0.045 at \(n=500\) to 0.057 at \(n=1000\), although this difference is modest
relative to the Monte Carlo standard error of approximately 0.007 for a
rejection probability near 0.05 based on 1,000 replications.
The selection strata help explain the finite-sample calibration pattern. The
proportion of replications in which both populations select the correctly
specified Gaussian family rises from 7.2\% at \(n=50\) to 69.8\% at \(n=500\)
and 83.4\% at \(n=1000\), indicating that model selection becomes progressively
more stable with sample size. For \(n\le500\), bootstrap rejection rates are
higher when at least one population selects a non-Gaussian family, while the
AIF procedure has lower rejection rates within this stratum. This ordering
disappears at \(n=1000\), when selection is substantially more stable. These
results are consistent with model-selection uncertainty being an important
source of finite-sample size distortion: when candidate copulas are not sharply
separated, sampling variation can alter the selected family or the associated
AIC weights, whereas this source of error diminishes as selection stabilizes.
The nonlinear AIF implementation is designed to mitigate this effect by
allowing relative model support to vary across multiplier draws and thereby
partially propagating model-selection uncertainty. This pattern is also
consistent with the maintained unique-best-family asymptotics, under which
selection uncertainty vanishes and the AIF and bootstrap procedures share the
same first-order limit.
\begin{table}[htbp]
\centering
\caption{P-DGP 1 empirical size by copula-selection event}
\label{tab:dgp1-selection-size}
\begin{threeparttable}
\footnotesize
\setlength{\tabcolsep}{3.2pt}
\begin{tabular}{@{}llccccc@{}}
\toprule
& &
\multicolumn{5}{c}{Sample size per group} \\
\cmidrule(lr){3-7}
Method
& Selection event
& \(n=50\)
& \(n=100\)
& \(n=200\)
& \(n=500\)
& \(n=1000\) \\
\midrule
\multirow{3}{*}{Bootstrap}
& Overall
&
\makecell[c]{0.093\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.090\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.066\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.043\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.042\\[-1pt]{\scriptsize(1000)}}
\\
& Both Gaussian
&
\makecell[c]{0.056\\[-1pt]{\scriptsize(72)}}
&
\makecell[c]{0.034\\[-1pt]{\scriptsize(178)}}
&
\makecell[c]{0.035\\[-1pt]{\scriptsize(423)}}
&
\makecell[c]{0.026\\[-1pt]{\scriptsize(698)}}
&
\makecell[c]{0.044\\[-1pt]{\scriptsize(834)}}
\\
& \makecell[l]{At least one\\non-Gaussian}
&
\makecell[c]{0.096\\[-1pt]{\scriptsize(928)}}
&
\makecell[c]{0.102\\[-1pt]{\scriptsize(822)}}
&
\makecell[c]{0.088\\[-1pt]{\scriptsize(577)}}
&
\makecell[c]{0.083\\[-1pt]{\scriptsize(302)}}
&
\makecell[c]{0.030\\[-1pt]{\scriptsize(166)}}
\\
\addlinespace[0.35em]
\multirow{3}{*}{AIF}
& Overall
&
\makecell[c]{0.059\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.059\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.058\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.045\\[-1pt]{\scriptsize(1000)}}
&
\makecell[c]{0.057\\[-1pt]{\scriptsize(1000)}}
\\
& Both Gaussian
&
\makecell[c]{0.042\\[-1pt]{\scriptsize(72)}}
&
\makecell[c]{0.039\\[-1pt]{\scriptsize(178)}}
&
\makecell[c]{0.043\\[-1pt]{\scriptsize(423)}}
&
\makecell[c]{0.033\\[-1pt]{\scriptsize(698)}}
&
\makecell[c]{0.061\\[-1pt]{\scriptsize(834)}}
\\
& \makecell[l]{At least one\\non-Gaussian}
&
\makecell[c]{0.060\\[-1pt]{\scriptsize(928)}}
&
\makecell[c]{0.063\\[-1pt]{\scriptsize(822)}}
&
\makecell[c]{0.069\\[-1pt]{\scriptsize(577)}}
&
\makecell[c]{0.073\\[-1pt]{\scriptsize(302)}}
&
\makecell[c]{0.036\\[-1pt]{\scriptsize(166)}}
\\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\scriptsize
\item \textit{Notes:}
Entries report rejection rates at the nominal 5\% level.
The number of Monte Carlo replications in each selection stratum is reported
in parentheses. ``Overall'' denotes the unconditional rejection rate across
all 1,000 replications.
\end{tablenotes}
\end{threeparttable}
\end{table}
Finally, the AIF method takes approximately 1--5 seconds per replication,
compared with 159--898 seconds for the bootstrap method. Figure~
\ref{fig:simulation-runtime} summarizes the difference across designs and
sample sizes.
\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{figures/runtime_comparison_dgp1_dgp3.pdf}
\caption{Mean computation time per Monte Carlo replication. The vertical
axis is logarithmic.}
\label{fig:simulation-runtime}
\end{figure}
\section{Empirical Application}
\label{sec:empirical}
We apply the test to intergenerational income mobility using the Panel Study of
Income Dynamics (PSID). The application asks whether the conditional
distribution of adult income rank differs across parental-education groups at
the same childhood parental-income rank. It is well suited to the method
because the marginal distribution of childhood income rank differs sharply
across education groups, while the estimand conditions on a common physical
rank value rather than requiring these margins to be equal.
\subsection{Sample and variables}
For each child, \(X\) is the average of the child's within-year parental-family-
income ranks observed at ages 13--17. The outcome \(Y\) is the average of the
child's within-year adult-family-income ranks observed at ages 30--36. We
require at least three valid annual observations in each age window and retain
the oldest eligible child from each origin family. Parental education is the
maximum completed schooling observed for the linked birth parents. We define
the low, middle, and high groups as at most 12 years, 13--15 years, and at
least 16 years, respectively. The main analysis is unweighted.
The final sample contains 1,047 individuals: 435 in the low group, 289 in the
middle group, and 323 in the high group. Table~\ref{tab:empirical-sample}
documents the pronounced differences in the \(X\) margins. Mean childhood
parental-income rank rises from 0.356 to 0.702 between the low and high
education groups; mean adult-income rank rises from 0.384 to 0.625. The
empirical CDFs in Figure~\ref{fig:empirical-x-margins} reinforce why a method
that permits population-specific conditioning-variable margins is useful.
\begin{table}[htbp]
\centering
\caption{PSID parental-education sample}
\label{tab:empirical-sample}
\footnotesize
\begin{tabular}{lrrrrr}
\toprule
Group & \(N\) & Mean \(X\) & Median \(X\) & Mean \(Y\) & Median \(Y\) \\
\midrule
Low (\(\leq12\) years) & 435 & 0.356 & 0.324 & 0.384 & 0.356 \\
Middle (13--15 years) & 289 & 0.495 & 0.490 & 0.486 & 0.486 \\
High (\(\geq16\) years) & 323 & 0.702 & 0.749 & 0.625 & 0.669 \\
\bottomrule
\end{tabular}
\end{table}
\begin{figure}[htbp]
\centering
\includegraphics[width=0.94\textwidth]{figures/parental_income_rank_ecdf_by_education.png}
\caption{Empirical distribution of childhood parental-income rank by parental-
education group.}
\label{fig:empirical-x-margins}
\end{figure}
As a conditional-mean benchmark, an OLS regression of \(Y\) on \(X\), group
indicators, and their interactions yields high- and middle-group indicator
\(p\)-values of 0.357 and 0.729 and interaction \(p\)-values of 0.190 and
0.325. The distributional comparison below therefore contains information
not summarized by this linear conditional-mean specification.
\subsection{Testing strategy and main results}
For every unordered pair, we test both directions. The predicted ordering is
that the higher-education group has a better adult-income distribution,
\(H_{\mathrm{higher}}(y\mid x)\leq H_{\mathrm{lower}}(y\mid x)\); the reverse
test interchanges the groups. We describe an ordering as supported only when the predicted null is not rejected, and the reverse null is rejected. This two-direction rule does not turn a failure to reject into proof of dominance,
but it rules out the opposite ordering. We adjust the six directional
\(p\)-values by the Holm procedure.
The baseline AIC model-average specification is based on pseudo-MLE,
\(\varepsilon=0.03\), a \(25\times25\) initial grid, and 1,999 multiplier
draws. For each unordered pair \((g,h)\), we set
\(a_N=\sqrt{\log(n_g+n_h)}\), as in the simulations. Let
\[
p_j=\varepsilon+\frac{j-1}{J-1}(1-2\varepsilon),
\qquad j=1,\ldots,J,\qquad J=25.
\]
For each pair, we construct a common grid as the Cartesian product of the
pooled-sample empirical \(p_j\)-quantiles of \(X\) and \(Y\), retaining only
points at which all four group-specific empirical margins lie in
\([\varepsilon,1-\varepsilon]\). Thus, the predicted and reverse tests for a given pair are evaluated at the
same \((x,y)\) values, although the evaluation grid may differ across pairs.
As a robustness check, we also use 1,999 replications of the full weighted
bootstrap. In each replication, the empirical margins are re-estimated, all
candidate copula families are refitted for each group, the AIC weights are
recomputed, and the conditional CDFs are reconstructed on the pair-specific
grid. The full bootstrap therefore propagates uncertainty from marginal
estimation, copula-parameter estimation, and model weighting through the
entire re-estimation procedure, rather than conditioning on the copula
families selected in the original sample.
Table~\ref{tab:empirical-tests} shows that only the high--low comparison
satisfies the two-direction criterion after multiplicity adjustment. The
predicted high-versus-low null is not rejected by either inference procedure
(raw \(p=1.0000\) for both). The reverse null is rejected by the
model-average multiplier procedure (raw \(p=0.0010\), Holm-adjusted
\(p=0.0060\)) and by the full weighted bootstrap (raw
\(p=0.0045\), Holm-adjusted \(p=0.0270\)). The high--middle and middle--low
comparisons remain inconclusive because neither direction is rejected after
Holm adjustment. The evidence therefore supports an endpoint separation
between high and low parental education, but not a complete stepwise
education gradient.
\begin{table}[htbp]
\centering
\caption{Directional conditional-dominance tests in the PSID application}
\label{tab:empirical-tests}
\footnotesize
\begin{tabular}{llcccc}
\toprule
Pair & Direction & \multicolumn{2}{c}{AIC-average multiplier} &
\multicolumn{2}{c}{Full AIC-average bootstrap} \\
\cmidrule(lr){3-4}\cmidrule(lr){5-6}
& & Raw \(p\) & Holm \(p\) & Raw \(p\) & Holm \(p\) \\
\midrule
High--Middle & Predicted & 0.9925 & 1.0000 & 0.9945 & 1.0000 \\
& Reverse & 0.2795 & 1.0000 & 0.3605 & 1.0000 \\
Middle--Low & Predicted & 0.9825 & 1.0000 & 0.9875 & 1.0000 \\
& Reverse & 0.0410 & 0.2050 & 0.0965 & 0.4825 \\
High--Low & Predicted & 1.0000 & 1.0000 & 1.0000 & 1.0000 \\
& Reverse & 0.0010 & 0.0060 & 0.0045 & 0.0270 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.94\textwidth}
\footnotesize\emph{Notes:} ``Predicted'' tests whether the higher parental-
education group dominates the lower group; ``Reverse'' tests the opposite
ordering. Holm adjustment is over all six directions within each inference
procedure. In every weighted-bootstrap draw, the empirical margins and all
13 candidate copulas are re-estimated and the AIC weights are recomputed.
Each inference procedure uses 1,999 draws.
\end{minipage}
\end{table}
The estimated AIC-averaged conditional-CDF difference for the high--low pair
ranges from \(-0.198\) to \(-0.033\) on its pair-specific common grid, has mean
\(-0.112\), and is nonpositive at every evaluated point. Thus, conditional on
the same observed parental-income rank, the high-education group's estimated
adult-income distribution is uniformly shifted to the right. The largest
gaps at parental-income ranks near 0.24, 0.50, and 0.76 are 0.137, 0.149, and
0.198 in CDF units. Figure~\ref{fig:empirical-heatmaps} displays all three
pairwise contrasts; blue cells are consistent with the predicted ordering.
The heatmaps show effect direction and magnitude, not pointwise
\(p\)-values---inference is based on the supremum statistic in
Table~\ref{tab:empirical-tests}.
\begin{figure}[p]
\centering
\includegraphics[width=\textwidth]{figures/education_gradient_difference_heatmaps.png}
\caption{AIC-averaged conditional-CDF differences: higher parental education
minus lower parental education. Negative values (blue) are consistent with
the predicted dominance ordering. The black curve marks estimated zero
crossings where present.}
\label{fig:empirical-heatmaps}
\end{figure}
\subsection{Robustness and interpretation}
The high--low conclusion remains unchanged when inference is based on the full
AIC-averaged weighted bootstrap. With 1,999 draws, the reverse null
has a raw \(p\)-value of 0.0045 and a Holm-adjusted \(p\)-value of 0.0270.
Multiplier-based sensitivity checks yield the same qualitative conclusion.
The Holm-adjusted \(p\)-value for the reverse test is 0.006 at trimming values
of 0.05 and 0.10; it is 0.009, 0.015, and 0.003 on 20-, 30-, and 40-point
grids, respectively, and 0.003 under the AIC single-best-family specification.
The full-sample AIC criterion selects Gaussian copulas for the low and middle
groups and a Frank copula for the high group. As a supplementary diagnostic
for these selected families, we apply White's information-matrix test to each
fitted copula pseudo-likelihood. The null hypothesis equates the negative
expected Hessian of the copula log-density with the expected outer product of
its score, as implied by correct specification of the selected copula. The
resulting \(p\)-values are 0.930, 0.120, and 0.970 for the low, middle, and high
groups, respectively, and the diagnostic therefore does not reject any of the
selected families. The test is applied to the full-sample AIC winner in each
group rather than to the AIC-averaged estimator. In addition, its
implementation treats the estimated margins as known and does not account for
the preceding model-selection step. We therefore view nonrejection only as
supplementary evidence against substantial misspecification, rather than as
evidence that the selected family is correctly specified or uniquely optimal.
The endpoint conclusion is more sensitive to changes in the target sample and
outcome definition. In the SRC-only subsample, the high--low reverse test has
a raw \(p\)-value of 0.013 and a Holm-adjusted \(p\)-value of 0.078. Replacing
adult family-income rank with individual labor-income rank yields corresponding
\(p\)-values of 0.024 and 0.144. The main finding should therefore be
interpreted as specific to the unweighted family-income analysis. It reflects
a descriptive conditional-distribution comparison rather than a causal effect
of parental education, which may also proxy for race, location, family
resources, and long-run selection. Population-representative inference would
further require survey weights to be incorporated consistently into the
empirical margins, copula estimation, and resampling procedure.
\section{Conclusion}
This paper makes a common-value, two-population conditional-distribution surface an operational object of inference when covariate margins differ across populations. The resulting comparison asks whether the entire outcome distributions are ordered along an empirically relevant continuum of physical covariate and outcome values. Its value is therefore research-level rather than estimator-specific: it permits a coherent region-wide comparison while preserving the distinction between a common physical value and a common percentile position. The accompanying inferential contribution turns the fitted surface into simultaneous evidence for a uniform ordering, rather than leaving it as a visualization or a collection of unadjusted pointwise contrasts.
The copula-derivative representation supplies the structure needed to implement this comparison. Population-specific margins locate the same physical covariate value within each group, and the fitted dependence model links conditional distributions across the region. Uniform process and contact-set arguments then propagate uncertainty from the estimated margins and dependence parameters through a one-sided statistic whose binding locations are unknown. The formal guarantees require correct specification within a finite copula class, smoothness and trimming conditions, and a uniquely best candidate family. Thus, the method offers an explicit structure--flexibility tradeoff rather than a claim that dependence modeling is costless.
The simulations reinforce both the value and the limits of the inferential
guarantee. Rejection probabilities generally increase as alternatives become
easier to distinguish, but small-sample calibration remains sensitive to
model selection and varies across designs. Evidence from the misspecified
design is informative about robustness without extending the theory beyond
its stated regime. In the descriptive PSID application, both inference
procedures support the high--low endpoint comparison. The predicted
high-versus-low null is not rejected, whereas the reverse null is rejected
after multiplicity adjustment by both the multiplier procedure and the full
weighted bootstrap; for the latter, the Holm-adjusted \(p\)-value
is 0.0270. The adjacent high--middle and middle--low comparisons remain
inconclusive. The empirical evidence therefore supports bounded endpoint
separation rather than a complete stepwise education gradient. This finding
is specific to the unweighted family-income analysis and does not establish a
causal effect of parental education. Future work can address model ties,
broader misspecification, and weighted population comparisons. More
generally, the paper shows how an explicit dependence structure can make
region-wide distributional comparisons feasible while keeping their
inferential and empirical boundaries visible.
\clearpage
\bibliographystyle{apalike}
\bibliography{sample}
\newpage