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.
135,998 characters
\raggedbottom
\widowpenalty=100
\clubpenalty=100
\thispagestyle{empty}
\begin{center}
\vspace{-0.3in}
{\large A Pairwise Differencing Distribution Regression Approach for Network Models\footnote{Acknowledgements: I would like to thank Frank Kleibergen,
Art\={u}ras Juodis and Bo Honor\'{e} for their advice and support.
I also thank Jad Beyhum, Otilia Boldea,
Pavel \v{C}\'\i\v{z}ek, Geert Dhaene,
Michal Koles\'{a}r, Louise Laage,
Elena Manresa, Konrad Menzel, Chris Muris, Cavit Pakel,
Mikkel Plagborg-M\o{}ller, Stephen Redding,
Sebastian Roelsgaard, Yassine Sbai Sassi, Timo Schenk,
Sami Stouli, Mark Watson,
Martin Weidner, Kaspar W\"{u}thrich, Andrei Zeleneev,
Lina Zhang and seminar participants at Oxford University,
University of Michigan, New York University, University
of Bristol, University of Exeter, University of
Melbourne, Erasmus University Rotterdam, University of
Groningen, Fordham University, University of Queensland,
Tilburg University, University of Manchester, Princeton
University, and participants at various workshops and
conferences for comments and discussions. Any errors are my own.}}
\vspace{0.1in}
{ Gabriela M. Miyazato Szini$^{\dag}$}\\
{\today}
\end{center}
\footnotetext{Tilburg University, e-mail: [email removed]}
\setcounter{footnote}{0}
\begin{center}\textbf{Abstract}\end{center}
\small\noindent
I develop an estimation and inference framework for
distribution regression in dyadic network settings with
two-way fixed effects that vary across thresholds of the outcome. I show that identification of the structural
parameters is achieved
through binarization of the outcome at each threshold,
and estimate the model by conditional maximum likelihood,
which ``differences out'' the fixed effects and circumvents
the incidental parameter problem. The estimator remains
asymptotically unbiased under sparsity, whether
from the network structure or binarization
at extreme thresholds. The second
novelty is to establish the joint asymptotic distribution of the
estimators across multiple thresholds with different convergence rates, and to develop simultaneous confidence bands and tests for equality of coefficients across thresholds. Monte Carlo simulations confirm small bias, valid inference, and correct simultaneous coverage under sparsity. An application to bilateral trade finds that coefficients vary substantially across the distribution, with equality rejected for key trade barriers.
\vspace{0.1in}
\noindent\textbf{Keywords:} Distribution regression, conditional maximum
likelihood, dyadic data, network models, sparsity, simultaneous confidence bands.
\noindent\textbf{JEL classification codes:} C14, C23, C24
\clearpage
\section{Introduction}
The vast majority of studies, especially for network models, estimate the effects of covariates on the mean of an outcome variable. However, in many settings, interest also lies in how the effects of covariates vary at different locations of the outcome distribution. For instance, in gravity models of international trade, the effects of trade barriers on bilateral exports may differ substantially across quantiles of trade flows. The distribution regression (DR) approach, initially proposed by \citet{foresi1995conditional} and developed further by \citet{chernozhukov2013inference}, addresses this by directly modelling the conditional distribution function through a sequence of binary choice models. At each threshold $y$, defined as a value of the outcome variable $y_{ij}$, the outcome is binarized into an indicator of whether it falls at or below $y$, and a binary response model is estimated. Varying the threshold characterizes how covariate effects differ across the distribution. Because the conditional distribution is approximated pointwise, the approach accommodates outcomes with point masses, such as variables bounded below at zero, without requiring smoothness of the conditional density.
A broad range of economic relationships involve bilateral
interactions between agents, naturally giving rise to network
data. Examples include models for international trade flows \citep{helpman2008estimating}, firm-level trade networks \citep{alfaro2023firm,bernard2022origins}, and earnings in employee-employer data \citep{bonhomme2019distributional}. In many of these settings, heterogeneity in the coefficients
of covariates across different locations of the outcome
distribution is of direct economic interest. Importantly, many of these networks are sparse, with only a small fraction of potential bilateral links realized.
Accordingly, I consider a directed network structure through a dyadic model (where the outcomes reflect pairwise interactions among the sampled units, \citeauthor{graham2020dyadic}, \citeyear{graham2020dyadic}) with additively separable two-way fixed effects that capture unobserved heterogeneity of senders and receivers. In many network settings, nodes (agents) differ substantially in their unobserved characteristics, and failing to account for this heterogeneity can bias the estimated effects of covariates. In addition, in network formation models, the fixed effects capture degree heterogeneity, that is, the tendency of some nodes to form more or stronger connections than others. In the distribution regression setting, the fixed effects play an analogous role: capturing node-specific shifts in the conditional
distribution at each threshold, allowing for systematic differences across senders and receivers in the level of the outcome. Both the structural parameters and the fixed effects are allowed to vary across levels of the outcome, providing
a flexible characterization of the conditional distribution. Moreover, the fixed effects are treated as unrestricted
parameters, with no distributional assumption imposed,
neither on their distribution across nodes nor on their
relationship with the covariates.
Each threshold of the distribution regression involves a nonlinear binary choice model with two-way fixed effects. Estimating these models by maximum likelihood gives rise to the incidental parameter problem \citep{neyman1948consistent}, since the number of fixed effects grows with the sample size, yielding asymptotically biased
estimates and invalid inference. This problem is compounded by network sparsity, which occurs through two channels in the distribution regression setting. First, in sparse networks with outcomes bounded below at, for instance,
zero, a large fraction of observations sit at the bound, so that the first threshold at which the binarized outcome varies
already corresponds to an extreme quantile, compressing
the estimable range into a sparse regime. At thresholds just above the bound, most binarized
outcomes equal one, leaving many units without
informative variation, rendering their fixed effects
unidentified. Second, even
in denser networks or with unbounded outcomes,
binarization at extreme thresholds generates sparsity,
since few observations fall in the tails.
To address these challenges, I develop a distribution regression framework based on the conditional maximum likelihood estimator (CMLE) of \citet{charbonneau2017multiple} and \citet{jochmans2018semiparametric}, originally proposed for network formation models. At each threshold, the approach relies on conditioning the likelihood function on a specific set of conditions for quadruples of nodes of the network, so that under a logistic specification the fixed effects are ``differenced out'' from the likelihood. Therefore, it yields estimates that are free of the incidental parameter problem and remain valid under sparsity.
The contributions of this paper are threefold. First, I
show that the structural parameters are identified in the
flexible specification where
both the coefficients and the fixed effects vary across
thresholds, and that identification remains valid under
network sparsity. The key mechanism is the binarization
of the outcome at each threshold, which serves not only
as an estimation tool but as an identification strategy.
The pointwise asymptotic properties at each threshold
follow from \citet{jochmans2018semiparametric}, and Monte
Carlo simulations confirm low bias and valid pointwise
inference across different degrees of network sparsity,
including at extreme thresholds.
Second, I develop methods for inference across thresholds. The object of interest in distribution regression is the entire profile of coefficients, and the empirical question is not only whether the effect of a covariate differs across the distribution, but where. Pointwise confidence intervals, plotted threshold by threshold, under-cover when interpreted jointly, while a joint equality test such as a Wald test delivers a single rejection decision without by itself indicating which parts of the distribution drive it. To address this, I derive the joint asymptotic distribution of the estimator across a finite number of thresholds, accommodating both the cross-threshold dependence (since the binary indicators at different thresholds are constructed from the same underlying outcome) and the varying convergence rates induced by the different degrees of sparsity at different locations of the distribution, without restricting these rates relative to one another. The main result is stated as a Gaussian approximation whose correlation matrix varies with the
sample size, which requires the dependence across thresholds to be non-degenerate but it does not impose it to have a fixed limit. Using this result, I construct simultaneous confidence bands via the sup-$t$ approach of \citet{montiel2019simultaneous}, which cover the entire coefficient path with prescribed probability and thereby indicate where in the distribution the effects differ; inverting the bands yields a test of equality of the coefficients across thresholds. The Wald test is also asymptotically valid and complements the sup-$t$ test, as the two have power against different alternatives. In finite samples, however, the sup-$t$ bands deliver correct simultaneous coverage and the equality test maintains size close to the nominal level across designs, while the Wald test shows size distortions for some covariates as the number of thresholds grows.
Third, I apply the method to gravity models of international trade, testing whether the effects of trade barriers vary
across the distribution of bilateral trade flows, a question with
direct implications for welfare analysis in the international trade
literature \citep{arkolakis2012new, melitz2015new,
bergstrand2025quantile}. The setting is a natural one for
distribution regression, since the outcome is bounded below at zero,
approximately 55\% of country pairs record no bilateral trade, and
the distribution among those with positive trade is heavily
right-skewed. Moreover, recent evidence suggests that the effects of
trade determinants vary substantially across the conditional
distribution of trade flows \citep{baltagi2016estimationmain,
bergstrand2025tails}, motivating methods that go beyond conditional
mean estimation. The estimated coefficients indeed vary substantially
across the distribution: the coefficient on distance increases from
approximately $1.05$ at the median to $1.78$ at the 99th percentile,
indicating that distance is a stronger barrier for the largest
bilateral trade relationships, and the sup-$t$ test rejects the null
of constant coefficients for distance, common legal system, and
border, providing formal evidence of heterogeneity in the structural
parameters across the distribution.
In general, conditional maximum likelihood methods eliminate the fixed effects by conditioning the likelihood on specific sets or sufficient statistics, avoiding their estimation entirely. In the network setting, \citet{charbonneau2017multiple} introduces the approach for directed networks, \citet{jochmans2018semiparametric} establishes its asymptotic theory under sparsity, and \citet{graham2017econometric} develops a related approach for undirected networks. Recent extensions include triadic network formation with dyad-level fixed effects \citep{muris2025triadic}, unified frameworks for static and dynamic binary choice in panels
and networks \citep{dano2025binary}, estimation for ordered outcomes in networks \citep{muris2025dyadic}, and two-way duration models \citep{roelsgaard2025essays}. Relatedly, \citet{bonhomme2024functional} derive moment
restrictions that eliminate fixed effects through functional
differencing.\footnote{Approaches that relax the logistic specification or
the additive structure of the fixed effects have also been
developed for network formation models.
\citet{candelaria2020semiparametric},
\citet{toth2017semiparametric}, and
\citet{gao2020nonparametric} study identification of the
structural parameters without a known parametric form for the
disturbance term, while \citet{zeleneev2026identification}
allows for nonparametric structure in the unobserved
heterogeneity. These methods are developed for single binary choice models,
primarily in undirected networks, and impose alternative
assumptions such as special regressors, continuously
distributed covariates, or compact support of the unobserved
heterogeneity.}
The present paper extends the framework from a single binary choice model to distribution regression, studying how structural parameters vary across the outcome distribution. The binarization at each threshold, which underpins
identification in this paper, is related to the strategy
used in the fixed-effects ordered logit literature to
identify a common coefficient under the proportional odds
assumption \citep{das1999panel, baetschmann2015consistent,
muris2017estimation}. More broadly, the identification argument of this paper extends to generalized ordered choice
models with fixed effects, including dyadic and network structures, where coefficients vary
across categories; to my knowledge no estimator is currently available for that setting.
The setting considered in this paper builds on
\citet{chernozhukov2020network}, who, to my knowledge, are the
first to propose distribution regression for network models with
two-way fixed effects. The key difference lies in the estimation
method employed at each threshold. They address the incidental
parameter problem at each threshold through analytical bias
corrections \citep{fernandez2016individual}. Several related approaches have
been proposed for single binary choice models in networks, including
analytical corrections for probit specifications
\citep{dzemski2019empirical} and jackknife bias
corrections \citep{hughes2026jackknife}. These methods
require dense network asymptotics, meaning that the link
probabilities must be bounded away from zero and one.\footnote{\citet{yan2019statistical} allow fixed
effects to grow at rate $\log n$, permitting some degree of
sparsity, but expected degrees must still grow nearly proportionally
to $n$.} In the distribution regression setting, this translates to requiring that the probability of the outcome falling below
a given threshold is bounded away from zero and one at each
threshold, a condition that is necessarily violated in the
tails of the distribution, where the binarized outcome has
little or no variation. In applications where, for instance, the outcome is bounded below at zero, and there is a high prevalence of zeros, the threshold at which the
binarized outcome first exhibits variation already corresponds
to an extreme quantile, so that this condition fails
broadly rather than only in the tails.
The conditional maximum likelihood approach does not require this assumption. The fixed effects are eliminated from the likelihood entirely, and the estimator remains consistent under sequences where the link probabilities can approach zero or one. The Monte Carlo simulations confirm this contrast in finite samples. An advantage of the bias correction approach, however, is that it delivers estimates of the fixed
effects, enabling the construction of counterfactual distributions
and average partial effects, objects that the conditional
likelihood approach cannot recover, since the fixed effects are
eliminated by conditioning. The two approaches are therefore complementary,
with the choice depending on whether the application requires
inference on structural parameters under sparsity or estimation of fixed effects for constructing
counterfactual distributions and average effects, which
requires a dense setting.
\vspace{0.25cm}
\noindent \textbf{Plan of the paper.} Section \ref{model_estimation} outlines the model and presents the estimation method; Section \ref{asymptotics} shows the pointwise asymptotic properties of the proposed estimator and discusses the sparsity conditions and the estimable quantile (threshold) range; Section \ref{joint_distribution} derives the joint distribution of the estimators across a finite number of thresholds, and constructs the sup-$t$ confidence bands and equality tests; Section \ref{simulations} presents the different settings for the Monte Carlo simulations and the obtained results; Section \ref{application} applies the method to gravity
models of international trade; and Section \ref{conclusion} concludes.
\section{Model and Estimation}\label{model_estimation}
\subsection{Distribution Regression Model for Networks}\label{subsection_model}
This section introduces the distribution regression model for a directed network structure formed through bilateral ties of units. To accommodate the network framework, a general dyadic setting is considered. More specifically, the conditional distribution function is parametrized as a function of dyad-specific characteristics and fixed effects for each unit in the observed pair of nodes.
The model follows \citet{chernozhukov2020network}, who
introduced the distribution regression framework with two-way
fixed effects for network data. Let $\{(y_{ij}, \boldsymbol{x}_{ij}) : (i,j) \in \mathcal{D} \}$ be the observed dataset, where $y_{ij}$ is a scalar outcome variable that can be discrete, continuous or mixed for a dyad $(i,j)$, and $\boldsymbol{x}_{ij}$ is a vector of dyad-specific covariates, such as measures of distance or similarity between units. Let $\mathcal{Y}$ be a region of interest in the support
of the outcome and $\mathcal{X} \subseteq \mathbb{R}^p$
the support of the covariates. The asymptotic theory is
developed for a finite collection of thresholds in
$\mathcal{Y}$; when
the outcome is discrete with finite support, $\mathcal{Y}$
can be taken as a subset of the support, while for
continuous outcomes, $\mathcal{Y}$ is a finite grid of
values in the support. The set of nodes\footnote{Thoughout the paper, I use units, nodes, or individuals interchangeably.} in the network is given by $\mathcal{N} = \{1,2, \dots, n\}$, and the set of observed
dyads is $\mathcal{D} = \{(i,j) : i \neq j,\; i,j \in
\mathcal{N}\}$, with $|\mathcal{D}| = n(n-1)$.\footnote{We consider that all the nodes are senders and receivers, but the method in this paper also allows for cases where the nodes that are senders differs from the nodes that are receivers, i.e., $i = 1, \dots, I$ and $j = 1, \dots J$, with $I \neq J$; and also for self-links to be formed.}
\enlargethispage{\baselineskip}
The unobserved heterogeneity of units $i$ and $j$ is captured by vectors $\boldsymbol{\nu}_i$ and $\boldsymbol{\omega}_j$ of unspecified dimension, whose relationship with the covariates
$\boldsymbol{x}_{ij}$ is left unrestricted. I assume that the conditional distribution of $y_{ij}$ given $(\boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)$ is given by:
\begin{align} \label{eq:model}
F_{y_{i j}}\left(y \mid \boldsymbol{x}_{i j}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j\right)=\Lambda \left(\boldsymbol{x}_{i j}^{\prime} \boldsymbol{\theta}_0(y)+\alpha\left(\boldsymbol{\nu}_{i}, y\right)+\gamma\left(\boldsymbol{\omega}_{j}, y\right)\right), \quad y \in \mathcal{Y}, \quad(i, j) \in \mathcal{D},
\end{align}
{\looseness=-1 \noindent where $\Lambda(\cdot)$ is a known link function, and $\boldsymbol{\theta}_0(y)$ is an unknown parameter vector of interest varying with $y$. $\alpha\left(\boldsymbol{\nu}_{i}, y\right)$ and $\gamma\left(\boldsymbol{\omega}_{j}, y\right)$ are unspecified measurable functions that can be seen as the unobserved individual fixed effects at a given level of $y$. This model is naturally semiparametric, not only because the parameters are allowed to vary with the output levels but also because it does not restrict how the individual unobserved effects correlate with the covariates. For notational simplicity, and without loss of generality, I denote $\boldsymbol{\theta}_0(y) = \boldsymbol{\theta}_{y,0}$, $\alpha_{i,y} = \alpha(\boldsymbol{\nu}_{i}, y)$ and $\gamma_{j,y} = \gamma(\boldsymbol{\omega}_{j}, y)$. When referring to a finite collection of thresholds
$\{y_1, \ldots, y_K\}$, the parameter vector at the
$k$-th threshold is denoted
$\boldsymbol{\theta}_{y_k}$, with $\theta_{y_k,d}$
its $d$-th element.\par}
Throughout this paper, $\Lambda(\cdot)$ is the logistic
distribution. While the distribution regression framework
accommodates general link functions
\citep{foresi1995conditional, chernozhukov2013inference},
the logistic specification is standard in fixed-effects
settings \citep{charbonneau2017multiple,
jochmans2018semiparametric, chernozhukov2020network}. This specification is essential for the conditional likelihood approach in Section \ref{subsection_estimation_method}, since it is the only link under which conditioning on sufficient statistics yields a likelihood free of nuisance parameters.
A key feature of this model is the two-way fixed effects,
which account for part of the dependence across
dyads. For instance, the outcomes for dyads $(i,j)$ and $(i,k)$ can be correlated through the shared sender effect
$\alpha_{i,y}$ and possible correlations in covariates
sharing index $i$. Allowing sender and receiver effects
to differ, together with $y_{ij}$ and $\boldsymbol{x}_{ij}$ not
necessarily equal to $y_{ji}$ and
$\boldsymbol{x}_{ji}$, accommodates directed
networks. However, notice that the model outlined in this section and the estimator proposed in the following section can be easily modified to accommodate undirected and bipartite networks. The additive separability of the fixed effects,
which is standard in the network formation literature
\citep{charbonneau2017multiple, jochmans2018semiparametric,
graham2017econometric}, is what enables the conditional
likelihood approach developed in the next section.
\enlargethispage{\baselineskip}
Finally, the conditional distribution $F_{y_{ij}} (y \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)$ can be written as:
\begin{align} \label{eq:binary}
F_{y_{ij}}\left(y \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j \right) &= \mathbb{E} [{1} \{y_{ij} \leq y \} \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j] \nonumber \\
&= \text{Pr}[\tilde{y}_{ij,y} = 1 \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j] \nonumber \\
&=\Lambda \left(\boldsymbol{x}_{i j}^{\prime} \boldsymbol{\theta}_{y,0}+\alpha_{i,y}+\gamma_{j,y}\right),
\end{align}
\noindent where $\tilde{y}_{ij,y} = 1\{y_{ij} \leq y \}$ is the binary indicator that the outcome falls at or below threshold $y$. The parameters $\boldsymbol{\theta}_{y,0}$ can therefore be estimated for each $y \in \mathcal{Y}$ as a sequence of binary (logistic) regressions with two-way fixed effects.\footnote{Distribution regression is preferred over
quantile regression (QR) in this setting because the
linear-in-parameters QR may poorly approximate the
conditional distribution when the outcome does not have a
smooth conditional density. In contrast, the DR
approximates the conditional distribution pointwise at
each threshold in the support of the outcome, making it
well-suited for outcomes with point masses without
requiring strong assumptions on how such masses are
generated \citep{chernozhukov2020network}. Moreover, to
my knowledge, no QR methods for two-way fixed effects in
dyadic network settings are currently available.} At each threshold, the maximum likelihood estimator is subject to the incidental parameter problem. Moreover, identification of individual fixed effects requires sufficient variation in the binary outcomes at each threshold, which may fail under sparsity. Both concerns motivate the conditional maximum likelihood approach developed in the next section, which eliminates the fixed effects from the likelihood by conditioning, and therefore does not require their identification or estimation.
\subsection{Conditional Maximum Likelihood
Estimation}\label{subsection_estimation_method}
The main challenge in the estimation of the sequence of binary regressions given by Equation \eqref{eq:binary} is that, even for a single binary regression, the incidental parameter problem \citep{neyman1948consistent} arises from estimating
models with two-way fixed effects by maximum likelihood. To circumvent this problem, I propose to estimate the parameters of the model $\boldsymbol{\theta}(y)$ for each threshold (for a given level $y$) independently, using the conditional maximum-likelihood method for network formation models suggested by \cite{charbonneau2017multiple} (for directed networks) and concurrently by \cite{graham2017econometric} (for undirected networks). This is applicable here because each binarized
threshold yields a binary choice model with two-way fixed
effects, identical in structure to a network formation
model. The approach extends the conditional maximum likelihood method for logistic panel data models with a single fixed effect \citep{rasch1960studies, chamberlain2010binary}\footnote{Also refer to \cite{arellano2001panel} for a survey.} to models with two-way fixed effects in dyadic structures, avoiding any
distributional assumption on the fixed effects.
The method relies on the existence of a set of conditions for quadruples of nodes in the observed network such that the fixed effects drop out of the conditional likelihood under a logistic link. The logistic specification is
essential: it is the only link under which the
conditional likelihood is entirely free of nuisance
parameters and takes a known closed
form.\footnote{In standard panel settings,
\citet{chamberlain2010binary} shows that for $T=2$,
parametric-rate estimation is only possible under
logistic errors. For $T \geq 3$,
\citet{davezies2023fixed} show that identification
beyond logistic is possible through conditional moment
restrictions estimated by GMM, though these do not yield
a closed-form conditional likelihood.}
From the conditional distribution $F_{y_{ij}}$ and the constructed binary variables $\tilde{y}_{ij,y}$, it follows that
$$ \tilde{y}_{ij,y} = 1\{\boldsymbol{x}_{ij}'\boldsymbol{\theta}_{y,0} + {\alpha}_{i,y} + {\gamma}_{j,y} + \varepsilon_{ij,y} \geq 0\}, \quad (i,j) \in \mathcal{D} $$
\noindent where $\varepsilon_{ij,y} \sim \text{i.i.d.}\
\text{Logistic}(0,1)$ across dyads for each threshold
$y$, with $\varepsilon_{ij,y} \perp
\{\boldsymbol{x}_{ij}, \alpha_{i,y}, \gamma_{j,y}\}$.
Under the logistic assumption, the conditional
probability takes the explicit form:
\begin{align} \label{eq:main_charbonneau}
\mathbb{E} [1\{y_{ij} \leq y\} \mid \boldsymbol{x}_{ij}, {\alpha}_{i,y}, {\gamma}_{j,y}] &= \text{Pr}[\tilde{y}_{ij,y} = 1 \mid \boldsymbol{x}_{ij}, {\alpha}_{i,y}, {\gamma}_{j,y}] \nonumber\\
&= \frac{\text{exp}(\boldsymbol{x}_{ij}'\boldsymbol{\theta}_{y,0} + {\alpha}_{i,y} + {\gamma}_{j,y})}{1+\text{exp}(\boldsymbol{x}_{ij}'\boldsymbol{\theta}_{y,0} + {\alpha}_{i,y} + {\gamma}_{j,y})}
\end{align}
The logistic specification ensures the existence of
sufficient statistics for the fixed effects, which is formalized in the Lemma below.
\begin{lemma}\label{lemma_sufficient}
Under the model specification given by Equation \eqref{eq:main_charbonneau}, the sums across each dimension of the pseudo panel, $\sum_{j=1}^n \tilde{y}_{ij,y}$ and $\sum_{i=1}^n \tilde{y}_{ij,y}$, are sufficient statistics for ${\alpha}_{i,y}$ and ${\gamma}_{j,y}$.
\end{lemma}
\noindent \textit{Proof.} See \ref{appendix_sufficient}.
This result is well-known for the standard panel case with one fixed effect, and \cite{graham2017econometric} establishes sufficiency of the degree sequence for undirected networks with a single set of node effects. Lemma \ref{lemma_sufficient} extends it to directed networks with two-way fixed effects, where both the row and column sums serve as sufficient statistics.
Even though one could construct a conditional maximum likelihood estimator based on the sufficient statistics, the maximization is intractable. \cite{charbonneau2017multiple} provides a tractable alternative by showing that it is possible to eliminate the two-way fixed effects by conditioning the above probability on the set of events $\{\tilde{y}_{ij,y} + \tilde{y}_{ik,y} = 1, \tilde{y}_{lj,y} + \tilde{y}_{lk,y} = 1, \tilde{y}_{ij,y} + \tilde{y}_{lk,y} = 1\}$ for different indices of senders and receivers $\{i,l;j,k\}$ (quadruples of nodes), such that:
\begin{align}
\text{Pr}[\tilde{y}_{ij,y} &= 1 \mid \boldsymbol{x}_{ij}, {\alpha}_{i,y}, {\gamma}_{j,y}, \tilde{y}_{ij,y} + \tilde{y}_{ik,y} = 1, \tilde{y}_{lj,y} + \tilde{y}_{lk,y} = 1, \tilde{y}_{ij,y} + \tilde{y}_{lk,y} = 1] \nonumber \\ &= \frac{\text{exp}(((\boldsymbol{x}_{ij} - \boldsymbol{x}_{ik}) - (\boldsymbol{x}_{lj} - \boldsymbol{x}_{lk})) '\boldsymbol{\theta}_{y,0})}{1+\text{exp}(((\boldsymbol{x}_{ij} - \boldsymbol{x}_{ik}) - (\boldsymbol{x}_{lj} - \boldsymbol{x}_{lk})) '\boldsymbol{\theta}_{y,0})},
\end{align}
\noindent which no longer depends on the fixed effects. This result follows from applying the classical conditional logit argument for static panel data models with a single fixed effect \citep{rasch1960studies} sequentially, first to eliminate the sender fixed effects, then the receiver fixed effects.
Summing over all quadruples that satisfy the conditioning
events, the CMLE maximizes:
\begin{align} \label{eq:main}
\sum_{i=1}^n \sum_{j=1, j \neq i}^n \sum_{l,k \in \mathbb{Z}_{ij,y}} \text{log} \left( \frac{\text{exp}(((\boldsymbol{x}_{ij} - \boldsymbol{x}_{ik}) - (\boldsymbol{x}_{lj} - \boldsymbol{x}_{lk})) '\boldsymbol{\theta}_y)}{1 + \text{exp}(((\boldsymbol{x}_{ij} - \boldsymbol{x}_{ik}) - (\boldsymbol{x}_{lj} - \boldsymbol{x}_{lk})) '\boldsymbol{\theta}_y)}\right) ,
\end{align}
where $\mathbb{Z}_{ij,y}$ is the set of nodes $k$ and
$l$ satisfying the conditioning events for pair $(i,j)$.
In practice, this reduces to a standard logit estimation
on pairwise-differenced outcomes and covariates, as
described in the next section.
\begin{figure}[h]
\centering
\captionsetup{skip=0pt}
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\begin{subfigure}[b]{0.40\textwidth}
\centering
\begin{tikzpicture}[scale=1.4]
\node[circle,draw,minimum size=0.5cm] (i) at (0,1.5) {\small i};
\node[circle,draw,minimum size=0.5cm] (j) at (1.5,1.5) {\small j};
\node[circle,draw,minimum size=0.5cm] (k) at (0,0) {\small k};
\node[circle,draw,minimum size=0.5cm] (l) at (1.5,0) {\small l};
\node[above=0.1cm of i, font=\tiny] {sender};
\node[above=0.1cm of j, font=\tiny] {receiver};
\node[below=0.1cm of k, font=\tiny] {receiver};
\node[below=0.1cm of l, font=\tiny] {sender};
\draw[->,thick,blue,line width=1.2pt] (i) -- (j);
\draw[->,thick,blue,line width=1.2pt] (l) -- (k);
\draw[->,dashed,red,line width=1.2pt] (i) -- (k);
\draw[->,dashed,red,line width=1.2pt] (l) -- (j);
\end{tikzpicture}
\caption{$z_\sigma = 1$}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.40\textwidth}
\centering
\begin{tikzpicture}[scale=1.4]
\node[circle,draw,minimum size=0.5cm] (i) at (0,1.5) {\small i};
\node[circle,draw,minimum size=0.5cm] (j) at (1.5,1.5) {\small j};
\node[circle,draw,minimum size=0.5cm] (k) at (0,0) {\small k};
\node[circle,draw,minimum size=0.5cm] (l) at (1.5,0) {\small l};
\node[above=0.1cm of i, font=\tiny] {sender};
\node[above=0.1cm of j, font=\tiny] {receiver};
\node[below=0.1cm of k, font=\tiny] {receiver};
\node[below=0.1cm of l, font=\tiny] {sender};
\draw[->,dashed,red,line width=1.2pt] (i) -- (j);
\draw[->,dashed,red,line width=1.2pt] (l) -- (k);
\draw[->,thick,blue,line width=1.2pt] (i) -- (k);
\draw[->,thick,blue,line width=1.2pt] (l) -- (j);
\end{tikzpicture}
\caption{$z_\sigma = -1$}
\end{subfigure}
\caption{Informative quadruple configurations with
\{i,l\} as senders and \{j,k\} as receivers. Blue solid
arrows: links present ($\tilde{y}=1$). Red dashed arrows:
links absent ($\tilde{y}=0$). The two senders connect to
opposite receivers, yielding $z_\sigma \in \{-1, 1\}$.
Values: (a) $\tilde{y}_{ij,y}=1, \tilde{y}_{ik,y}=0,
\tilde{y}_{lj,y}=0, \tilde{y}_{lk,y}=1$; (b)
$\tilde{y}_{ij,y}=0, \tilde{y}_{ik,y}=1,
\tilde{y}_{lj,y}=1, \tilde{y}_{lk,y}=0$.}
\label{figure1:correct}
\end{figure}
\begin{figure}[h]
\centering
\captionsetup{skip=0pt}
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\begin{subfigure}[b]{0.30\textwidth}
\centering
\begin{tikzpicture}[scale=1.4]
\node[circle,draw,minimum size=0.5cm] (i) at (0,1.5) {\small i};
\node[circle,draw,minimum size=0.5cm] (j) at (1.5,1.5) {\small j};
\node[circle,draw,minimum size=0.5cm] (k) at (0,0) {\small k};
\node[circle,draw,minimum size=0.5cm] (l) at (1.5,0) {\small l};
\node[above=0.1cm of i, font=\tiny] {sender};
\node[above=0.1cm of j, font=\tiny] {receiver};
\node[below=0.1cm of k, font=\tiny] {receiver};
\node[below=0.1cm of l, font=\tiny] {sender};
\draw[->,dashed,red,line width=1.2pt] (i) -- (j);
\draw[->,dashed,red,line width=1.2pt] (l) -- (k);
\draw[->,dashed,red,line width=1.2pt] (i) -- (k);
\draw[->,dashed,red,line width=1.2pt] (l) -- (j);
\end{tikzpicture}
\caption{All absent ($z_\sigma = 0$)}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.30\textwidth}
\centering
\begin{tikzpicture}[scale=1.4]
\node[circle,draw,minimum size=0.5cm] (i) at (0,1.5) {\small i};
\node[circle,draw,minimum size=0.5cm] (j) at (1.5,1.5) {\small j};
\node[circle,draw,minimum size=0.5cm] (k) at (0,0) {\small k};
\node[circle,draw,minimum size=0.5cm] (l) at (1.5,0) {\small l};
\node[above=0.1cm of i, font=\tiny] {sender};
\node[above=0.1cm of j, font=\tiny] {receiver};
\node[below=0.1cm of k, font=\tiny] {receiver};
\node[below=0.1cm of l, font=\tiny] {sender};
\draw[->,thick,red,line width=1.2pt] (i) -- (j);
\draw[->,thick,red,line width=1.2pt] (l) -- (k);
\draw[->,thick,red,line width=1.2pt] (i) -- (k);
\draw[->,thick,red,line width=1.2pt] (l) -- (j);
\end{tikzpicture}
\caption{All present ($z_\sigma = 0$)}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.30\textwidth}
\centering
\begin{tikzpicture}[scale=1.4]
\node[circle,draw,minimum size=0.5cm] (i) at (0,1.5) {\small i};
\node[circle,draw,minimum size=0.5cm] (j) at (1.5,1.5) {\small j};
\node[circle,draw,minimum size=0.5cm] (k) at (0,0) {\small k};
\node[circle,draw,minimum size=0.5cm] (l) at (1.5,0) {\small l};
\node[above=0.1cm of i, font=\tiny] {sender};
\node[above=0.1cm of j, font=\tiny] {receiver};
\node[below=0.1cm of k, font=\tiny] {receiver};
\node[below=0.1cm of l, font=\tiny] {sender};
\draw[->,thick,blue,line width=1.2pt] (i) -- (j);
\draw[->,thick,blue,line width=1.2pt] (l) -- (j);
\draw[->,dashed,red,line width=1.2pt] (i) -- (k);
\draw[->,dashed,red,line width=1.2pt] (l) -- (k);
\end{tikzpicture}
\caption{Same direction ($z_\sigma = 0$)}
\end{subfigure}
\caption{Examples of non-informative quadruple configurations with
\{i,l\} as senders and \{j,k\} as receivers. In (a), all
links are absent; in (b), all links are present; in (c),
both senders connect to the same receiver. In all three
cases $z_\sigma = 0$, so the quadruple provides no
information for estimation.}
\label{Figure3:exampletetrads}
\end{figure}
As in the static panel logit considered in \cite{rasch1960studies} and \cite{chamberlain2010binary}, only
quadruples whose outcomes exhibit variation (the
analogues of movers in the panel data literature)
contribute to the likelihood: a quadruple $\{i,l;j,k\}$ is informative only if each
node's outcomes vary across its two potential links. To illustrate this argument, Figure~\ref{figure1:correct} shows the two informative
configurations for a fixed sender-receiver assignment
$\{i,l;j,k\}$ ($i$ and $l$ are fixed to be senders; and $j$ and $k$ are fixed to be receivers). Note that there is variation in the outcomes for each node, i.e., each node has exactly one link
present and one absent, yielding $z_\sigma \in \{-1, 1\}$. In a directed network, the sender-receiver roles can
also be reversed; Figure~\ref{figure2:correct} in
\ref{identification} shows all twelve informative
configurations across all possible assignments. Figure~\ref{Figure3:exampletetrads} illustrates
three non-informative configurations, maintaining the same
sender-receiver assignment. In Subfigure (a), all links are
absent; in Subfigure (b), all links are present; in both
cases, no node's outcomes vary. In Subfigure (c), both
senders connect to the same receiver, so that while each
sender's outcomes vary, the receivers' do not, and the
pairwise difference yields $z_\sigma = 0$. None of these
quadruples contribute to the likelihood. For further
intuition on the identification of the common parameters, I
refer to \ref{identification}.
\begin{remark}[Conditioning across thresholds]
If the fixed effects were constant across different levels of the outcome $y$ (thresholds),
one could also condition on events combining information across thresholds, such as $\tilde{y}_{ij,y_1} + \tilde{y}_{ij,y_2} = 1$. However, this approach is not pursued in this paper. Allowing the fixed effects to vary flexibly across thresholds, as done in the current framework, provides a more general specification that avoids possible misspecification.
\end{remark}
A property of the model central to the asymptotic theory of this estimator is that the binary indicators are conditionally independent across dyads:
$$
\Pr(\tilde{y}_{ij,y}, \tilde{y}_{kl,y} \mid
\{\boldsymbol{x}_{ij}\}_{n,n}, \{\alpha_{i,y}, \gamma_{j,y}\}_n) =
\Pr(\tilde{y}_{ij,y} \mid \boldsymbol{x}_{ij}, \alpha_{i,y},
\gamma_{j,y}) \cdot \Pr(\tilde{y}_{kl,y} \mid \boldsymbol{x}_{kl},
\alpha_{k,y}, \gamma_{l,y})
$$
for all $(i,j) \neq (k,l)$, where $\{\boldsymbol{x}_{ij}\}_{n,n}$ denotes the entire set of covariates, and $\{\alpha_{i,y}, \gamma_{j,y}\}_{n}$ denotes the full set of fixed effects. The conditional independence follows from the errors $\varepsilon_{ij,y}$ being independent across dyads for each threshold $y$, even when
dyads share a node (with $\varepsilon_{ij,y}$ and
$\varepsilon_{ji,y}$ treated as distinct and independent).\footnote{One drawback is that, in a network formation model context, transitivity across the probabilities is not taken into account by this model. It rules out interdependent link preferences, where individuals' preferences over a link may vary with the presence or absence of links elsewhere in the network. However, \cite{dzemski2019empirical} shows that this dyadic structure can still reproduce the transitivity patterns observed in some datasets.} The model thus belongs to the class of conditionally independent dyad (CID) models \citep{fafchamps2007formation,
graham2020dyadic}.
While the conditional independence holds across dyads within each threshold, the model allows for correlation of $\varepsilon_{ij,y}$ and $\varepsilon_{ij,y'}$ across thresholds $y \neq y'$ for the same dyad, since estimation is performed
separately at each threshold. This
cross-threshold dependence becomes relevant when
deriving the joint distribution in Section
\ref{joint_distribution}.
\section{Pointwise Asymptotics}\label{asymptotics}
The pointwise asymptotic properties of the estimator $\boldsymbol{\theta}_{n,y}$ for a single threshold value $y$ of the conditional distribution follow from results in \cite{jochmans2018semiparametric} for the estimator of \cite{charbonneau2017multiple}. The proofs in \ref{appendix_asymptotics} adapt the same broad structure, with modifications that facilitate the extension to the joint asymptotic distribution across thresholds in Section \ref{joint_distribution}. Subsection \ref{subsection_sparsity} discusses their implications for the estimable range of thresholds in the distribution regression setting.
Throughout this section, the sequence of individual
effects $\{\alpha_{i,y}, \gamma_{j,y}\}_n$ are treated
as fixed, since the analysis conditions on them. The
asymptotic framework considers $n \to \infty$, so that
both dimensions of the network grow at the same rate. For an ordered quadruple of distinct nodes $\{i,l;j,k\}$ from $\mathcal{N}$, define:
$$ z_y(\sigma\{i,l;j,k\}) = \frac{(\tilde{y}_{ij,y} - \tilde{y}_{ik,y}) - (\tilde{y}_{lj,y} - \tilde{y}_{lk,y})}{2} $$
$$ \boldsymbol{r}(\sigma\{i,l;j,k\}) = (\boldsymbol{x}_{ij} - \boldsymbol{x}_{ik}) - (\boldsymbol{x}_{lj} - \boldsymbol{x}_{lk}) ,$$
\noindent where $\sigma(\cdot)$ maps an ordered quadruple to the index set $\mathcal{N}_{m_n} = \{1,2,\dots, m_n \}$, with $m_n$ denoting the number of distinct ordered quadruples from $\mathcal{N}$, i.e, $m_n = n (n-1) (n-2) (n-3)$.\footnote{\citet{jochmans2018semiparametric} exploits the permutation invariance of the score contributions in senders $(i,l)$ and receivers $(j,k)$, thus, considering combinations of senders and receivers and defining $m_n = n (n-1) (n-2) (n-3)/4$. I depart from this and consider ordered quadruples of nodes throughout. Since the score contributions are invariant under these permutations, the two formulations yield identical estimators. This choice is made for consistency with the projections of the scores defined later in this Section and with the proofs, which, in this paper, follow the ordered structure of indices.} Hereafter, I use the shortcut notation $z_{\sigma,y}$ and $\boldsymbol{r}_\sigma$.
The transformed variable $z_{\sigma,y}$ takes values in
$\{-1,-1/2,0,1/2,1\}$, with $z_{\sigma,y} \in \{-1,1\}$
corresponding to the conditioning set
$\{\tilde{y}_{ij,y} + \tilde{y}_{ik,y} = 1,\;
\tilde{y}_{lj,y} + \tilde{y}_{lk,y} = 1,\;
\tilde{y}_{ij,y} + \tilde{y}_{lk,y} = 1\}$ from
Section~\ref{subsection_estimation_method}. Collecting $\boldsymbol{x} = (\boldsymbol{x}_{ij},
\boldsymbol{x}_{ik}, \boldsymbol{x}_{lj},
\boldsymbol{x}_{lk})$, the conditioning argument from the previous section yields:
\begin{lemma}\label{lemma:sufficiency}(Sufficiency) $$ \operatorname{Pr}[z_{\sigma,y} = 1 \mid \boldsymbol{x}, z_{\sigma,y} \in \{-1,1\}] = \frac{\exp(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_{y,0})}{1 + \exp(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_{y,0})}$$
\end{lemma}
\textit{Proof.} Follows from the conditioning argument in
Subsection~\ref{subsection_estimation_method}. In particular, it follows immediately from the exponential
family structure of the logistic specification in
Equation~\eqref{eq:main_charbonneau}.
Lemma \ref{lemma:sufficiency} implies that, after conditioning, the fixed effects are eliminated from the likelihood, circumventing the incidental parameter problem, and each
informative quadruple contributes a standard logistic
term. Summing over all quadruples in $\mathcal{N}_{m_n}$, for which $z_{\sigma,y} \in \{-1,1\}$, the estimator is defined as
$$ {\boldsymbol{\theta}}_{n,y} = \operatorname*{arg\,max}_{\boldsymbol{\theta}_y \in \Theta} L_n (\boldsymbol{\theta}_y),$$
\noindent where $\Theta$ is the parameter space searched over, and the conditional log-likelihood takes the form
$$ L_n (\boldsymbol{\theta}_y) = \sum_{\sigma \in \mathcal{N}_{m_n}} 1 \{z_{\sigma,y} = 1 \} \text{log} \Lambda(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_y) + 1 \{z_{\sigma,y} = -1 \} \text{log} (1 -\Lambda(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_y)).$$
The number of quadruples with $z_{\sigma,y} \in \{-1,1\}$ is denoted by $m_{n,y}^* = \sum_{\sigma \in \mathcal{N}_{m_n}} 1 \{z_{\sigma,y} \in \{-1,1\}\}$.
The following standard assumptions are needed to establish consistency of the estimator:
\begin{assumption}(Sampling) The n nodes in $\mathcal{N}$ are sampled independently.
\label{assumption1}
\end{assumption}
\begin{assumption}(Parameter space) $\boldsymbol{\theta}_{y,0}$ is interior to $\Theta$, a compact subset of $\mathbb{R}^{dim (\boldsymbol{\theta}_y)}$.
\label{assumption2}
\end{assumption}
\begin{assumption}(Moments) For all $(i,j) \in \mathcal{D}$, $\mathbb{E} (||\boldsymbol{x}_{ij}||^2) < C_1$, where $C_1$ is a finite constant.
\label{assumption3}
\end{assumption}
Define the expected fraction of quadruples that contribute to the log-likelihood at each threshold as:
$$ p_{n,y} = \frac{\mathbb{E}(m^*_{n,y})}{m_n} = \frac{\sum_{\sigma \in \mathcal{N}_{m_n}} \Pr \{z_{\sigma,y} \in \{-1,1\}\}}{m_n}. $$
\begin{assumption} (Identification) $n p_{n,y} \xrightarrow[]{} \infty$ as $n \xrightarrow[]{} \infty$ and the matrix
$$ \lim_{n \xrightarrow[]{} \infty} (m_n p_{n,y})^{-1} \sum_{\sigma \in \mathcal{N}_{m_n}} \mathbb{E}(-\boldsymbol{r}_\sigma\boldsymbol{r}_\sigma' f(\boldsymbol{r}_\sigma'\boldsymbol{\theta}_{y,0}) 1\{z_{\sigma,y}\in \{-1,1\}\}), $$ where $f$ is the logistic density function, has maximal rank.
\label{assumption4}
\end{assumption}
Assumption \ref{assumption1} is a standard sampling scheme for network data, requiring independent sampling of nodes. Importantly, it does not impose independence across dyads, since covariates for dyads sharing a node may be correlated, which, together with the same unrestricted fixed effects appearing across different pairs sharing a node, generates the network dependence. This accommodates for settings where covariates take the form $x_{ij} = g(x_i, x_j)$ \citep{graham2017econometric}, where $g(\cdot)$ is a measurable function, but it is more general than that. Moreover, the assumption does not require that the nodes are identically distributed.
Assumptions \ref{assumption2}, \ref{assumption3} and the rank condition in Assumption \ref{assumption4} are standard for establishing consistency in non-linear models \citep{newey1994chapter}. The matrix in Assumption \ref{assumption4} is the limit of the normalized expected Hessian, so that the rank condition is equivalent to requiring the limiting Hessian to be negative definite. The key condition in Assumption \ref{assumption4} is that $np_{n,y} \to \infty$, where $p_{n,y}$ measures the expected fraction of informative quadruples at the fixed threshold $y$. In the network formation setting, for which this estimator was originally proposed, this allows $p_{n,y}$ to shrink as $n$ grows, so that the probability of link formation can approach zero or one. In the DR context, this condition allows the probabilities $\Pr(y_{ij} \leq y)$ to approach zero or one, accommodating
sparsity both from the underlying network structure and from
binarization at extreme thresholds. Since $p_{n,y}$ is the
average of $\Pr(z_{\sigma,y} \in \{-1,1\})$ over all
quadruples $\sigma$, it depends on the fixed effects of each
quadruple, so that sequences of fixed effects growing without
bound are allowed as long as $p_{n,y}$ does not shrink faster
than $n^{-1}$. The condition $np_{n,y} \to \infty$ thus
requires that the expected number of informative quadruples
continues to grow with $n$, even as the fraction of
informative quadruples may vanish. The formal
characterization of these two sources of sparsity and their
implications for the estimable range of thresholds are
developed in Subsection~\ref{subsection_sparsity}.
The following theorem establishes pointwise consistency at
each threshold.
\begin{theorem} \label{theorem1} (Consistency) Let Assumptions \ref{assumption1}-\ref{assumption4} hold. Then ${\boldsymbol{\theta}}_{n,y} \overset{p}{\to} \boldsymbol{\theta}_{y,0}$ as $n \xrightarrow[]{} \infty$ for each fixed $y \in \mathcal{Y}$.
\end{theorem}
\noindent \textit{Proof.} See \ref{appendix_asymptotics}.
Conventional logit standard errors are not valid for the estimated ${\boldsymbol{\theta}}_{n,y}$. The score vector sums over quadruples of nodes, and quadruples sharing common nodes induce dependence across summands (such that each node participates in $O(n^3)$ quadruples), so the information matrix equality does not hold, and a sandwich-type variance estimator is needed. To derive the asymptotic distribution of the estimator, the moment requirements must be
strengthened:
\begin{assumption} \label{assumption5} (Moments) For all $(i,j) \in \mathcal{D}$, $\mathbb{E} (||\boldsymbol{x}_{ij}||^6) < C_2$, where $C_2$ is a finite constant.
\end{assumption}
Each summand of the score vector takes the form
$$ \boldsymbol{s}_y(\sigma, \boldsymbol{\theta}_y) = \boldsymbol{r}_\sigma \{ 1 \{z_{\sigma,y} = 1 \} (1-\Lambda(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_y)) - 1 \{z_{\sigma,y} = -1 \} \Lambda(\boldsymbol{r}_\sigma' \boldsymbol{\theta}_y) \},$$
and the score vector is
$$ \boldsymbol{S}_{n,y} (\boldsymbol{\theta}_y) = \sum_{i}^n \sum_{j \neq i} \sum_{\substack{l \neq i,j}} \sum_{\substack{k \neq i,j,l}} \boldsymbol{s}_y(\sigma\{i,l;j,k\}; \boldsymbol{\theta}_y).$$
The asymptotic distribution of the estimator is characterized by the result that $$\boldsymbol{\Upsilon}_{n,y} (\boldsymbol{\theta}_{y,0})^{-1/2} \boldsymbol{S}_{n,y}(\boldsymbol{\theta}_{y,0}) \overset{d}{\to} N(\boldsymbol{0},\boldsymbol{I}),$$
\noindent where $\boldsymbol{\Upsilon}_{n,y}(\boldsymbol{\theta}_y)$
is defined as follows
\begin{align}
\boldsymbol{\Upsilon}_{n,y}({\boldsymbol{\theta}}_y) = \sum_{i} \sum_{j \neq i} \sum_{i' \neq i,j} \sum_{j' \neq i,j,i'} \sum_{i'' \neq i,j,i'} \sum_{j'' \neq i,j,j',i''} 16 \times \left[ \boldsymbol{s}_y(\sigma\{i,i';j,j'\}; {\boldsymbol{\theta}}_y) \boldsymbol{s}_y(\sigma\{i,i'';j,j''\}; {\boldsymbol{\theta}}_y)'\right].
\end{align}
Its expectation at
$\boldsymbol{\theta}_{y,0}$ gives the leading term of the
variance of the score vector, comprising the $O(n^6)$ terms
corresponding to pairs of quadruples that share exactly one
dyad.\footnote{The factor $16$ arises from fixing the shared
dyad $(i,j)$ to be the first sender-receiver pair in both
quadruples: there are four positions in each quadruple where
$(i,j)$ can appear, giving $4 \times 4 = 16$ ordered pairs
that share $(i,j)$ in the same position.
\citet[Supplement, p.~28]{jochmans2018semiparametric} obtains
the same result by multiplying each score contribution by
$4$.}
Theorem \ref{theorem2} also requires that the dyad-clustered score outer-product matrix be
nondegenerate at the scale of its leading term, ruling out first-order degeneracy of the score in
any parameter direction.
\begin{assumption} \label{assumption_scorevar} (Score-covariance nondegeneracy) For each fixed $y\in\mathcal Y$, there exists a constant $c_y>0$ such that
\[
\Pr\!\left\{
\lambda_{\min}\!\left[(n^6p_{n,y})^{-1}
\boldsymbol\Upsilon_{n,y}(\boldsymbol\theta_{y,0})\right]\ge c_y
\right\}\longrightarrow1.
\]
\end{assumption}
An analogous condition is required for the corresponding result in
\citet{jochmans2018semiparametric}, and \citet{muris2025dyadic} impose an analogous eigenvalue
condition on their dyad-clustered score outer-product matrix. Defining the Hessian
$$ \boldsymbol{H}_{n,y}(\boldsymbol{\theta}_y) = - \sum_{\sigma \in \mathcal{N}_{m_n}} \boldsymbol{r}_\sigma\boldsymbol{r}_\sigma' f(\boldsymbol{r}_\sigma'\boldsymbol{\theta}_y) 1\{z_{\sigma,y} \in \{-1,1\}\}$$
and the sandwich variance estimator
$$ \boldsymbol{\Omega}_{n,y}(\boldsymbol{\theta}_{n,y}) = \boldsymbol{H}_{n,y}({\boldsymbol{\theta}}_{n,y})^{-1} \boldsymbol{\Upsilon}_{n,y} ({\boldsymbol{\theta}}_{n,y})\boldsymbol{H}_{n,y}({\boldsymbol{\theta}}_{n,y})^{-1},$$
the following result holds, which establishes that pointwise (for each threshold $y$), the estimator converges to the true parameter value, and the sandwich estimator for the asymptotic variance delivers valid inference.
\begin{theorem} (Asymptotic distribution) Let Assumptions \ref{assumption1}-\ref{assumption_scorevar} hold. Then
$$|| {\boldsymbol{\theta}}_{n,y} - \boldsymbol{\theta}_{y,0} || = O_p (1/ \sqrt{n(n-1)p_{n,y}})$$ and
$$ \boldsymbol{\Omega}_{n,y}(\boldsymbol{\theta}_{n,y})^{-1/2} ({\boldsymbol{\theta}}_{n,y} - \boldsymbol{\theta}_{y,0}) \overset{d}{\to} N({0},\boldsymbol{I}) $$ as $n \xrightarrow[]{} \infty$, for each fixed $y \in \mathcal{Y}$.
\label{theorem2}
\end{theorem}
\noindent \textit{Proof.} See \ref{appendix_derivation}.
The proof follows a four-step structure, with steps akin to those used when establishing the limit distribution of a U-statistic \citep{graham2017econometric}: (i) a projection of the score vector is proposed\footnote{The proposed projection resembles a H\'{a}jek projection, but it is not formally one, since the kernel of this projection is not symmetric, and we do not only condition on observable and unobservable attributes of a specific dyad $(i,j)$, but also on node attributes.} and its asymptotic equivalence to the score evaluated at the true parameter $\boldsymbol{\theta}_{y,0}$ is established; (ii) the limit distribution of the projection is derived via a conditional CLT; (iii) the uniform convergence of the Hessian is established; and (iv) the results are combined via a mean-value expansion.
The main departure from \citet{jochmans2018semiparametric} is in step (i): while their approach explicitly computes the projection of the scores, I leverage the conditional independence structure, following \citet{graham2017econometric}, to show that the scores and projections are asymptotically equivalent. This approach extends naturally to the cross-threshold setting of Section \ref{joint_distribution}, avoiding the need to compute explicit cross-threshold joint probabilities. Moreover, I derive that the score variance is $O(n^6p_{n,y})$, which determines the convergence rate of the estimator and is used to establish the joint distribution across thresholds.
\begin{remark}[Computation]
Evaluating the objective function, score, and Hessian requires
summation over all informative quadruples, which can be
computationally costly for large networks. Appendix~2.A of
\citet{roelsgaard2025essays} describes implementation strategies for
a related conditional likelihood estimator, more specifically, a two-way duration model. The proposed implementation includes storing the covariates in contiguous memory, pre-computing the covariate differences and storing only the relevant ones (only the ones referring to informative quadruples) in a vector, speeding up access, and reusing intermediate calculations within
each Newton step, all of which apply directly to the present setting. Moreover, since
estimation at each threshold is performed independently, the
$K$ threshold-level optimizations can be run in parallel.
\end{remark}
\subsection{Sparsity Conditions and Estimable Threshold Range}\label{subsection_sparsity}
The identification condition
$np_{n,y} \to \infty$ accommodates sparsity from two distinct
sources. This subsection characterizes these sources formally,
translates the conditions of
\citet{jochmans2018semiparametric} to the DR setting, and
proposes a parameterization of fixed effects across thresholds
that captures the variation in sparsity across the
distribution.
\begin{figure}[htbp]
\centering
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\captionsetup{skip=0pt}
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{density_mild.png}
\caption{Density, mild first-degree sparsity}
\label{fig:sparsity_density_mild}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{density_strong.png}
\caption{Density, strong first-degree sparsity}
\label{fig:sparsity_density_strong}
\end{subfigure}
\vspace{0.5em}
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{cdf_mild.png}
\caption{CDF, mild first-degree sparsity}
\label{fig:sparsity_cdf_mild}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[width=\textwidth]{cdf_strong.png}
\caption{CDF, strong first-degree sparsity}
\label{fig:sparsity_cdf_strong}
\end{subfigure}
\caption{Two types of sparsity in distribution regression, for an outcome
bounded below at zero, under mild (left) and strong (right) concentration
at the bound. The mass point at zero is the source of first-degree
sparsity. Second-degree
sparsity, shown by the shading intensity in the density panels and the band
at $F(y) \approx 1$ in the CDF panels, arises where the binary indicators
$\tilde{y}_{ij,y}$ vary little across dyads. The mild and strong columns show
how a larger mass point places $F(y)$ near one over a wider range of thresholds.}
\label{fig:two_types_sparsity}
\end{figure}
\paragraph{Two Types of Sparsity in Distribution Regression} The distribution
regression framework for network data involves two distinct but related
notions of sparsity, illustrated in Figure~\ref{fig:two_types_sparsity}.
\textit{First-degree sparsity} arises when the outcome variable is bounded
and a positive fraction of observations sit at the bound, creating a mass
point at the boundary of the support. Estimation begins
at the bound, the lowest threshold at which the binary indicators
$\tilde{y}_{ij,y} = \mathbf{1}\{y_{ij} \leq y\}$ vary across dyads. When the
mass point is large, as in firm-level trade networks, patent citation
networks, or venture capital flows, where 90--99\% of dyads take the
boundary value, this first threshold already places most indicators at one,
so the bound itself lies in a sparse regime and the estimable range is
compressed to a narrow interval.
\textit{Second-degree sparsity} refers to
the proportion of zeros or ones in the binary indicators at each threshold,
which determines whether sufficient variation exists for identification. It
arises at thresholds in either tail of the conditional distribution, as
$\Pr(y_{ij} \leq y)$ approaches zero or one, including in outcomes with unbounded support, which exhibit no first-degree sparsity.\footnote{Analogously, sparsity
can be induced through few links or few non-links in a network formation
setting.} The two are related: when first-degree sparsity is extreme, it
places $\Pr(y_{ij} \leq y)$ near one across most thresholds, inducing
second-degree sparsity at nearly all thresholds. Whether sparsity arises from a mass point or from the
tails of the conditional distribution, bias correction
methods require $\Pr(y_{ij} \leq y)$ bounded away from
zero and one, a condition violated under both regimes.
The approach in this paper accommodates both through
Assumption~\ref{assumption4}: $p_{n,y} \to 0$ is
allowed provided $np_{n,y} \to \infty$.
\paragraph{Identification condition and the estimable quantile range} The
identification condition $np_{n,y} \to \infty$ can be translated into
restrictions on (i) the average probability of falling below threshold $y$,
$q_{n,y} = \sum_{i=1}^n \sum_{j \neq i} \Pr\{y_{ij} \leq y\} / n(n-1)$,
providing a direct characterization of the estimable threshold range; and
(ii) on the growth rate of the fixed effects. Following the argument in
\citet{jochmans2018semiparametric}, the exponential tails of the logistic
distribution imply $q_{n,y} \sim \sqrt{p_{n,y}}$ in the left tail and
$1 - q_{n,y} \sim \sqrt{p_{n,y}}$ in the right tail. Combined with the
condition $np_{n,y} \to \infty$, this yields rate restrictions on the
estimable thresholds:
$$
\sqrt{n}\, q_{n,y} \to \infty \text{ in the left tail,}
\qquad \sqrt{n}\,(1-q_{n,y}) \to \infty \text{ in the
right tail.}
$$
Accessing very extreme quantiles therefore requires
correspondingly larger samples, regardless of the distribution of node
heterogeneity. \ref{appendix_sparsity_rates} gives the full derivation.
\paragraph{Fixed effect parameterization across thresholds} To understand
the growth rates of fixed effects that satisfy Assumption~\ref{assumption4}
and their connection to sparsity across thresholds, I adopt the
parameterization
$$
\alpha_{i,y} = a_{n,i} \cdot g_y(n), \qquad
\gamma_{j,y} = b_{n,j} \cdot g_y(n),
$$
where $a_{n,i}, b_{n,j} \in [0,1]$ are bounded sequences in $n$ capturing
node-level heterogeneity, and $g_y(n)$ is a real-valued sequence that scales
all fixed effects proportionately, governing how sparsity varies with the
threshold $y$. Specifically, $g_y(n) < 0$ at low thresholds (sparsity through many zeros),
$g_y(n) \approx 0$ at intermediate thresholds (dense network), and
$g_y(n) > 0$ at high thresholds (sparsity through many ones). The larger
$|g_y(n)|$, the more the binary indicators concentrate near zero or one.
This parameterization describes the fixed effects at each threshold, not a
data generating process for the outcome $y_{ij}$ itself. When the outcome
is bounded below, for instance, at zero, the parameterization is defined only
for thresholds from the bound onward, where sparsity varies across
thresholds as governed by $g_y(n)$.
The sequences $a_{n,i}$ and $b_{n,j}$ capture each node's sensitivity to the changes in sparsity generated by $g_y(n)$: nodes with high values
have fixed effects that change substantially across thresholds, causing
sharp transitions from mostly zeros to mostly ones; nodes with low values
have more stable fixed effects, maintaining more stable probabilities
across thresholds. This parameterization extends the framework of
\citet{jochmans2018semiparametric} by allowing $g_y(n)$ to vary across
thresholds, capturing how sparsity changes across the distribution. This
cross-threshold structure underlies the joint asymptotic analysis of
Section~\ref{joint_distribution} and the specific parameterizations
explored in the Supplemental \ref{sim_single_threshold}.
The choice of $a_{n,i}$ and $b_{n,j}$ also affects how far into the
sparse region estimation remains feasible. The largest admissible
$|g_y(n)|$ depends on the heterogeneity pattern, and so, through the
induced link probability, does the set of estimable thresholds. Holding
the mean of $a_{n,i}$ and $b_{n,j}$ fixed, more dispersed sequences
admit larger $|g_y(n)|$. The
identification condition admits, however, a single characterization in
terms of the link probability that holds whatever the heterogeneity
pattern: a threshold is estimable when $\sqrt{n}\, q_{n,y} \to \infty$
in the left tail and $\sqrt{n}\,(1 - q_{n,y}) \to \infty$
in the right tail. The
heterogeneity pattern enters this criterion only through the value of
$q_{n,y}$ it induces at each threshold; the criterion itself is
unchanged. \ref{appendix_sparsity_heterogeneity} develops this
distinction.
\section{Joint Inference Across Thresholds}\label{joint_distribution}
This section extends the pointwise results of
Section~\ref{asymptotics} to the joint distribution of
the estimators across multiple thresholds, enabling
simultaneous confidence bands and formal tests of
coefficient equality. Subsection~\ref{joint} establishes
the joint asymptotic distribution.
Subsections~\ref{bands} and \ref{testing} develop the
inference tools.
\subsection{Joint Asymptotic Distribution Across Thresholds}\label{joint}
Since the estimators at different thresholds are computed from the same underlying outcome variable $y_{ij}$, they are not independent: informative quadruples at different thresholds are correlated since they potentially share the same dyadic outcomes, and the sequences of fixed effects and idiosyncratic errors are allowed to be correlated across thresholds. Moreover, the convergence rates differ across thresholds through the threshold-specific parameters $p_{n,y_k}$, reflecting the expected fraction of informative quadruples at each threshold and the different degrees of sparsity at different thresholds. The joint distribution established below accounts for both the cross-threshold dependence and the different convergence rates.
Let $\mathbf{y} = (y_1, ..., y_K)'$ denote the vector of $K$ ordered thresholds with $y_1 < y_2 < ... < y_K$, and each $y_k \in \mathcal{Y}$. For each threshold $y_k$, the parameter vector $\boldsymbol{\theta}_{y_k} \in \mathbb{R}^p$ represents the $p$-dimensional coefficient vector at that threshold. Define the stacked parameter vector $\boldsymbol{\theta}_{\mathbf{y}} = (\boldsymbol{\theta}_{y_1}', ..., \boldsymbol{\theta}_{y_K}')' \in \mathbb{R}^{Kp}$. Stacking the threshold-specific scores gives
$$\boldsymbol{S}_{n,\mathbf{y}}(\boldsymbol{\theta}_{\mathbf{y},0}) =
\big(\boldsymbol{S}_{n,1}(\boldsymbol{\theta}_{y_{1,0}})', ...,
\boldsymbol{S}_{n,K}(\boldsymbol{\theta}_{y_{K,0}})'\big)'.$$
Define the $(Kp) \times (Kp)$ dyad-clustered outer-product matrix
$\boldsymbol{\Upsilon}_{n, \mathbf{y}}(\boldsymbol{\theta}_{\mathbf{y}})$ with blocks
$\left[\boldsymbol{\Upsilon}_{n, \mathbf{y}}\right]_{kl} =
\boldsymbol{\Upsilon}_{n, kl}\left(\boldsymbol{\theta}_{y_k},
\boldsymbol{\theta}_{y_l}\right)$, where
\begin{adjustwidth}{-\oddsidemargin}{-\oddsidemargin}
\centering
$$\boldsymbol{\Upsilon}_{n, k l}
(\boldsymbol{\theta}_{y_k}, \boldsymbol{\theta}_{y_l})=\sum_{i=1}^n \sum_{j \neq i} \sum_{i^{\prime} \neq i, j} \sum_{j^{\prime} \neq i, j, i^{\prime}} \sum_{i^{\prime \prime} \neq i, j, i^{\prime}} \sum_{j^{\prime \prime} \neq i, j, j^{\prime}, i^{\prime \prime}} 16 \times \left[\boldsymbol{s}_k\left(\sigma\left\{i, i^{\prime} ; j, j^{\prime}\right\}; \boldsymbol{\theta}_{y_{k}}\right) \boldsymbol{s}_l\left(\sigma\left\{i, i^{\prime \prime} ; j, j^{\prime \prime}\right\}; \boldsymbol{\theta}_{y_l}\right)^{\prime}\right]$$
\end{adjustwidth}
whose expectation at $\boldsymbol{\theta}_{\mathbf{y},0}$ gives the leading term of the score
covariance for thresholds $y_k$ and $y_l$. The diagonal blocks coincide with the pointwise matrices $\boldsymbol{\Upsilon}_{n, y_k}(\boldsymbol{\theta}_{y_k})$, and the off-diagonal blocks capture the cross-threshold dependence.
The joint limit distribution is driven by that of the stacked score vector, and the first two
steps of Theorem \ref{theorem2} carry over. In step (i), the stacked projections are shown to be
asymptotically equivalent to the stacked scores by comparing their leading covariance terms,
which is where the conditional independence structure replaces the explicit cross-threshold joint
probabilities of the binary indicators, as anticipated in Section \ref{asymptotics}. In step
(ii), conditioning on the joint information set containing the covariates and all
threshold-specific fixed effects preserves cross-dyad independence, and a conditional central
limit theorem applied to the projections delivers joint asymptotic normality of
$\boldsymbol{S}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y},0})$ after normalization by its
own covariance.
Normalizing by the covariance itself absorbs the threshold-specific rates automatically, but it
requires that the covariance stays invertible as $n$ grows, and this cannot be stated without a reference scale: $\boldsymbol{\Upsilon}_{n,\mathbf{y}}$
grows without bound, at rates that differ across blocks. Define the score-rate matrix
$$\boldsymbol{D}_{n,\mathbf{y}} := \operatorname{diag}\left((n^6p_{n,y_1})^{1/2}\boldsymbol{I}_p,
\ldots, (n^6p_{n,y_K})^{1/2}\boldsymbol{I}_p\right),$$
so that the $(k,l)$ block of
$\boldsymbol{D}_{n,\mathbf{y}}^{-1}\boldsymbol{\Upsilon}_{n,\mathbf{y}}
\boldsymbol{D}_{n,\mathbf{y}}^{-1}$ is
$\boldsymbol{\Upsilon}_{n,kl}/(n^6\sqrt{p_{n,y_k}p_{n,y_l}})$, which is $O_p(1)$ for every
pair of thresholds.
A scalar normalization must be matched to one threshold or another: chosen for the sparsest, the
densest block diverges; chosen for the densest, the sparsest block vanishes and the bound below
fails for reasons of rate alone. Dividing each block by its own scale instead leaves only the
shape of the covariance, so that the condition below is a restriction on how the score behaves
across directions and not on how the $p_{n,y_k}$ compare across thresholds. The same blockwise scaling is used throughout the proof, so that no condition is required on the
relative convergence rates across thresholds.
\begin{assumption} (Joint score-covariance nondegeneracy) \label{assumption_jointscorevar}
For every fixed finite collection $\mathbf{y} = (y_1, \ldots, y_K)'$ of distinct thresholds,
there exists a constant $c_{\mathbf{y}} > 0$ such that
$$\Pr\left\{\lambda_{\min}\left[\boldsymbol{D}_{n,\mathbf{y}}^{-1}
\boldsymbol{\Upsilon}_{n,\mathbf{y}}(\boldsymbol{\theta}_{\mathbf{y},0})
\boldsymbol{D}_{n,\mathbf{y}}^{-1}\right] \geq c_{\mathbf{y}}\right\} \longrightarrow 1.$$
\end{assumption}
Assumption \ref{assumption_jointscorevar} requires that no nonzero linear combination of the
threshold-specific scores become asymptotically degenerate after this normalization. For $K = 1$ it reduces to Assumption \ref{assumption_scorevar}. Assumption \ref{assumption_jointscorevar} fails, for instance, when two thresholds are close
enough that the corresponding binary indicators nearly coincide, which suggests that the
intervals between selected thresholds should contain a non-negligible fraction of observations.
The $K p \times K p$ joint Hessian inverse
$\boldsymbol{H}_{n, \mathbf{y}}^{-1}\left(\boldsymbol{\theta}_{\mathbf{y}}\right) =
\text{diag} \left(\boldsymbol{H}_{n, 1}\left(\boldsymbol{\theta}_{y_1}\right)^{-1}, \ldots,
\boldsymbol{H}_{n, K}\left(\boldsymbol{\theta}_{y_K}\right)^{-1}\right)$ is block-diagonal since each threshold is estimated separately. The sandwich variance estimator is
$$\boldsymbol{\Omega}_{n, \mathbf{y}}\left(\boldsymbol{\theta}_{n,\mathbf{y}}\right)=\boldsymbol{H}_{n, \mathbf{y}}^{-1}\left(\boldsymbol{\theta}_{n,\mathbf{y}}\right) \boldsymbol{\Upsilon}_{n, \mathbf{y}}\left(\boldsymbol{\theta}_{n,\mathbf{y}}\right) \boldsymbol{H}_{n, \mathbf{y}}^{-1}\left(\boldsymbol{\theta}_{n,\mathbf{y}}\right),$$
where each block is given by $
\left[\boldsymbol{\Omega}_{n, \mathbf{y}}\left(\boldsymbol{\theta}_{n,\mathbf{y}}\right)\right]_{k l}=\boldsymbol{H}_{n, k}\left(\boldsymbol{\theta}_{n,y_k}\right)^{-1} \boldsymbol{\Upsilon}_{n, k l}\left(\boldsymbol{\theta}_{n,y_k}, \boldsymbol{\theta}_{n,y_l}\right) \boldsymbol{H}_{n, l}\left(\boldsymbol{\theta}_{n,y_l}\right)^{-1}
$.
For any $\boldsymbol\theta$, define
\[
\boldsymbol\sigma_{n,\mathbf y}(\boldsymbol\theta)
:=\operatorname{diag}\!\left(
\sqrt{[\boldsymbol\Omega_{n,\mathbf y}(\boldsymbol\theta)]_{11}},
\ldots,
\sqrt{[\boldsymbol\Omega_{n,\mathbf y}(\boldsymbol\theta)]_{Kp,Kp}}
\right),
\]
\[
\boldsymbol P_{n,\mathbf y}(\boldsymbol\theta)
:=\boldsymbol\sigma_{n,\mathbf y}(\boldsymbol\theta)^{-1}
\boldsymbol\Omega_{n,\mathbf y}(\boldsymbol\theta)
\boldsymbol\sigma_{n,\mathbf y}(\boldsymbol\theta)^{-1},
\]
and
\[
\boldsymbol T_{n,\mathbf y}
:=\boldsymbol\sigma_{n,\mathbf y}
(\boldsymbol\theta_{n,\mathbf y})^{-1}
(\boldsymbol\theta_{n,\mathbf y}-\boldsymbol\theta_{\mathbf y,0}).
\]
Thus $\boldsymbol P_{n,\mathbf y}(\boldsymbol\theta_{n,\mathbf y})$ is the feasible plug-in version of the sample correlation matrix. Under Assumption \ref{assumption_jointscorevar}, the diagonal entries of $\boldsymbol{\sigma}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$ are positive and $\boldsymbol{P}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$ is nondegenerate with probability approaching one (see \ref{joint_distribution_appendix}).
For probability measures $\mu, \nu$ on $\mathbb{R}^{Kp}$, let
$$d_{\mathrm{BL}}(\mu,\nu) := \sup_{f}\left|\int f \, d\mu - \int f \, d\nu\right|,$$
where the supremum is over all $f: \mathbb{R}^{Kp} \to \mathbb{R}$ with $\|f\|_\infty \leq 1$ and
$|f(\boldsymbol{u}) - f(\boldsymbol{v})| \leq \|\boldsymbol{u} - \boldsymbol{v}\|$ for all
$\boldsymbol{u}, \boldsymbol{v}$. This is the bounded--Lipschitz distance, which metrizes weak
convergence and remains meaningful when the approximating law varies with $n$. Write
$\mathcal{F}_{n,\mathbf{y}}$ for the covariates together with the fixed effects at all $K$
thresholds, and $\mathcal{L}(\boldsymbol{X} \mid \mathcal{F}_{n,\mathbf{y}})$ for the conditional
law of $\boldsymbol{X}$ given these.
\begin{theorem} (Joint asymptotic distribution) \label{theorem3} Let Assumptions
\ref{assumption1}--\ref{assumption_scorevar} and \ref{assumption_jointscorevar} hold, and fix a
finite collection $\mathbf y=(y_1,\ldots,y_K)'$ of distinct thresholds in $\mathcal Y$. Then, as
$n\to\infty$:
\begin{enumerate}
\item[(i)] The coordinatewise studentized vector $\boldsymbol T_{n,\mathbf y}$ satisfies
\[
d_{\mathrm{BL}}\!\left(
\mathcal{L}\!\left(\boldsymbol T_{n,\mathbf y} \mid \mathcal{F}_{n,\mathbf{y}}\right),\;
N\!\left(\boldsymbol 0, \boldsymbol P_{n,\mathbf y}(\boldsymbol\theta_{n,\mathbf y})\right)
\right) \overset p\longrightarrow 0.
\]
\item[(ii)] For every fixed nonzero $\boldsymbol a\in\mathbb R^{Kp}$,
\[
\frac{\boldsymbol a'
(\boldsymbol\theta_{n,\mathbf y}-\boldsymbol\theta_{\mathbf y,0})}
{\sqrt{\boldsymbol a'
\boldsymbol\Omega_{n,\mathbf y}(\boldsymbol\theta_{n,\mathbf y})
\boldsymbol a}}
\overset d\longrightarrow N(0,1).
\]
\item[(iii)] For every fixed $q\times Kp$ matrix $\boldsymbol R$ of full row rank,
\[
\begin{aligned}
&\big[\boldsymbol R(\boldsymbol\theta_{n,\mathbf y}-\boldsymbol\theta_{\mathbf y,0})\big]'
\big[\boldsymbol R\boldsymbol\Omega_{n,\mathbf y}
(\boldsymbol\theta_{n,\mathbf y})\boldsymbol R'\big]^{-1}\\
&\hspace{3.5cm}\times
\big[\boldsymbol R(\boldsymbol\theta_{n,\mathbf y}-\boldsymbol\theta_{\mathbf y,0})\big]
\overset d\longrightarrow\chi_q^2.
\end{aligned}
\]
\end{enumerate}
\end{theorem}
\noindent \textit{Proof.} See \ref{joint_distribution_appendix}.
The block structure of the sandwich covariance is what accommodates the different rates: the
Hessian blocks satisfy $\boldsymbol{H}_{n,k}(\boldsymbol{\theta}_{n,y_k}) = O_p(n^4 p_{n,y_k})$
and the score covariance blocks
$\boldsymbol{\Upsilon}_{n,kl}(\boldsymbol{\theta}_{n,y_k}, \boldsymbol{\theta}_{n,y_l}) =
O_p(n^6 \sqrt{p_{n,y_k} p_{n,y_l}})$, so that
$$\left[\boldsymbol{\Omega}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})\right]_{kl} =
O_p\left(\frac{1}{n^2\sqrt{p_{n,y_k}\, p_{n,y_l}}}\right).$$
The diagonal blocks recover the threshold-specific rates $(n(n-1)p_{n,y_k})^{-1/2}$ of Theorem
\ref{theorem2}, and the off-diagonal blocks scale with the geometric mean of the two
informativeness parameters, so that the implied cross-threshold correlations are unaffected by
the relative magnitudes of $p_{n,y_k}$ and $p_{n,y_l}$. Dividing each coefficient by its own
standard error, as in $\boldsymbol{T}_{n,\mathbf{y}}$, normalizes these rates away. What it does
not remove is the dependence across thresholds, which stays in the correlation matrix
$\boldsymbol{P}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$, and no limit is imposed on
that matrix. Doing so would be a condition on the joint informativeness structure across
thresholds, which determines the off-diagonal blocks of
$\boldsymbol{P}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$ and concerns which quadruples
are informative at two thresholds at once: a restriction on the joint design rather than on each
threshold separately.
Part (i) of Theorem \ref{theorem3} is therefore stated as a Gaussian approximation, with a
correlation matrix that varies with $n$, rather than as convergence to a fixed law. Normalizing
the whole vector instead would give a fixed limit without any such restriction, but it would
replace the individual coefficients by linear combinations of them across thresholds, which is
not the form the bands of Section \ref{bands} require; parts (ii) and (iii) take that route, and
their limits are correspondingly fixed.
\subsection{Sup-$t$ Joint Confidence Bands}\label{bands}
Simultaneous confidence bands across thresholds can be constructed via the sup-$t$ method \citep{montiel2019simultaneous}. Since typically the question of interest is on inference on each separate covariate (whether the effect of covariate $d$ varies across the distribution of the outcome), the bands are constructed per covariate.\footnote{Joint bands across all covariates and thresholds can be constructed analogously by replacing $\boldsymbol{\Omega}_{n,d}$ with $\boldsymbol{\Omega}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$ in Algorithm \ref{alg:supt_bands}; the sup-$t$ critical value adjusts automatically to the dimension of the chosen index set.} Let $\boldsymbol{\Omega}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})$ denote the $K \times K$
submatrix of $\boldsymbol{\Omega}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$
corresponding to the $d$-th coefficient across all $K$ thresholds, and let
$\sigma_{n,y_k,d} = \sqrt{[\boldsymbol{\Omega}_{n,d}]_{kk}}$, so that the
$\sigma_{n,y_k,d}$ are the corresponding diagonal entries of
$\boldsymbol{\sigma}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$ and the correlation
matrix implied by $\boldsymbol{\Omega}_{n,d}$ is the corresponding submatrix of
$\boldsymbol{P}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$, denoted by $\boldsymbol{P}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})$.
The sup-$t$ confidence band is the set
$${C}_{n,d} = \left\{ \boldsymbol{\theta}_d \in \mathbb{R}^K :
\max_{k=1,\ldots,K}
\frac{|{\theta}_{n,y_k,d} - \theta_{y_k,d}|}
{{\sigma}_{n,y_k,d}} \leq {c}_{n,d,1-\alpha} \right\}$$
where ${c}_{n,d,1-\alpha}$ is the empirical $(1-\alpha)$ quantile
of $\max_{k=1,\ldots,K} |{Z}_{n,y_k,d}| / {\sigma}_{n,y_k,d}$
with ${\boldsymbol{Z}}_{n,d} \sim N(0, {\boldsymbol{\Omega}}_{n,d})$. This is equivalent to the Cartesian product of intervals $C_{n,y_k,d} = [\theta_{n,y_k,d} \pm
c_{n,d,1-\alpha} \sigma_{n,y_k,d}]$. Unlike pointwise confidence intervals, which cover each $\theta_{y_k,d}$ individually at level
$(1-\alpha)$, the joint confidence bands cover the entire parameter vector $\boldsymbol{\theta}_d$ simultaneously at level $(1-\alpha)$. This allows for valid confidence statements that compare across thresholds, accounting for the cross-threshold dependence.
The critical value ${c}_{n,d,1-\alpha}$
generally exceeds the pointwise critical value $z_{1-\alpha/2}$.
However, as shown in \citet{montiel2019simultaneous}, the sup-$t$ band achieves the smallest critical value among the one-parameter class of simultaneous confidence bands with coverage at least $(1 - \alpha)$, which includes Bonferroni and Wald projection bands as special cases. The gains are particularly sizeable when estimates are highly correlated, as is the case for nearby thresholds in the DR setting, and the sup-$t$ critical value increases slowly with the number of thresholds $K$. Moreover, unlike confidence ellipsoids, the rectangular structure of the bands is easily visualized and permits visual hypothesis testing, regardless of the dimension of the parameter vector.
The critical value ${c}_{n,d,1-\alpha}$ is obtained by simulation, as described in Algorithm \ref{alg:supt_bands}.
\begin{algorithm}[h]
\setlength{\interspacetitleruled}{0pt}
\setlength{\interspacealgoruled}{2pt}
\setlength{\algomargin}{1em}
\fontsize{9}{11}\selectfont
\linespread{1.2}\selectfont
\caption{sup-$t$ Confidence Bands for Covariate $d$}
\label{alg:supt_bands}
\KwIn{Estimates $\theta_{n,y_k,d}$ for $k=1,\ldots,K$;
covariance matrix $\boldsymbol{\Omega}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})$;
significance level $\alpha$;
number of draws $B$}
\KwOut{Simultaneous confidence band}
\For{$b = 1, \ldots, B$}{
Draw ${\boldsymbol{Z}}^{(b)}_{n,d} \sim N(0, {\boldsymbol{\Omega}}_{n,d})$\;
Compute $T^{(b)}_{n,d} = \max_{k=1,\ldots,K}
|{Z}^{(b)}_{n,y_k,d}| / {\sigma}_{n,y_k,d}$\;
}
Set ${c}_{n,d,1-\alpha} = (1-\alpha)$ empirical quantile of
$\{T^{(1)}_{n,d}, \ldots, T^{(B)}_{n,d}\}$\;
\For{$k = 1, \ldots, K$}{
${C}_{n,y_k,d} = [\theta_{n,y_k,d} -
{c}_{n,d,1-\alpha} {\sigma}_{n,y_k,d}, \;
\theta_{n,y_k,d} +
{c}_{n,d,1-\alpha} {\sigma}_{n,y_k,d}]$\;
}
\end{algorithm}
The sup-$t$ band requires (i) a joint Gaussian approximation for the coordinatewise studentized
estimates, with (ii) the correlation structure of the approximating law consistently estimated in the
relative sense that its difference from the correlation matrix conditional on the covariates and fixed effects at all thresholds vanishes, and (iii)
strictly positive marginal variances. Theorem \ref{theorem3}(i) delivers the first two directly,
since the approximating law is built from
$\boldsymbol{P}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$, which the proof shows to
satisfy this, and the third holds because the feasible covariance matrix is positive definite with probability approaching one under
Assumption \ref{assumption_jointscorevar}; see \ref{joint_distribution_appendix}. Algorithm
\ref{alg:supt_bands} is a direct application of Algorithm 1 in
\citet{montiel2019simultaneous}. The draws enter only through the standardized maxima
$\max_k |{Z}^{(b)}_{n,y_k,d}|/{\sigma}_{n,y_k,d}$, and since the ${\sigma}_{n,y_k,d}$ are the
square roots of the diagonal entries of $\boldsymbol{\Omega}_{n,d}$, drawing from
$\boldsymbol{\Omega}_{n,d}$ and dividing by ${\sigma}_{n,y_k,d}$ is equivalent to drawing from the
implied correlation matrix $\boldsymbol{P}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})$ directly; the
former is retained to match Algorithm 1 of \citet{montiel2019simultaneous}.
Because the draws are standardized, each coordinate has unit variance and the different
convergence rates across thresholds are absorbed without restricting how they compare. The standardization also makes the anti-concentration bound for the maximum of Gaussian
coordinates \citep{chernozhukov2015comparison} hold with a constant that does not depend on the
correlation matrix, which is what converts the bounded--Lipschitz approximation of Theorem
\ref{theorem3}(i) into a statement about the probability of the rectangle defining ${C}_{n,d}$.
This step is not needed in \citet{montiel2019simultaneous}, where the correlation matrix has a
fixed limit; here it is what allows the correlation matrix to vary with $n$. The resulting coverage statement holds
conditionally and hence unconditionally, since coverage probabilities are bounded.
\subsection{Joint Testing}\label{testing}
While simultaneous confidence bands provide visual evidence on whether coefficients vary across thresholds, formal tests of the null hypothesis of coefficient equality are of interest. For a given covariate $d$, the null hypothesis is
$$H_{0,d}: \theta_{y_1,d} = \theta_{y_2,d} = \cdots
= \theta_{y_K,d},$$
which can be expressed as $\boldsymbol{R} {\boldsymbol{\theta}}_d =
{0}$, where $\boldsymbol{R}$ is any $(K-1) \times K$
matrix of full row rank satisfying $\boldsymbol{R} \mathbf{1}_K =
\boldsymbol{0}$. The choice of $\boldsymbol{R}$ does not affect the
null hypothesis, but determines the parameterization
of the deviations from equality. Common choices include adjacent differences
($\boldsymbol{R}{\boldsymbol{\theta}}_d = (\theta_{y_2,d} -
\theta_{y_1,d}, \ldots, \theta_{y_K,d} -
\theta_{y_{K-1},d})'$) and differences from a
reference threshold ($\boldsymbol{R}{\boldsymbol{\theta}}_d =
(\theta_{y_2,d} - \theta_{y_1,d}, \ldots,
\theta_{y_K,d} - \theta_{y_1,d})'$). I focus on two possible tests: the standard Wald test and a sup-$t$ test obtained by band inversion. The Wald test is
invariant to the choice of $\boldsymbol{R}$, while
the sup-$t$ test depends on the specific contrasts
used, which may affect power but not validity.
The sup-$t$ equality test is obtained by inverting a sup-$t$ confidence band for the contrast
vector $\boldsymbol{R}{\boldsymbol{\theta}}_d$ \citep{montiel2019simultaneous}. Let
$\boldsymbol{R}$ select the $d$th coefficient across thresholds and form the contrasts. The
argument of Theorem \ref{theorem3}(iii), applied to this $\boldsymbol{R}$, gives a joint Gaussian
approximation for $\boldsymbol{R}({\boldsymbol{\theta}}_{n,d} - \boldsymbol{\theta}_{d,0})$ with
covariance $\boldsymbol{R}\boldsymbol{\Omega}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})
\boldsymbol{R}'$. Studentizing each contrast by its own standard error (the square root of the
corresponding diagonal entry of that matrix), then puts the vector on the scale used by the
sup-$t$ construction, with correlation structure implied by the same matrix. Under $H_{0,d}$ the
centering vanishes, since $\boldsymbol{R}\boldsymbol{\theta}_{d,0} = \boldsymbol{0}$, so the same
approximation applies to $\boldsymbol{R}{\boldsymbol{\theta}}_{n,d}$ itself, and the test
statistic is
$$T_{n,d}^{\text{eq}} = \max_{k=1,\ldots,K-1}
\frac{|(\boldsymbol{R}{\boldsymbol{\theta}}_{n,d})_k|}
{{\sigma}_{n,\delta,k}}$$
where ${\sigma}_{n,\delta,k} =
\sqrt{[\boldsymbol{R}{\boldsymbol{\Omega}}_{n,d}
\boldsymbol{R}']_{kk}}$, and $H_{0,d}$ is rejected
if $T_{n,d}^{\text{eq}} > {c}^{\text{eq}}_{n,d,1-\alpha}$. The critical value ${c}^{\text{eq}}_{n,d,1-\alpha}$
is computed as in Algorithm \ref{alg:supt_bands},
replacing ${\boldsymbol{\Omega}}_{n,d}$ with
$\boldsymbol{R}{\boldsymbol{\Omega}}_{n,d}
\boldsymbol{R}'$. Since the contrast vector is studentized by the diagonal of
$\boldsymbol{R}\boldsymbol{\Omega}_{n,d}(\boldsymbol{\theta}_{n,\mathbf{y}})\boldsymbol{R}'$, its
approximating law again has unit diagonal, and the argument of Section \ref{bands} applies
unchanged.
As a complement, the standard Wald test for $H_{0,d}$ is valid by Theorem \ref{theorem3}(iii),
with $\boldsymbol{R}$ there replaced by the $(K-1) \times Kp$ matrix mapping
$\boldsymbol{\theta}_{\mathbf{y}}$ to the $K-1$ contrasts of the $d$th coefficient across
thresholds. The test does not depend on which contrast matrix is used, since any two full-rank matrices with
the same null space differ by a nonsingular transformation, which cancels in the quadratic form. The Wald test rejects when $(\boldsymbol{R}{\boldsymbol{\theta}}_{n,d})'
(\boldsymbol{R}{\boldsymbol{\Omega}}_{n,d}\boldsymbol{R}')^{-1}(\boldsymbol{R}
{\boldsymbol{\theta}}_{n,d})$, which under $H_{0,d}$ coincides with the centered quadratic form of
Theorem \ref{theorem3}(iii), exceeds the $(1-\alpha)$-quantile of the $\chi^2_{(K-1)}$
distribution. The two tests have complementary power properties: the sup-$t$ test has higher power against sparse alternatives (when a small number of thresholds have different coefficients), while the Wald test has higher power against diffuse alternatives (where all coefficients differ slightly across thresholds). In
particular, neither test uniformly dominates the
other in terms of local power \citep{montiel2019simultaneous}. However, the sup-$t$ framework allows immediate
identification of which thresholds drive the
rejection, either through the simultaneous bands for
$\boldsymbol{\theta}_d$ or through bands for the
contrasts $\boldsymbol{R}\boldsymbol{\theta}_d$,
whereas the Wald test provides only a joint rejection
decision.\footnote{The joint distribution established in
Theorem \ref{theorem3} also permits testing more
refined hypotheses, such as the equality of adjacent
coefficients ($H_k: \theta_{y_k,d} =
\theta_{y_{k+1},d}$) or monotonicity restrictions
($\theta_{y_1,d} \leq \cdots \leq
\theta_{y_K,d}$). I focus on the global test
and joint confidence bands, which suffice to
characterize heterogeneity across the distribution.}
\begin{remark}[Extension to a continuum of thresholds]
\label{rem:continuum}
The asymptotic results in this paper are established for
finitely many thresholds. The main payoff of extending them
to a continuum $y \in \mathcal{Y}$ would be to obtain uniform confidence
bands for counterfactual distributions, and by inversion, for quantile
functions and effects \citep{chernozhukov2013inference,
chernozhukov2020network}. In the sparse network setting, where consistent estimation of the fixed effects is not possible, the construction of these objects is unavailable. Additionally, average effects that aggregate over
the fixed effects are at best set-identified under sparsity. What remains possible is inference on the structural parameter process $y \mapsto \boldsymbol{\theta}_0(y)$ as a functional object, but the sup-$t$ test developed in
Section~\ref{joint_distribution} already provides
simultaneous inference across multiple thresholds, which suffices for the
question of whether covariate effects vary across the
distribution. Moreover, existing uniform results for dyadic and network data do not directly apply to the setting in this paper.\footnote{The empirical process theory for exchangeable arrays of \citet{davezies2021empirical} does not accommodate triangular array structures with the varying rates of convergence across thresholds. The functional central limit theorem in \citet{chernozhukov2020network} is developed for bias-corrected fixed effects estimators in dense networks where the rate of convergence is uniform across thresholds, and does not apply to the conditional maximum likelihood estimator, whose asymptotic structure (involving projections over conditionally independent quadruples with threshold-dependent rates) differs fundamentally from theirs. Among existing approaches, the
strong approximation methods of
\citet{cattaneo2024uniform} for dyadic kernel density
estimation share key structural features with the present
setting: the dyad-level contributions are conditionally
independent given node-specific attributes (and in the case of this paper, dyad-specific attributes as well), and the rate of convergence is not uniform across evaluation points (due to degeneracy of the H\'{a}jek projection in their setting, and due to varying network sparsity across thresholds in the present one). Adapting
their approach to the present setting requires establishing
that the linearization of the CMLE holds uniformly over
the continuum of thresholds (in particular, uniform
convergence of the Hessian over $\mathcal{Y}$ and uniform negligibility of
the higher-order remainder term) which is not required in
\citet{cattaneo2024uniform} because their estimator is a sample average that can be directly expressed as a sum over dyads, whereas the CMLE is defined implicitly as the solution to a score equation. One should also verify that the strong
approximation error vanishes uniformly over $\mathcal{Y}$
when the sparsity parameter $p_{n,y}$ varies continuously
across thresholds.} Developing this extension is left for
future research.
\end{remark}
\section{Monte Carlo Simulations}\label{simulations}
This section presents Monte Carlo simulation studies evaluating the finite-sample performance of the CMLE, relative to the maximum likelihood estimator (MLE) and the analytical bias correction method (BC, \citet{chernozhukov2020network}). I consider two sets of simulations: (i) for the full distribution regression setting calibrated to the international trade application; and (ii) for joint inference across thresholds, including sup-$t$ confidence bands and equality tests. Appendix ... additionally provides a Monte Carlo simulation for a single threshold (a network formation model) with varying degrees of sparsity, extending the DGP of \citet{jochmans2018semiparametric} with additional node heterogeneity specifications and right-tail sparsity (the left-tail case is studied in their paper). The results confirm that the CMLE maintains smaller bias than BC under extreme sparsity, consistent with the theoretical predictions of Section~\ref{asymptotics}.
\subsection{Monte Carlo Simulations for the Distribution Regression}\label{monte_carlo_DR}
Using the bilateral trade dataset \citep{helpman2008estimating,jochmans2018semiparametric,chernozhukov2020network} described in Section \ref{application}, which is standard in the gravity model literature, a logit model with two-way fixed effects is estimated by MLE at each threshold $y_k$, and the estimates are set as the true parameters in the simulation. Following the notation of Section~\ref{model_estimation},
the parameter vector at the $k$-th threshold is denoted
$\boldsymbol{\theta}_{y_k}$. When thresholds correspond
to empirical quantiles, $\tau$ denotes the quantile level
associated with threshold $y$, so that
$\theta_{y,d} = \theta_d(\tau)$ for the $d$-th covariate. At each threshold, the true DGP is
$$\Pr(\tilde{y}_{ij,k} = 1 \mid \boldsymbol{x}_{ij},
\alpha_{i,y_k}^{\text{MLE}}, \gamma_{j,y_k}^{\text{MLE}}) =
\Lambda\big(\boldsymbol{x}_{ij}'
\boldsymbol{\theta}_{y_k}^{\text{MLE}}
+ \alpha_{i,y_k}^{\text{MLE}}
+ \gamma_{j,y_k}^{\text{MLE}}\big), \quad
k = 1, \dots K, \quad (i,j) \in \mathcal{D},$$
where $\boldsymbol{x}_{i j}$ are the values of the covariates for the observational unit $(i, j)$ in the trade data set, $\tilde{y}_{ij,k} = 1 (y_{ij} \leq y_k)$, and the thresholds $y_1 < \dots < y_K$ correspond to the
empirical quantile levels $\tau_1 < \dots < \tau_K$ of the trade outcome at
$\tau_k \in \{0.545, 0.550, \ldots, 0.990\}$.
The true parameter vector includes all coefficients and fixed effects:
$$\boldsymbol{\beta}_k^{\text{MLE}}=\left(
\boldsymbol{\theta}_{y_k}^{\text{MLE}},
\alpha_{1,y_k}^{\text{MLE}}, \ldots,
\alpha_{n,y_k}^{\text{MLE}},
\gamma_{2,y_k}^{\text{MLE}}, \ldots,
\gamma_{n,y_k}^{\text{MLE}}\right).$$
The simulated data are generated by drawing a single vector of $n(n-1)$
logistic errors $\varepsilon_{ij} \sim \text{Logistic}(0,1)$, shared
across all thresholds, so that the binary indicators at different thresholds are generated from
common shocks, as they are in real distribution regression data. The parameters remain
threshold-specific: $\boldsymbol{\theta}_{y_k}^{\text{MLE}}$, $\alpha_{i,y_k}^{\text{MLE}}$ and
$\gamma_{j,y_k}^{\text{MLE}}$ are estimated separately at each threshold, so the design imposes
no restriction on how they vary across the distribution, which is what DR is intended to
accommodate. This cross-threshold dependence is accounted for by the joint inference procedures
in Section~\ref{joint_distribution}. The collection of binary variables at the different thresholds is generated as
$$\tilde{y}_{ij,k}^{sim} = {1}\big(\boldsymbol{W}_{ij}'\boldsymbol{\beta}_k^{\text{MLE}} + \varepsilon_{ij} \geq 0\big), \quad k = 1, \ldots, K, \quad(i, j) \in \mathcal{D},$$
where $\boldsymbol{W}_{ij}$ includes covariates and dummy variables for the dimensions $i$ and $j$ (fixed effects indicators).\footnote{This design choice differs from the Tobit-based DGP of \citet{chernozhukov2020network}. In their DGP, a single latent variable, together with a single sequence of fixed effects (estimated by Tobit on the underlying continuous outcome), generates the binarized outcomes at all thresholds. Since the fixed effects do not vary with $y$,
link probabilities become increasingly homogeneous
across nodes at extreme thresholds. The threshold-specific design adopted here allows the fixed effects to reflect the data at each threshold separately, capturing the cross-threshold
heterogeneity patterns observed in the trade data.}
The dataset contains information on $n=157$ countries, and the covariates included are log of distance, common legal system, contiguous border, common language, and common religion.\footnote{Three additional covariates (colonial ties, currency union, and regional free trade agreement) are available in the dataset but excluded from the analysis because each takes the value one for fewer than 2\% of country pairs (Table~\ref{tab:descriptive} in the Supplemental Appendix), resulting in insufficient variation in the pairwise-differenced covariates at many thresholds.} Results are based on 500 simulations.
\paragraph{Results}
The average probability $\Pr(\tilde{y}_{ij,k} = 1)$ at each
threshold increases approximately linearly from around 0.55
at the 54.5th percentile to nearly 1.0 at the 99th
percentile, recovering the empirical quantile structure of
the trade data on which the DGP is calibrated
(Figure~\ref{fig:sparsity} in the Supplemental \ref{simulation_appendix}).
\begin{figure}
\centering
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\captionsetup{skip=0pt}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_bias_plot1.png}
\caption{Mean bias: Log distance}
\label{fig:sim_bias_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_bias_plot2.png}
\caption{Mean bias: Legal system}
\label{fig:sim_bias_legal}
\end{subfigure}
\vspace{0.15cm}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_biasmed_plot1.png}
\caption{Median bias: Log distance}
\label{fig:sim_medbias_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_biasmed_plot2.png}
\caption{Median bias: Legal system}
\label{fig:sim_medbias_legal}
\end{subfigure}
\caption{Mean bias (top) and median bias (bottom) of the
bias-corrected estimator (BC) and the conditional maximum
likelihood estimator (CMLE) across quantiles of the trade
distribution, based on 500 replications using the empirically
calibrated DGP. Thresholds correspond to empirical quantiles
in [0.545, 0.990] at intervals of 0.005.}
\label{fig:sim_bias}
\label{fig:sim_medbias}
\end{figure}
Figure~\ref{fig:sim_medbias}
plots the mean and median bias of BC and CMLE across
quantiles for log distance and common legal system; results
for contiguous border, common language, and common religion
are reported in Figures \ref{fig:sim_bias_appendix} and \ref{fig:sim_medbias_appendix} in the Supplemental Appendix. For most of the quantile range, both estimators show negligible bias. In a narrow window before the more extreme quantiles, the BC estimator shows a slightly smaller bias for most covariates. This pattern reflects the fact that the BC uses all observations to estimate the bias correction, while the CMLE uses only the informative quadruples, which constitute a smaller effective sample at moderate sparsity levels. At the most extreme thresholds, in particular beyond the 98th percentile, the BC estimator shows sharply increasing bias for most covariates, while the CMLE's bias remains substantially smaller. On median bias, which is robust to outlier replications, CMLE has smaller bias than BC at extreme thresholds across all covariates. This pattern is consistent with the single-threshold
simulations in Tables~\ref{tab:est_jochmans_right_uniform_full} and \ref{tab:est_jochmans_right_beta_full} in the Supplemental \ref{sim_single_threshold}.
\begin{figure}
\centering
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\captionsetup{skip=0pt}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_RMSE_plot1.png}
\caption{RMSE: Log distance}
\label{fig:sim_rmse_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_RMSE_plot2.png}
\caption{RMSE: Legal system}
\label{fig:sim_rmse_legal}
\end{subfigure}
\vspace{0.15cm}
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_Size_plot1.png}
\caption{Size: Log distance}
\label{fig:sim_size_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{plots_pdr_colonyfta_nonan/pdr_colonyfta_Size_plot2.png}
\caption{Size: Legal system}
\label{fig:sim_size_legal}
\end{subfigure}
\caption{RMSE (top) and rejection frequency of the two-sided
$t$-test at the 5\% nominal level (bottom) for the
bias-corrected estimator (BC) and the conditional maximum
likelihood estimator (CMLE) across quantiles of the trade
distribution, based on 500 replications using the empirically
calibrated DGP. Thresholds correspond to empirical quantiles
in [0.545, 0.990] at intervals of 0.005.}
\label{fig:sim_rmse}
\label{fig:sim_size}
\end{figure}
Figure~\ref{fig:sim_rmse} reports the RMSE (top) and
rejection frequencies (bottom) across quantiles. Through
most of the quantile range, the RMSE of both estimators is
nearly identical. In the narrow range before the extremes, the CMLE's RMSE increases somewhat more than BC's,
reflecting the higher variance from relying on only the
informative quadruples. The RMSE comparison at the extreme
quantiles depends on which effect dominates. For log
distance, where BC's bias is large, the CMLE has
substantially smaller RMSE. For common legal system, where
BC's bias is more moderate, the two estimators have
comparable RMSE. For some covariates for which results are in the Appendix (notably
contiguous border), the CMLE's variance from the smaller
effective sample outweighs its bias advantage, giving BC
smaller RMSE. This bias-variance trade-off is structural:
the CMLE uses only observations that identify $\boldsymbol{\theta}_{y_k}$,
which gives smaller bias at the cost of higher variance.
Turning to size control, for log distance, both estimators
maintain adequate rejection frequencies through most of the
quantile range. At the extreme quantiles, the BC rejection
rate rises sharply, reaching over 0.50 at the 99th
percentile, as the bias dominates the test statistic. The
CMLE remains close to the nominal level throughout. For
common legal system, both estimators fluctuate around the
nominal level, with both becoming conservative at some of
the highest quantiles. Across covariates, BC displays over-rejection at the extremes
when the bias dominates (log distance, contiguous border) and
under-rejection when the standard error estimate becomes
inflated (common language at the extreme). The CMLE generally maintains better size control,
though it too becomes conservative for some covariates at
the most extreme quantiles.
\subsection{Monte Carlo Simulations for Joint Confidence Bands and Testing}\label{monte_carlo_simultaneous}
The simultaneous confidence bands and equality tests across
thresholds developed in Section~\ref{joint_distribution} are
evaluated using the empirically calibrated DGP described in
the previous Subsection, with parameters varied across
specifications to test specific null and alternative
hypotheses.
The simulations vary along three dimensions: the number of thresholds $K \in \{5, 15, 50\}$, the location of the thresholds in the distribution (capturing different sequences of convergence rates of threshold-specific estimators), and the values of the parameters. For the location of the thresholds, three configurations are considered: (i) \textit{equally spaced}, with $K$ thresholds spread across the full range of available quantiles, from $\tau = 0.545$ to $\tau = 0.95$ (with intervals of 0.005), capturing variation across different sparsity levels; (ii) \textit{first $K$}, which considers the $K$ lowest thresholds, where the network is relatively dense ($\Pr(\tilde{y}_{ij,k} = 1)$ is moderate); and (iii) \textit{last $K$}, which considers the $K$ highest thresholds, where the network is sparse ($\Pr(\tilde{y}_{ij,k} = 1)$ is very high).
The different designs for the values of the parameters are
distinguished by the specification of the structural
parameter vector $\boldsymbol{\theta}_{y_k}$ across thresholds;
in both designs, the fixed effects $\alpha_{i,y_k},
\gamma_{j,y_k}$ are as in the previous Subsection and remain
threshold-specific. Under \textit{constant
parameters}, the parameter vector is held constant
across thresholds at $\boldsymbol{\theta}_{y_k} = \boldsymbol{\theta}_{y_1}^{\text{MLE}}$ for all $k$, so that the null hypothesis $H_{0,d}: \theta_{y_1,d} = \cdots = \theta_{y_K,d}$ holds and the empirical rejection rates should be close to the
nominal $5\%$ level. Under \textit{varying parameters},
the coefficients are set to the empirical estimates
$\boldsymbol{\theta}_{y_k} = \boldsymbol{\theta}_{y_k}^{\text{MLE}}$, which
vary across thresholds, so that $H_{0,d}$ is false.
For each simulation and configuration, I compute the CMLE at each of
the $K$ thresholds, estimate the joint covariance matrix
$\boldsymbol{\Omega}_{n,\mathbf{y}}(\boldsymbol{\theta}_{n,\mathbf{y}})$, and obtain sup-$t$ critical values from 100000 draws from the estimated
Gaussian distribution. The equality tests use deviations from the first
threshold as the contrast specification, i.e.,
$\delta_{y_k,d} = \theta_{y_k,d} - \theta_{y_1,d}$ for $k = 2, \ldots, K$ and all covariates $d$. Results are based on $500$ simulations.
\paragraph{Results}
Table \ref{tab:main_coverage_095} reports the empirical coverage of the sup-$t$ joint confidence bands and the empirical simultaneous coverage of the pointwise confidence intervals (the fraction of simulations where all $K$ pointwise intervals simultaneously contain the true values) for the configuration considering equally spaced thresholds, for constant and varying parameters, with the highest threshold considered at
$\tau_{\max} = 0.95$. Supplemental \ref{simulation_appendix} provides the simulation results for the remaining designs and specifications considered above and extends the threshold range
to $\tau_{\max} = 0.99$, assessing whether the sup-$t$ bands maintain correct coverage
when thresholds reach into the sparse tail where fewer
informative quadruples are available.
\begin{table}[htbp]
\centering
\caption{Coverage Comparison: Pointwise vs Sup-$t$ Bands ($\tau_{\max} = 0.95$, Equally Spaced)}
\label{tab:main_coverage_095}
\begin{threeparttable}
\fontsize{9}{11}\selectfont
\begin{tabular}{ll cc cc cc cc cc}
\toprule
& & \multicolumn{2}{c}{Distance} & \multicolumn{2}{c}{Legal} & \multicolumn{2}{c}{Border} & \multicolumn{2}{c}{Language} & \multicolumn{2}{c}{Religion} \\
\cmidrule(lr){3-4} \cmidrule(lr){5-6} \cmidrule(lr){7-8} \cmidrule(lr){9-10} \cmidrule(lr){11-12}
$K$ & DGP & PW & Sup-$t$ & PW & Sup-$t$ & PW & Sup-$t$ & PW & Sup-$t$ & PW & Sup-$t$ \\
\midrule
5 & Varying & 83.60 & 97.40 & 81.80 & 96.80 & 85.40 & 96.00 & 81.20 & 96.60 & 83.20 & 95.20 \\
5 & Constant & 81.40 & 95.20 & 83.60 & 97.40 & 82.00 & 94.80 & 82.80 & 94.40 & 85.60 & 97.00 \\
\addlinespace
15 & Varying & 69.00 & 96.60 & 69.00 & 96.80 & 73.20 & 95.40 & 66.80 & 95.60 & 69.60 & 96.20 \\
15 & Constant & 68.00 & 94.80 & 71.40 & 96.20 & 71.80 & 95.80 & 68.00 & 95.20 & 70.80 & 97.60 \\
\addlinespace
50 & Varying & 56.20 & 96.60 & 56.20 & 97.00 & 61.60 & 96.00 & 54.40 & 95.40 & 56.20 & 96.00 \\
50 & Constant & 56.40 & 96.20 & 58.60 & 97.40 & 57.80 & 95.00 & 55.00 & 95.20 & 59.80 & 97.20 \\
\bottomrule
\end{tabular}
\begin{tablenotes}
\smallskip\fontsize{9}{11}\selectfont
\item \textit{Notes:} Empirical coverage rates (in \%) of 95\% simultaneous confidence bands. ``PW'' (Pointwise) uses $z_{0.975} = 1.96$ at each threshold; ``Sup-$t$'' uses simulated critical values. Target coverage is 95\%. ``Constant'' = coefficients identical across thresholds ($H_0$ true); ``Varying'' = coefficients follow empirical trade data pattern ($H_0$ false). Based on 500 Monte Carlo simulations.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table \ref{tab:main_coverage_095} shows that while the results for the pointwise confidence intervals show under-coverage that worsens with $K$, the sup-$t$ bands maintain coverage at or near the $95\%$
level. The Supplemental \ref{app_monte_carlo_simultaneous} shows that this
pattern holds across $K$, covariates, $\tau_{\max}$, and DGP
design, with sup-$t$ coverage remaining near the nominal level
throughout. This pattern shows the need for simultaneous
inference when constructing coverage statements across multiple
thresholds: pointwise intervals constructed at the marginal
$95\%$ level under-cover when interpreted simultaneously,
while the sup-$t$ critical value uses the joint distribution
to deliver correct simultaneous coverage. For pointwise
coverage, extending to $\tau_{\max} = 0.99$ further reduces
simultaneous coverage substantially.
\begin{table}[htbp]
\centering
\caption{Size and Power by Covariate: Sup-$t$ vs Wald ($\tau_{\max} = 0.95$, Equally Spaced)}
\label{tab:main_by_covariate_095}
\begin{threeparttable}
\fontsize{9}{11}\selectfont
\begin{tabular}{l cc cc cc cc cc}
\toprule
& \multicolumn{2}{c}{Distance} & \multicolumn{2}{c}{Legal} & \multicolumn{2}{c}{Border} & \multicolumn{2}{c}{Language} & \multicolumn{2}{c}{Religion} \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5} \cmidrule(lr){6-7} \cmidrule(lr){8-9} \cmidrule(lr){10-11}
$K$ & Sup-$t$ & Wald & Sup-$t$ & Wald & Sup-$t$ & Wald & Sup-$t$ & Wald & Sup-$t$ & Wald \\
\midrule
\multicolumn{11}{l}{\textit{Panel A: Constant Coefficients (Size, target: 5\%)}} \\
\addlinespace
5 & 4.60 & 3.80 & 2.80 & 2.60 & 3.60 & 3.00 & 4.20 & 4.20 & 3.60 & 3.80 \\
15 & 3.20 & 3.20 & 3.60 & 1.60 & 4.00 & 5.00 & 3.60 & 3.40 & 3.20 & 3.80 \\
50 & 2.80 & 2.60 & 2.80 & 1.60 & 5.80 & 13.80 & 3.60 & 2.20 & 2.40 & 2.40 \\
\midrule
\multicolumn{11}{l}{\textit{Panel B: Varying Coefficients (Power)}} \\
\addlinespace
5 & 100.00 & 100.00 & 100.00 & 99.80 & 99.40 & 99.00 & 39.80 & 40.00 & 52.80 & 60.40 \\
15 & 100.00 & 100.00 & 100.00 & 100.00 & 98.80 & 97.80 & 58.80 & 53.40 & 77.80 & 91.20 \\
50 & 100.00 & 99.60 & 100.00 & 100.00 & 99.00 & 98.20 & 54.40 & 79.00 & 79.20 & 92.40 \\
\addlinespace
\bottomrule
\end{tabular}
\begin{tablenotes}
\smallskip\fontsize{9}{11}\selectfont
\item \textit{Notes:} Rejection rates (in \%) for testing $H_0: \theta_d(\tau_1) = \cdots = \theta_d(\tau_K)$ at the 5\% level. Sup-$t$ = Sup-$t$ test; Wald = Wald test. ``Constant'' = coefficients identical across thresholds ($H_0$ true); ``Varying'' = coefficients follow empirical trade data pattern ($H_0$ false). Based on 500 Monte Carlo simulations.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table~\ref{tab:main_by_covariate_095} reports size and power
of the sup-$t$ and Wald equality tests for the equally
spaced configuration at $\tau_{\max} = 0.95$. Under the null, the sup-$t$ test stays at or near the nominal
$5\%$ level across all values of $K$, with rejection rates
between $2.4\%$ and $5.8\%$, indicating size control without
substantial over-rejection but some conservativeness in
particular cases. The Wald
test, by contrast, shows size distortion that varies by
covariate and worsens with $K$: for Border at $K = 50$, it
over-rejects at 13.8\%, while for Legal it under-rejects at
1.6\%. The sup-$t$ test exhibits high power for Distance,
Legal, and Border, consistent with the large heterogeneity in the coefficients for these covariates, shown in
Figure~\ref{fig:coef_paths_main} in the Supplemental \ref{app_monte_carlo_simultaneous}. Language and Religion show
more moderate power, reflecting smaller heterogeneity in their coefficients across thresholds. The Wald test shows higher power than the sup-$t$ test for Language and Religion
at large $K$, but at the cost of substantial size distortion
for Border at the same $K$; the sup-$t$ test provides more
uniform performance across covariates. The Supplemental \ref{app_monte_carlo_simultaneous}
shows that these patterns hold at $\tau_{\max} = 0.99$.
\section{Application to gravity models of international trade}\label{application}
There are two important features that models for bilateral international trade should take into account. First, the outcome of interest (the volume of bilateral trade) is bounded below at zero, with a mass at zero: in the dataset used in this application, approximately 55\% of country pairs have zero bilateral trade in a given year. Second, consistent estimation typically requires controlling for country-specific terms that capture unobservable barriers each country faces with all its trading partners. Such terms are known in the international trade literature as multilateral resistance terms \citep{anderson2003gravity}, and are typically treated as two-way (exporter and importer) fixed effects.
Several approaches to handling these features focus on the conditional mean of trade. The Poisson pseudo-maximum likelihood (PPML) estimator \citep{silva2006log} retains zero observations and delivers consistent estimates of the gravity coefficients in two-way fixed effects gravity specifications \citep{fernandez2016individual}. Another approach is the Heckman-type sample selection model of \cite{helpman2008estimating}, which models the decision to trade and the volume of trade jointly via a two-stage procedure with a first-stage selection equation and a second-stage equation for log positive trade.\footnote{\cite{helpman2008estimating} estimate the first-stage selection equation by standard probit, which is subject to the incidental parameter problem in two-way fixed effects settings. Bias-corrected estimators for binary outcome models with two-way fixed effects are available \citep{fernandez2016individual}, and in the dataset used here (where approximately 45\% of country pairs have positive trade), where there is no first-degree sparsity, they could be applied. In settings with sparser first-stage outcomes, the fixed effects cannot be consistently estimated even with bias correction, making the Heckman-type selection correction infeasible. \cite{sakamoto2024dyadic} develops a semiparametric sample selection estimator for dyadic data extending
\cite{kyriazidou1997estimation}, but this approach requires panel data with time-varying covariates, which is not available in this application since key gravity determinants (such as distance) are time-invariant.} However, both methods estimate effects on the conditional mean, and neither characterizes how covariates might affect different parts of the trade distribution differently.
There are theoretical reasons to expect that the effects of trade determinants are not constant across the distribution of bilateral trade flows \citep{novy2013international, bas2017micro, carrere2020gravity}. Several econometric methods have been used to study this heterogeneity. Quantile regression on log positive trade \citep{baltagi2016estimationmain} documents that gravity-covariate effects vary across quantiles but drops zero-trade observations entirely. More recently, \cite{bergstrand2025quantile} use censored quantile regression that retains zeros. However, they address the incidental parameter problem from the high-dimensional fixed effects by imposing a parametric restriction on the unobserved heterogeneity via a Chamberlain-Mundlak correlated random effects parameterization \citep{abrevaya2008effects}. Expectile regression provides an alternative distributional approach \citep{bergstrand2025tails} that extends Poisson pseudo-maximum likelihood to recover expectile-specific coefficients across the conditional distribution. This approach handles zeros and documents heterogeneous effects, but does not inherit the robustness of Poisson-based estimation to the incidental parameter problem, with no bias correction currently available.\footnote{Their Monte Carlo simulations show coverage
degradation at extreme expectiles, particularly with shorter time
dimensions. The authors leave the development of bias corrections
for future work.}
The distribution regression framework developed in this paper is used to study how the effects of standard gravity covariates vary across the conditional distribution of bilateral trade. At each threshold $y_k$, the binary indicator $\tilde{y}_{ij,k} = \mathbf{1}\{y_{ij} \leq y_k\}$ is well-defined whether or not country pair $(i,j)$ has positive trade, and the estimator is consistent without imposing a parametric specification on the unobserved heterogeneity. The joint inference framework developed in Section~\ref{joint_distribution} allows formal tests of whether covariate effects are constant across the full threshold range, rather than pairwise comparisons at selected points.
The dataset, used previously by \cite{helpman2008estimating}, \cite{jochmans2018semiparametric}, and \cite{chernozhukov2020network}, contains information on bilateral trade flows and covariates for 157 countries in 1986, where $i$ and $j$ index each country as an exporter and an importer, respectively. Descriptive statistics are reported in Table~\ref{tab:descriptive} in the Supplemental Appendix. The outcome $y_{ij}$ is the volume of trade in thousands of constant 2000 U.S. dollars from country $i$ to country $j$.\footnote{The bilateral trade flows data are from Feenstra's ``World Trade Flows, 1970--1992,'' transformed to constant 2000 U.S. dollars using the U.S. CPI by \cite{helpman2008estimating}.} The covariates $\boldsymbol{x}_{ij}$ include the following bilateral determinants of trade flows: the logarithm of distance between capitals, and binary indicators for shared legal system, contiguous border, common language, and common religion.
\begin{figure}[htbp]
\centering
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\captionsetup{skip=0pt}
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{application_plots_tau099/band_ldist.png}
\caption{Log distance}
\label{fig:band_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{application_plots_tau099/band_legal.png}
\caption{Legal system}
\label{fig:band_legal}
\end{subfigure}
\caption{Distribution regression estimates of the effect of log distance and
legal system on bilateral trade. Solid blue: CMLE; dotted red: BC; dash-dotted
gray: MLE. Shaded region and dashed lines show the 95\% simultaneous sup-$t$
confidence band and pointwise confidence intervals for the CMLE, respectively.
Thresholds correspond to empirical quantiles in $[0.545, 0.990]$ at intervals
of $0.005$.}
\label{fig:bands_main}
\end{figure}
\begin{figure}[htbp]
\centering
\captionsetup[subfigure]{font=footnotesize, skip=0pt}
\captionsetup{skip=0pt}
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{application_plots_tau099/differences/diff_h40_ldist.png}
\caption{Log distance}
\label{fig:diff_ldist}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\textwidth]{application_plots_tau099/differences/diff_h40_legal.png}
\caption{Legal system}
\label{fig:diff_legal}
\end{subfigure}
\caption{Differences $\theta_{n,d}(\tau) - \theta_{n,d}(\tau - 0.20)$ from the CMLE
estimator for log distance and legal system. Shaded region and dashed lines show
the 95\% simultaneous sup-$t$ confidence band and pointwise confidence intervals
for the differences, respectively. The horizontal dashed line marks zero.}
\label{fig:diffs_main}
\end{figure}
\begin{table}[htbp]\label{table_application}
\centering
\caption{Joint tests of coefficient equality across thresholds}
\label{tab:equality_tests}
\fontsize{9}{11}\selectfont
\begin{tabular}{l cc c cc c cc}
\toprule
& \multicolumn{2}{c}{$\tau_{\max} = 0.95$} & & \multicolumn{2}{c}{$\tau_{\max} = 0.97$} & & \multicolumn{2}{c}{$\tau_{\max} = 0.99$} \\
\cmidrule{2-3} \cmidrule{5-6} \cmidrule{8-9}
Covariate & Sup-$t$ & Wald & & Sup-$t$ & Wald & & Sup-$t$ & Wald \\
\midrule
Log distance & $3.83{}^{**}$ & $82.7$ & & $3.83{}^{**}$ & $88.8$ & & $3.83{}^{**}$ & $98.2$ \\
Legal system & $6.75{}^{**}$ & $141.3{}^{**}$ & & $6.75{}^{**}$ & $150.4{}^{**}$ & & $6.75{}^{**}$ & $156.5{}^{**}$ \\
Border & $3.95{}^{**}$ & $91.2$ & & $3.95{}^{**}$ & $95.2$ & & $3.95{}^{**}$ & $104.9$ \\
Language & $2.91$ & $79.1$ & & $2.91$ & $80.3$ & & $2.91$ & $81.1$ \\
Religion & $2.68$ & $90.7$ & & $2.68$ & $95.6$ & & $2.68$ & $96.5$ \\
\midrule
\textit{K} & \multicolumn{2}{c}{82} & & \multicolumn{2}{c}{86} & & \multicolumn{2}{c}{90} \\
\bottomrule
\end{tabular}
\begin{minipage}{\textwidth}
\smallskip\fontsize{9}{11}\selectfont
\textit{Notes:} The sup-$t$ statistic tests $H_0\colon \theta_{y_1,d} = \cdots = \theta_{y_K,d}$ using the
supremum of standardized deviations from the first threshold, with critical values
obtained from $100{,}000$ draws of the estimated Gaussian distribution implied by the joint
asymptotic distribution (Theorem~\ref{theorem3}). The Wald statistic uses
adjacent differences $R = R_{\text{adj}}$ and is distributed as $\chi^2_{K-1}$ under $H_0$.
${}^{**}$~Rejection at the $5\%$ level. $\tau_{\max}$ denotes the largest quantile included in the test.
\end{minipage}
\end{table}
Figure~\ref{fig:bands_main} displays the DR coefficient estimates
for log distance and legal system. A positive DR
coefficient indicates a negative effect on trade volume: a
covariate that increases the probability of trade falling below
a threshold $y$ reduces the likelihood of large bilateral flows. The CMLE coefficient $\theta_d(y)$ has a direct interpretation as a log-odds ratio for $\Pr(Y_{ij} \leq y \mid \boldsymbol{x}_{ij}, \nu_i, \omega_j)$, such that variation across thresholds characterizes how this log-odds effect varies across the distribution.
The CMLE estimates show substantial heterogeneity in the
structural parameters across the distribution. For log
distance, the coefficient increases from approximately $1.05$
at the 55th percentile to around $1.26$ at the 90th percentile, and $1.78$ at the 99th percentile. The increasing coefficient indicates that distance is a stronger barrier at the upper
end of the trade distribution.
For legal system, the coefficient becomes increasingly negative
across the distribution, suggesting that shared legal origins
facilitate trade more strongly among the largest bilateral
relationships.
Table~\ref{tab:equality_tests} confirms that the null of
constant coefficients is rejected for both covariates.
Turning to the comparison across estimators, for log distance, the CMLE estimates are systematically lower than both the BC and uncorrected MLE across most of the distribution, with a persistent gap of approximately $0.3$ in the area where the network is denser (until approximately the 95th percentile). Part of this gap may reflect higher-order incidental parameter bias in the BC estimator, or differences in pseudo-true parameters under misspecification of the logistic link \citep{white1982maximum, hughes2026jackknife}. In the upper tail of log distance, the uncorrected MLE diverges from the other two, while both BC and CMLE become noisier.
The pattern across the other four covariates is qualitatively different: BC and CMLE track each other closely in the interior of the distribution, then diverge meaningfully in the upper tail, where the estimates also become considerably noisier. For legal system, the divergence emerges from approximately the 93rd percentile; for border, language, and religion (Supplemental \ref{application_appendix}), it emerges earlier, from approximately the 85th-90th percentile, with larger gaps at extreme thresholds. This pattern across four of five covariates is consistent with the BC bias correction relying on the assumption of a dense network, which fails at extreme thresholds. As expected, the joint sup-$t$ confidence bands are wider than the pointwise intervals, and both widen toward the extremes as the effective sample size for estimation decreases.
Table~\ref{tab:equality_tests} reports the sup-$t$ and Wald
test statistics for $H_0\colon \theta_{y_1,d} = \cdots = \theta_{y_K,d}$, evaluated at three choices of the maximum considered
quantile: $\tau_{\max} \in \{0.95, 0.97, 0.99\}$. These vary how far the analysis extends into the sparse tail. The sup-$t$
test rejects at 5\% level for log distance, legal system, and border for all specifications, and does not reject the null for language and religion. The Wald test does not reject for any covariate apart from legal system. With $K - 1$ ranging from 81 to 89 restrictions, localized departures from equality contribute relatively little to the joint Wald statistic, while the sup-$t$ test is driven by the largest standardized deviation across thresholds. The rejection decisions are stable across $\tau_{\max}$: for sup-$t$, the supremum is attained at or below $\tau = 0.95$ for all covariates, so no rejection is overturned despite the larger critical value at higher $\tau_{\max}$; for the Wald test, the test statistic itself changes with $\tau_{\max}$ but the rejection decisions are unaffected.
To characterize where in the distribution the variation in the distribution regression coefficients occurs, Figure~\ref{fig:diffs_main} displays simultaneous 95\% confidence bands for the differences $\theta_{n,d}(\tau) - \theta_{n,d}(\tau - 0.20)$ for log distance and legal system. These bands use the feasible joint covariance of the CMLE and the Gaussian approximation of Theorem~\ref{theorem3}(i), applied to the contrasts. Since the bands have simultaneous coverage at the 95\% level, any threshold $\tau$ where the band excludes zero refers to a rejection of $H_0\colon \theta_d(\tau) = \theta_d(\tau - 0.20)$, with family-wise error rate controlled across all thresholds jointly. For log distance, the differences are mostly positive but the band excludes zero only at the two lowest threshold pairs, indicating that few 20-percentile intervals show a significant change. Combined with the rejection in Table~\ref{tab:equality_tests}, this indicates that thresholds far apart in the distribution differ significantly, but neighboring thresholds do not: the effect of distance varies gradually across the distribution. For legal system, by contrast, the band excludes zero across much of the upper range, consistent with shared legal origins becoming increasingly important for the largest bilateral relationships.
Supplemental \ref{application_appendix} presents additional results. The coefficient paths with sup-$t$ bands and pointwise confidence intervals are shown for the remaining covariates. Pointwise comparisons between the CMLE and BC for all covariates show that CMLE intervals are wider throughout, reflecting the efficiency cost of conditioning that yields the CMLE's robustness in sparse settings. Sup-$t$ bands evaluated at $\tau_{\max} = 0.97$ yield similar conclusions, and the simultaneous bands for the coefficient differences are reported for all covariates at both $\tau_{\max} = 0.97$ and $0.99$.
The CMLE does not deliver estimates of the fixed effects, which limits the recovery of the conditional distribution and the computation of counterfactual distributions or related functionals that aggregate over the fixed effects. Moreover, under sparsity, average (marginal) effects are at best set-identified, analogous to short panel data settings.\footnote{Methods for bounding partially identified average effects have been developed for panel data models with fixed effects \citep{honore2006bounds, chernozhukov2013average, davezies2025identification, botosaru2024adversarial, dobronyi2024identification, pakel2023bounds, aguirregabiria2024identification}; extending these methods to network models (under a dyadic structure and sparsity) remains an open question.} Nevertheless, beyond the log-odds interpretation of individual coefficients, ratios of coefficients at a given threshold correspond to ratios of partial derivatives of the conditional quantile function $Q(u \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)$, where $u \in (0,1)$ denotes the quantile index \citep{chernozhukov2020network}:
\begin{align}
\left.\frac{\theta_{\ell,y}}{\theta_{k,y}}\right|_{y=Q(u \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)} = \frac{\partial_{x_{ij}^{\ell}} Q(u \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)}{\partial_{x_{ij}^k} Q(u \mid \boldsymbol{x}_{ij}, \boldsymbol{\nu}_i, \boldsymbol{\omega}_j)}.
\end{align}
This relationship applies only at thresholds above the zero mass, where the conditional distribution is continuous. The ratios measure the relative importance of different covariates at each point in the distribution, even when absolute marginal effects are not point-identified. For instance, the ratio $\theta_{n,\text{ldist},y} / \theta_{n,\text{legal},y}$ indicates how much larger the marginal effect of distance is relative to that of legal system on the conditional quantile of trade at threshold $y$.
The findings in this section provide evidence of heterogeneity in covariate effects across the conditional distribution of trade. The directional patterns differ from earlier evidence based on quantile regression \citep{baltagi2016estimationmain, carrere2020gravity} and expectile regression \citep{bergstrand2025tails}, which find covariate effects typically larger at the lower part of the conditional distribution of trade flows. The bias-corrected distribution regression estimates of \cite{chernozhukov2020network} on the same data, by contrast, deliver coefficient paths qualitatively similar to those obtained with the approach proposed in this paper. The contrast across methods reflects that DR, QR, and expectile regression target different distributional objects, although differences in dataset and time period across studies may also contribute.
\section{Conclusion}\label{conclusion}
I develop a framework for distribution regression in network settings with two-way fixed effects that are treated as incidental parameters. The conditional maximum likelihood method of \citet{charbonneau2017multiple} and \citet{jochmans2018semiparametric}, originally proposed for network formation models to address the incidental parameter problem, is extended to a distribution regression setting. The estimator is applied at multiple thresholds of the outcome distribution (after a binarization of the outcome variable). The proposed method provides asymptotically unbiased pointwise estimates, in particular in sparse settings and regions of the distribution (extreme quantiles), filling a gap in the literature. I establish joint asymptotic results across multiple thresholds with heterogeneous convergence rates, and construct simultaneous sup-$t$ confidence bands and equality tests. Monte Carlo simulations confirm that the proposed estimator has reduced bias and valid inference in sparse settings, in particular at the extreme quantiles of the distribution, where methods based on analytical bias corrections are not well suited. The simulations also
confirm that the sup-$t$ bands provide correct simultaneous coverage across all configurations considered, while the Wald test exhibits size distortions that worsen with the
number of considered thresholds. An empirical application to international trade documents heterogeneity in the structural parameters of the gravity model across the distribution of trade flows and identifies where it is most pronounced.
The method proposed in this paper complements the bias-corrected distribution regression of \citet{chernozhukov2020network}. The two approaches target different objects: in dense settings, the bias-corrected DR recovers counterfactual distributions and functionals that aggregate over the fixed effects, while the DR with CMLE delivers asymptotically unbiased estimates of the structural parameter $\theta_0(y)$ in sparse settings (both in terms of first and second degree) where the fixed effects cannot be consistently estimated. Developing informative bounds on partially identified average effects in this setting, which are typically only set-identified under sparsity, is a natural direction for future research.
While the application in this paper is to international trade, the framework applies broadly to dyadic network settings with two-way fixed effects. Examples include bilateral migration \citep{anderson2011gravity, beine2016practitioners, grogger2011income}, foreign direct investment, bilateral patent flows, and bilateral trade in services. In these settings, researchers have been interested in the relative importance of different covariates; for instance, in migration, the relative importance of geographic, cultural, and policy-related barriers \citep{ortega2013effect, grogger2011income}; in trade, the tariff equivalents of various frictions. The ratio interpretation introduced in Section~\ref{application} provides a formal tool for such questions. More specifically, in bilateral migration, the ratio of the coefficient on a visa waiver indicator to the coefficient on log distance can be interpreted as the geographic distance offset by a visa waiver between origin and destination countries. Moreover, the framework can also be adapted to bipartite settings, such as worker-firm matching with wage outcomes. Notably, these applications typically feature sparse
networks with substantial mass at zero, precisely the
setting where the framework developed in this paper
provides valid inference.
\begin{comment}
Nonetheless, extending the asymptotic results to hold uniformly over a continuum
of thresholds remains a technically interesting question. Two natural approaches
from the empirical process literature face obstacles in this setting. The empirical
process results for exchangeable arrays developed by \citet{davezies2021empirical}
establish Donsker-type functional central limit theorems under conditions on the
function class that mirror those for i.i.d.\ data. However, their framework
assumes a fixed jointly exchangeable distribution, whereas the present setting
involves triangular arrays with heterogeneous incidental parameters whose
distributions change with $n$. Moreover, their results yield either a
$\sqrt{n}$~rate (non-degenerate case) or an $n$~rate (degenerate case), with no
provision for the intermediate, threshold-dependent rate $\sqrt{n^2 p_{n,y}}$ that
arises under sparsity. The functional central limit theorem approach of
\citet{chernozhukov2024network}, which operates in $\ell^\infty(\mathcal{Y})$ and
handles discrete outcomes, requires the network of binarized outcomes to be dense
at every threshold --- that is, $p_{n,y}$ bounded away from zero for all $y \in
\mathcal{Y}$ --- which is precisely the condition relaxed in this paper.
A more promising route would adapt the strong approximation methods of
\citet{cattaneo2024inference}, who establish Yurinskii couplings for dyadic kernel
density estimators allowing for varying degeneracy across evaluation points. Their
conditional-on-node-attributes coupling strategy is well suited to the present
setting, where the dyad-level score contributions are conditionally independent
given the covariates and fixed effects. The degenerate Haj\'ek projection in the
pairwise differencing estimator further simplifies the coupling architecture
relative to their general framework, eliminating the need for the KMT
approximation and the Vorob'ev--Berkes--Philipp gluing step. The main technical
challenge lies in establishing that the linearization of the conditional
maximum likelihood estimator --- specifically, the uniform convergence of the
Hessian and the uniform negligibility of the Taylor remainder --- holds over the
continuum of thresholds under sparsity, a step that is not required in
\citet{cattaneo2024inference} because their estimator is an explicit linear
functional of the data. Developing this extension is left for future research.
2 other alternatives:
A strong approximation approach similar to that of \citet{cattaneo2024inference} seems to more naturally suits this setting. In spite of the different object of interest, since \citet{cattaneo2024inference} considers the problem of uniformly estimating a dyadic Lebesgue density function using nonparametric kernel-based estimators taking the form of dyadic empirical processes, the setting is simular to the one in this paper in the sense that the dyad-level contributions are conditionally independent and the rate of
convergence varies across evaluation points. However, adapting the strong approximation approach of \citet{cattaneo2024inference} from nonparametric kernel density estimation to
conditional maximum likelihood estimation in the present setting requires
substantial additional work and is left for future research.
\end{comment}
\bibliography{library}