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.
58,006 characters
Robust Inference for Dyadic Data with Dependent Ordered Nodes
\maketitle
\begin{abstract}
{Dyadic regression models are commonly analyzed under the conventional dyadic dependence framework, where two observations may be dependent only if the corresponding dyads share a node. This paper studies inference when nodes are ordered and nearby nodes are exposed to common latent shocks, so that dyads with no shared endpoint may still be dependent. Although each additional covariance term may be weak, the number of nearby-node dyad pairs grows with the sample size, making their aggregate contribution asymptotically non-negligible. We develop an inferential framework for dyadic arrays with ordered-node dependence and propose two variance estimators: a dependent-node dyadic cluster-robust variance estimator that retains covariance terms between dyads with nearby endpoints, and a row-column moving-block jackknife method that deletes adjacent blocks of nodes together with all dyads touching those nodes. We establish the asymptotic validity of both procedures under weak dependence along the ordered node index. Monte Carlo evidence shows improvements in size control, with the jackknife procedure displaying comparatively stable finite-sample performance. An application to international trade gravity regressions shows that accounting for ordered-node dependence substantially weakens the statistical evidence for free trade agreement effects.}
\end{abstract}
\noindent\textbf{Keywords:} dyadic data, locally dependent, cluster-robust variance estimation, jackknife.\\
\textbf{JEL Classification:} C12, C15, C21, C31.
\newpage
\section{Introduction}
Dyadic data arise when an observation is attached to a pair of units. Examples
include trade between two countries, conflict between two states, financial
exposure between two banks, collaboration between two firms, and links in a
social network. A central feature of such data is that observations sharing a
node are generally dependent. For example, trade flows involving the same
country may be correlated because of country-specific shocks, and links
involving the same individual may be correlated because of individual
heterogeneity. This observation motivates the conventional dyadic
cluster-robust variance estimator, which keeps covariance terms between dyads
that share at least one endpoint.
This conventional dyadic asymptotic paradigm implicitly imposes a sparse
dependency graph: two dyads may be dependent when they share a node, but dyads
with no common endpoint are treated as asymptotically independent. This
restriction is natural in dissociated or exchangeable dyadic arrays, but it can
be too restrictive when nodes are ordered and nearby nodes are themselves
dependent. Suppose, for example, that bilateral trade flows are analyzed using a dyadic regression.
Standard dyadic inference allows dependence between Saudi Arabia-Japan and Saudi Arabia-South Korea because the two dyads share Saudi Arabia, but it treats Saudi Arabia-Japan and Kuwait-South Korea as asymptotically independent because they share no endpoint. This restriction can be implausible in many applications. The two dyads may both be affected by common oil-market shocks, global energy demand, shipping disruptions, or changes in the macroeconomic conditions of high-income importing economies. More generally, dyads that do not share a country may still be dependent when their endpoint countries are close along an economically meaningful dimension.
The ordering of the nodes is treated as given throughout the paper. {The ordering need not correspond to a physical ordering such as time or geography. It only needs to represent a meaningful one-dimensional proximity structure along which node-level dependence decays.} This is
appropriate in applications where the ordering is determined by an exogenous
and observable characteristic. For example, in gravity applications, countries can be ordered by GDP per capita, market size, or trade exposure, so that countries at similar levels of development or global-market integration are allowed to have more strongly dependent dyadic shocks. More generally, nodes may be ordered by a substantive dimension that governs dependence: firms by technological proximity, banks by balance-sheet characteristics, and individuals by cohort, location, or network position. When the ordering is estimated from the same data used
for inference, additional first-stage uncertainty may arise. Extending the
theory to estimated orderings is an important topic for future work.
This paper studies dyadic regression inference under ordered-node dependence.
We model node-level shocks as a weakly dependent process indexed by the ordered
node labels. Consequently, two dyads may be dependent not only when they share a
node, but also when one endpoint of the first dyad is close to one endpoint of
the second dyad. The resulting dependence graph is substantially denser than the
standard dyadic dependency graph. The key asymptotic phenomenon is that
conventional dyadic clustering omits an entire class of covariance terms.
Although individual omitted covariance terms may be weak, their aggregate
contribution is asymptotically non-negligible because the number of
nearby-node dyad pairs diverges with the sample size. Hence, the asymptotic
variance is no longer representable by conventional dyadic clustering.
We show that the conventional dyadic asymptotic framework fundamentally breaks down under ordered-node dependence and develop a new inferential framework for
this broader class of dyadic arrays. After a first-order projection, the
leading component of the dyadic score behaves like a weakly dependent sequence
indexed by nodes. Valid inference must therefore account simultaneously for
shared-node dependence and local dependence along the ordered node index.
We propose two variance estimators. The first is a dependent-node dyadic
cluster-robust variance estimator, abbreviated as the DN-Dyadic CRVE. It
retains covariance terms between dyads whose endpoint nodes are close in the
ordered-node metric. The second is a row-column moving-block jackknife procedure,
abbreviated as the JK-DN-Dyadic CRVE. It deletes adjacent blocks of nodes and
removes all dyads touching the deleted block. This deletion rule provides a natural dyadic analog of a moving-block jackknife because each dyadic observation is
attached to two endpoint nodes.
The paper contributes to the literature on dyadic and network inference.
Important contributions to dyadic, multiway clustered, and exchangeable-array
inference include, e.g., \citet{cameron2011robust}, \citet{thompson2011simple},
\citet{aronow2015cluster}, \citet{tabord2019inference},
\citet{menzel2021bootstrap}, \citet{davezies2021empirical}, and
\citet{davezies2025analytic}. Related work on network formation and sparse
network asymptotics includes \citet{fafchamps2007formation} and
\citet{graham2024sparse}. These papers provide tools for important dyadic and
network settings, but the conventional dyadic clustering logic is based on
exact node overlap. Our setting differs because the node labels carry an
ordering, and nearby nodes can generate additional dependence between dyads that
do not share an endpoint.
The paper is also closely related to \citet{jochmans2026two}, who studies
non-exchangeable dyadic data with dependence that decays over an ordered index
distance. The distinction is useful to make explicit. \citet{jochmans2026two}
constructs an estimator using an estimated first-order node projection. By
contrast, our DN-Dyadic CRVE is written directly in terms of dyadic regression
scores and dyad-pair covariance terms, which makes explicit which covariance
terms are added relative to conventional dyadic clustering. We also develop a
row-column moving-block jackknife procedure, motivated by the two-endpoint structure of
dyadic observations, and show that it provides significantly improved finite-sample performance
in the simulations. In addition, our theory covers a degenerate Gaussian case
in which the first-order node projection does not contribute. {The ordered-node framework gives rise to two distinct asymptotic regimes. When the first-order node projection is nondegenerate, the estimator converges at the \(\sqrt{n}\) rate, where \(n\) denotes the number of nodes. This rate reflects the effective node-level dependence induced by the ordered-node structure. In contrast, when the first-order projection is degenerate, the node-level component vanishes, and the convergence rate increases to \(n\), with the leading stochastic variation driven by the dyad-level component.
}
The paper is connected more broadly to recent work on clustered inference with
serial or local dependence. \citet{chiang2023standard},
\citet{chen2023fixed}, and \citet{hounyo2024wild} study two-way clustered regressions with serially
correlated time effects. Although their setting is not dyadic, the motivation
is related: exact cluster membership may not fully capture dependence when one
dimension is ordered. Our method also builds on the literature on jackknife
cluster-robust inference, including \citet{hansen2022jackknife},
\citet{mackinnon2023leverage}, \citet{mackinnon2024jackknife}, and \citet{hounyo2025jackknife}. The distinctive feature here is the row–column deletion rule, which removes a block of nodes together with all dyads attached to those nodes. Regression estimators based on dyadic data with dependent ordered nodes naturally lend themselves to this novel jackknife procedure.
We establish the asymptotic validity of the DN-Dyadic and JK-DN-Dyadic CRVEs
under standard moment and weak-dependence conditions. In the nondegenerate case,
both estimators consistently estimate the long-run variance generated by the
ordered node-level projection. In the degenerate Gaussian case, they adapt to
the dyad-level source of variation. Monte Carlo evidence illustrates that
conventional dyadic clustering can over-reject when ordered-node dependence is
present, while the proposed methods, especially the jackknife version, deliver
more reliable size control. {An empirical application to international trade gravity regressions further shows that accounting for ordered-node dependence can substantially weaken the statistical evidence for free trade agreement effects on bilateral manufacturing trade flows.}
The remainder of the paper is organized as follows. Section \ref{sec:model}
introduces the dyadic regression model and ordered-node dependence. Section
\ref{sec:estimators} defines the DN-Dyadic CRVE and the JK-DN-Dyadic CRVE.
Section \ref{sec:theory} presents the asymptotic validity results. Section
\ref{sec:simulation} reports the simulation evidence. Section \ref{sec:empirical_trade} illustrates the practical relevance of the
proposed approach through an empirical application. Section
\ref{sec:conclusion} concludes. Proofs are collected in the Appendix.
\section{Model and Dependence Structure}
\label{sec:model}
\subsection{Dyadic regression}
Let $i,j\in\{1,\ldots,n\}$ index nodes. We observe undirected dyadic data, with one observation for each unordered pair
\[
\mathcal{D}_{n}
=
\{(i,j):1\leq i<j\leq n\},
\qquad
M_n=\left\lvert \mathcal{D}_{n}\right\rvert=\frac{n(n-1)}{2}.
\]
For each dyad $(i,j)\in\mathcal{D}_{n}$, consider the linear regression
model
\begin{equation}
y_{ij}=x_{ij}'\beta+u_{ij},
\label{eq:model}
\end{equation}
where $x_{ij}\in\mathbb{R}^{K}$ includes a constant, $\beta\in\mathbb{R}^{K}$
is the parameter of interest, and $u_{ij}$ is the regression disturbance.
The dimension $K$ is fixed. Stacking observations over
$(i,j)\in\mathcal{D}_{n}$ gives
\[
y=X\beta+u,
\]
where $X$ is the $M_n\times K$ matrix of regressors. The OLS estimator is
\begin{equation}
\widehat{\beta}=(X'X)^{-1}X'y.
\label{eq:ols}
\end{equation}
Let
\(
\widehat{u}_{ij}=y_{ij}-x_{ij}'\widehat{\beta},
\
\widehat{s}_{ij}=x_{ij}\widehat{u}_{ij},
\)
and, for the population score, write
\(
s_{ij}=x_{ij}u_{ij}.
\)
Throughout the paper, the score $s_{ij}$ is a $K$-dimensional vector. We
focus on inference for a fixed scalar contrast $a'\beta$, where
$a\in\mathbb{R}^{K}$ is nonzero and does not depend on $n$. Given a
variance estimator $\widehat{V}$ for $\widehat{\beta}$, the corresponding
$t$ statistic is
\(
\widehat{t}
=
\frac{a'(\widehat{\beta}-\beta_{0})}
{\sqrt{a'\widehat{V}a}}.
\)
The OLS estimator satisfies the usual score expansion
\begin{equation}
\widehat{\beta}-\beta
=
(X'X)^{-1}\sum_{(i,j)\in\mathcal{D}_{n}}x_{ij}u_{ij}
=
(X'X)^{-1}\sum_{(i,j)\in\mathcal{D}_{n}}s_{ij}.
\label{eq:ols_score_expansion}
\end{equation}
Thus, the dependence structure relevant for inference is the dependence
structure of the dyadic score array $\{s_{ij}:(i,j)\in\mathcal{D}_{n}\}$.
\subsection{Ordered-node dependence}
The conventional dyadic dependence assumption allows two dyadic scores
$s_{ij}$ and $s_{pq}$ to be dependent only when the two dyads share at
least one endpoint, that is, when
\[
\{i,j\}\cap\{p,q\}\neq \varnothing .
\]
This assumption is natural when the dyadic observations are dissociated
after conditioning on independent node-specific latent variables. In many
applications, however, nodes have a meaningful order. For example, the node
index may represent time, geography along a line, birth cohort, firm rank,
or another ordering along which nearby nodes are more similar than distant
nodes.\footnote{Throughout the paper, the ordering is treated as given. This covers
settings in which the ordering is determined by an exogenous observable
characteristic, such as GDP per capita, trade exposure, time, cohort, or a pre-specified ranking.
If the ordering is estimated from the same data used for inference, additional
first-stage uncertainty may affect the limiting distribution. We leave a formal
treatment of estimated orderings to future work.} In such settings, two dyads may be dependent even when they do not
share a node.
To accommodate this feature, we allow the latent node variables to be
weakly dependent over the ordered node index. We describe the dependence structure using a latent-variable representation
in the spirit of the Aldous-Hoover-Kallenberg (AHK, \citet{aldous1981representations};
\citet{hoover1979relations}; \citet{kallenberg1989representation}) representation, but adapted
to ordered weakly dependent node variables.
\begin{assumption}[Ordered-node dyadic representation]
\label{ass:representation}
For each $(i,j)\in\mathcal{D}_{n}$,
\begin{equation}
(y_{ij},x_{ij},u_{ij})
=
h(Z_i,Z_j,Q_{ij}),
\label{eq:kernel_rep}
\end{equation}
where $\{Z_i:i\geq 1\}$ is a strictly stationary weakly dependent sequence,
$\{Q_{ij}:1\leq i<j\}$ are i.i.d. dyad-level shocks, and
$\{Q_{ij}:1\leq i<j\}$ is independent of $\{Z_i:i\geq 1\}$.
The function $h$ is symmetric in its first two arguments in the sense
needed for undirected dyadic observations.
\end{assumption}
Assumption~\ref{ass:representation} is an ordered-node version of the
usual latent-variable representation for dyadic data. The difference is
that the node-level variables $\{Z_i\}$ are not required to be independent.
If $\{Z_i\}$ is independent across $i$, then dyads with no common endpoint
are independent conditional on their node variables, and the model reduces
to the usual dissociated dyadic setting. If instead $\{Z_i\}$ is locally dependent over the ordered node index, then two dyads can be
dependent even when they do not share a node. For example, the scores
$s_{ij}$ and $s_{pq}$ may be correlated through dependence between $Z_i$
and $Z_p$, between $Z_i$ and $Z_q$, between $Z_j$ and $Z_p$, or between
$Z_j$ and $Z_q$.
The relevant notion of distance between two dyads is therefore the minimum
distance between their endpoint nodes. For dyads
$d=(i,j)$ and $d'=(p,q)$, define
\begin{equation}
\Delta(d,d')
=
\Delta\bigl((i,j),(p,q)\bigr)
=
\min\{
|i-p|,\ |i-q|,\ |j-p|,\ |j-q|
\}.
\label{eq:endpoint_distance}
\end{equation}
Conventional dyadic clustering keeps
only dyad pairs with $\Delta(d,d')=0$. Ordered-node dependence also generates
covariance terms for dyad pairs with $0<\Delta(d,d')\leq L$. For any fixed
local neighborhood, the number of such dyad pairs grows with $n$. Thus, even
when each individual covariance is small, the aggregate contribution of these
nearby-endpoint covariance terms can remain first order. This is why the
standard dyadic variance formula is not generally valid under ordered-node
dependence.
The ordered-node representation implies a useful decomposition of the
dyadic score. Define
\(
\mu = E[s_{ij}],
\)
where stationarity makes the expectation independent of $(i,j)$. Let \(F\) denote the common marginal distribution of the node variable. The first-order node projection is defined as \begin{equation} \gamma_i = \int E[s_{ij}\mid Z_i,Z_j=z]\,dF(z)-\mu, \label{eq:gamma_def} \end{equation} where \(j\ne i\) denotes a generic node index and the integral is taken with respect to the marginal law \(F\), rather than the conditional law of \(Z_j\) given \(Z_i\).
Define the second-order node interaction
\begin{equation}
\xi_{ij}=\xi(Z_i,Z_j)
=
E[s_{ij}\mid Z_i,Z_j]
-
\gamma_i
-
\gamma_j
-
\mu,
\label{eq:xi_def}
\end{equation}
and the dyad-level residual component
\begin{equation}
\zeta_{ij}
=
s_{ij}-E[s_{ij}\mid Z_i,Z_j].
\label{eq:zeta_def}
\end{equation}
Then, for each $(i,j)\in\mathcal D_n$,
\begin{equation}
s_{ij}
=
\mu+\gamma_i+\gamma_j+\xi_{ij}+\zeta_{ij}.
\label{eq:score_decomp}
\end{equation}
This decomposition is the node-dependent dyadic analogue of a Hoeffding projection, but its
interpretation differs from the conventional dyadic or two-way clustered case.
The first-order node component \(\gamma_i+\gamma_j\) captures the contribution
of node-level heterogeneity to the score. Because the ordered nodes may be
dependent, \(\gamma_i\) and \(\gamma_j\) are not independent in general.
Moreover, \(\xi_{ij}\) is the second-order component associated with the pair
of node variables \((Z_i,Z_j)\). Under ordered-node dependence, this component
may remain correlated with the first-order node component, unlike in the
standard independent-node Hoeffding decomposition. Finally, \(\zeta_{ij}\)
denotes the residual dyad-specific component after conditioning on
\((Z_i,Z_j)\). The components satisfy the following properties:
\[
E[\gamma_i]=0,\qquad
\int\xi(Z_i,Z_j=z)dF(z)=0,\qquad
\int\xi(Z_i=z,Z_j)dF(z)=0,\qquad
E[\zeta_{ij}\mid Z_i,Z_j]=0.
\]
\paragraph{Example 1:} Consider a simple example with \(E[Z_i]=0\), \(E[Z_i^2]=1\), \(E[Z_i^3]\neq 0\), and ordered-node dependence satisfying \(E[Z_j\mid Z_i]=\rho^{|i-j|} Z_i\) with \(\rho\neq 0\). Let \(s_{ij}=Z_i+Z_j+Z_iZ_j\). Then \(\mu=0\), and the marginal-projection definition gives \(\gamma_i=\int (Z_i+z+Z_i z)\,dF(z)=Z_i\) and \(\gamma_j=Z_j\). The second-order component is \(\xi_{ij}=Z_iZ_j\). It is degenerate with respect to marginal integration because, for fixed \(Z_i\), \(\int \xi_{ij}\,dF(Z_j)=Z_i\int z\,dF(z)=0\). However, under the true dependent joint law, \(E[\gamma_i\xi_{ij}]=E[Z_i^2Z_j] =E[Z_i^2E(Z_j\mid Z_i)]=\rho^{|i-j|} E[Z_i^3]\neq 0\). Thus, although \(\xi_{ij}\) is marginally degenerate, it need not be orthogonal to the first-order node component under ordered-node dependence.
Therefore, the first-order node projection is the leading component of the
average score whenever it is nondegenerate. Summing
\eqref{eq:score_decomp} over all dyads gives
\begin{align}
\frac{1}{M_n}\sum_{(i,j)\in\mathcal{D}_{n}}s_{ij}
&=
\mu
+
\frac{1}{M_n}\sum_{(i,j)\in\mathcal{D}_{n}}(\gamma_i+\gamma_j)
+
\frac{1}{M_n}\sum_{(i,j)\in\mathcal{D}_{n}}(\xi_{ij}+\zeta_{ij})
\nonumber \\
&=
\mu
+
\frac{2}{n}\sum_{i=1}^{n}\gamma_i
+
\frac{1}{M_n}\sum_{(i,j)\in\mathcal{D}_{n}}(\xi_{ij}+\zeta_{ij}).
\label{eq:average_score_decomp}
\end{align}
The second equality follows because each node appears in exactly $n-1$
dyads and
\[
\frac{1}{M_n}\sum_{(i,j)\in\mathcal{D}_{n}}(\gamma_i+\gamma_j)
=
\frac{n-1}{M_n}\sum_{i=1}^{n}\gamma_i
=
\frac{2}{n}\sum_{i=1}^{n}\gamma_i.
\]
Equation~\eqref{eq:average_score_decomp} is central. It shows that the
average dyadic score behaves, to first order, like an average of the
ordered node-level process $\{\gamma_i\}$. Therefore, if $\{\gamma_i\}$ is
locally dependent over $i$, the asymptotic variance of the OLS estimator
depends on the long-run covariance of the node projection. The remaining
terms $\xi_{ij}$ and $\zeta_{ij}$ are of smaller order under the
nondegenerate first-order projection condition imposed below.
\section{Variance Estimators}
\label{sec:estimators}
This section defines the two variance estimators studied in the paper.
Both estimators are designed for dyadic data with ordered-node dependence.
The first estimator is a sandwich-form variance estimator that keeps covariance
terms between dyads whose endpoint nodes are close. We call it the
dependent-node dyadic CRVE, abbreviated as DN-Dyadic CRVE. The second
estimator is a row-column moving-block jackknife analog. We call it the
dependent-node dyadic jackknife CRVE, abbreviated as JK-DN-Dyadic CRVE.
\subsection{Dependent-node dyadic CRVE}
\label{subsec:dn_dyadic_crve}
The usual dyadic CRVE is based on the sparse dyadic dependency graph in which
two dyads are neighbors only when they share a node. Equivalently, it assigns a
nonzero weight to the pair of dyads $(i,j)$ and $(p,q)$ only when
\(
\{i,j\}\cap\{p,q\}\neq\varnothing .
\)
Under ordered-node dependence, this graph is misspecified. Dyads with no common
endpoint may still have correlated scores when one endpoint of the first dyad is
close to one endpoint of the second dyad. Therefore, the variance estimator must
enlarge the dyadic neighborhood from exact endpoint overlap to nearby endpoint
overlap. This is not only a finite-sample correction: the omitted covariance
terms accumulate asymptotically because the number of nearby-endpoint dyad
pairs diverges with the number of nodes.
Under the dependent-node dyadic framework, $(i,j)$ and $(p,q)$ can be dependent when $i$ is
close to $p$ or $q$, or when $j$ is close to $p$ or $q$, even if the two dyads have no
endpoint in common. For dyads $d=(i,j)$ and $d'=(p,q)$, recall that $\Delta(d,d')$ in \eqref{eq:endpoint_distance} extends the usual dyadic-neighborhood relation to the
ordered-node setting.
Let
\begin{equation}
k_L(h)=\left(1-\frac{|h|}{L}\right)_{+}
=
\begin{cases}
1-|h|/L, & |h|<L,\\
0, & |h|\geq L,
\end{cases}
\label{eq:bartlett}
\end{equation}
be the Bartlett kernel, where $L$ denotes the bandwidth or block length. The DN-Dyadic CRVE meat is
\begin{equation}
\widehat{\Sigma}_{\mathrm{DN}}
=
\sum_{(i,j)\in\mathcal{D}_n}
\sum_{(p,q)\in\mathcal{D}_n}
k_L\!\left(\Delta\bigl((i,j),(p,q)\bigr)\right)
\widehat{s}_{ij}\widehat{s}_{pq}'.
\label{eq:dn_dyadic_meat}
\end{equation}
The corresponding variance estimator for $\widehat{\beta}$ is
\begin{equation}
\widehat{V}_{\mathrm{DN}}
=
(X'X)^{-1}\widehat{\Sigma}_{\mathrm{DN}}(X'X)^{-1}.
\label{eq:dn_dyadic_crve}
\end{equation}
Although \eqref{eq:dn_dyadic_meat} does not require the subtraction term
appearing in conventional dyadic CRVE, this does not by itself guarantee positive semidefiniteness. The reason is that $\Delta\bigl((i,j),(p,q)\bigr)$ is not a usual linear distance on one index. It is a minimum over four endpoint distances. Such a minimum distance over endpoints can destroy positive semidefiniteness.\footnote{{In the simulations and empirical application, non-positive semidefiniteness occurs rarely and does not materially affect inference. Standard eigenvalue-adjustment techniques may nevertheless be applied in practice if desired.}}
A useful way to interpret the DN-Dyadic estimator is through the projection decomposition of the dyadic score. In the nondegenerate case, the leading term is the first-order node projection. {So conceptually, the DN-Dyadic estimator extends HAC variance estimation to dyadic arrays by replacing temporal distance with an endpoint-distance metric between dyads.} In the degenerate case, the first-order node projection is absent, and the leading variation comes from the residual dyad-level component. {The same estimator therefore accommodates two asymptotic regimes:} it estimates a long-run variance over ordered nodes in the nondegenerate regime, while reducing to a variance estimator for residual dyad shocks in the degenerate regime.
A related but not identical construction is obtained by forming node-level
scores and applying a standard HAC estimator to the ordered sequence
\begin{align} \sum_{r=1}^{n}\sum_{s=1}^{n} k_L(r-s)\widehat{G}_{r}\widehat{G}_{s}',\qquad\widehat{G}_{r} = \sum_{(i,j)\in\mathcal{D}_{n}:r\in\{i,j\}} \widehat{s}_{ij}, \qquad r=1,\ldots,n. \label{eq:HAC node}\end{align}
It can be expanded as
\[
\sum_{(i,j)\in\mathcal{D}_n}\sum_{(p,q)\in\mathcal{D}_n}
\{k_L(|i-p|)+k_L(|i-q|)+k_L(|j-p|)+k_L(|j-q|)\}
\widehat s_{ij}\widehat s_{pq}'.
\]
This expression is not algebraically identical to \eqref{eq:dn_dyadic_meat}. The estimator in \eqref{eq:HAC node} assigns a separate kernel weight to each close endpoint pairing, whereas \eqref{eq:dn_dyadic_meat} assigns a single kernel weight according to the closest endpoint distance. Hence, the two weighting schemes differ for dyad pairs with more than one close endpoint pairing. Another related distinction is that the node-level HAC representation in \eqref{eq:HAC node} counts the same dyad through both of its endpoints. In particular, for the self-pair \((i,j)=(p,q)\), the two zero-distance endpoint pairings \((i,p)\) and \((j,q)\) both contribute, producing a double-counting term. Therefore, the dyad-pair representation in \eqref{eq:dn_dyadic_meat} is the preferred definition.\footnote{Under the bandwidth conditions imposed below, and after properly accounting for the corresponding double-counting term, this difference is asymptotically negligible under the normalization used for the dyadic meat in the nondegenerate node-dependence case.
}
The bandwidth $L$ controls how far the estimator looks along the ordered
node index. A larger $L$ includes more covariance terms and is appropriate
when dependence between nearby nodes is stronger or more persistent. A
smaller $L$ reduces variability when dependence decays quickly. In the implementation, $L$ is selected from the node-score process
$\{\widehat{G}_{r}\}_{r=1}^{n}$. The detailed process is available in Appendix \ref{app:implementation}.
\subsection{Row-column moving-block jackknife}
\label{subsec:jk}
We next define the jackknife analog of the DN-Dyadic CRVE. The key idea is
to delete a moving block of nodes and remove all dyads touching that block.
This is a row-column deletion: deleting node block $B_\ell$ removes both the
rows and the columns associated with those nodes in the dyadic array.
For $\ell=1,\ldots,n-L+1$, define the overlapping
node block\footnote{The JK-DN-Dyadic CRVE uses the same bandwidth \(L\) as the DN-Dyadic CRVE, which ensures consistency across the two implementations. One could instead recompute the bandwidth separately for each jackknife-deleted sample, after removing all dyads that touch the deleted block of nodes. We find that this alternative implementation produces no significant change in the simulation results.}
\begin{equation}
B_\ell=\{\ell,\ell+1,\ldots,\ell+L-1\}.
\label{eq:block}
\end{equation}
The set of dyads touching $B_\ell$ is
\begin{equation}
\mathcal{A}_\ell
=
\{(i,j)\in\mathcal{D}_n:
i\in B_\ell \text{ or } j\in B_\ell\}.
\label{eq:touch_block}
\end{equation}
The delete-block dyadic sample is
\(
\mathcal{D}_{n,-\ell}
=
\mathcal{D}_{n}\setminus\mathcal{A}_\ell.
\)
The corresponding delete-block estimator is
\begin{equation}
\widetilde{\beta}_{(-\ell)}
=
\left(
\sum_{(i,j)\in\mathcal{D}_{n,-\ell}}x_{ij}x_{ij}'
\right)^{+}
\sum_{(i,j)\in\mathcal{D}_{n,-\ell}}x_{ij}y_{ij},
\label{eq:beta_minus_block}
\end{equation}
where $A^{+}$ denotes the Moore-Penrose inverse. In regular cases,
$A^{+}$ equals the usual inverse; it is used here only to make the definition well-defined when a delete-block design matrix is nearly singular in finite samples.
The uncorrected row-column moving-block jackknife variance estimator is
\begin{equation}
\widehat{V}^{\mathrm{JK}}_0
=
\frac{1}{L}
\sum_{\ell=1}^{n-L+1}
\bigl(\widetilde{\beta}_{(-\ell)}-\widehat{\beta}\bigr)
\bigl(\widetilde{\beta}_{(-\ell)}-\widehat{\beta}\bigr)'.
\label{eq:jk_block_uncorrected}
\end{equation}
The normalization $1/L$ is the moving-block jackknife normalization. When
$L=1$, the estimator deletes one node at a time and removes all dyads
involving that node. When $L>1$, it deletes a local block of ordered nodes
and removes all dyads attached to that block.
Because each dyadic observation is attached to two endpoint nodes, the row-column jackknife contains a double-counting component. We therefore use the corrected JK-DN-Dyadic CRVE
\begin{equation}
\widehat{V}^{\mathrm{JK}}_{\mathrm{DN}}
=
\widehat{V}^{\mathrm{JK}}_0
-
(X'X)^{-1}
\left(
\sum_{(i,j)\in\mathcal{D}_{n}}
\widehat{s}_{ij}\widehat{s}_{ij}'
\right)
(X'X)^{-1}.
\label{eq:jk_block}
\end{equation}
The correction subtracts the White component computed from the
full-sample residual scores. We do not recompute this double-counting component
inside each jackknife deletion. This keeps the correction simple and stable
and improves finite-sample behavior.
\section{Asymptotic validity}
\label{sec:theory}
This section states the asymptotic validity of the DN-Dyadic CRVE and
the JK-DN-Dyadic CRVE. Throughout, for a matrix $A$, we write $A>0$ to denote that the matrix $A$ is positive definite. Let
$ Q_n
=
\frac{1}{M_n}
\sum_{(i,j)\in\mathcal D_n}x_{ij}x_{ij}',
\
Q=\lim_{n\to\infty}Q_n .$
The OLS expansion is
\begin{equation}
\widehat\beta-\beta
=
Q_n^{-1}
\frac{1}{M_n}
\sum_{(i,j)\in\mathcal D_n}s_{ij}.
\label{eq:ols_expansion_raw}
\end{equation}
We use the projection notation from Section~\ref{sec:model}. In the
nondegenerate case, the leading term is the first-order node projection $\{\gamma_i\}$. In the degenerate case considered below, the first-order node projection disappears, and the leading term is the
dyad-level residual component.
\begin{assumption}[Node projection and moments]
\label{ass:projection}
For some $\delta>0$ and $\lambda>1$,
\(
E(x_{ij}u_{ij})=0,\
E(x_{ij}x_{ij}')>0,\
E\|x_{ij}\|^{8(\lambda+\delta)}<\infty,
\
E|u_{ij}|^{8(\lambda+\delta)}<\infty .
\)
\end{assumption}
\begin{assumption}[Weak dependence]
\label{ass:mixing}
The sequence $\{Z_i\}$ is strictly stationary with
mixing coefficients $\beta(h)$ satisfying, for $\lambda$ defined in Assumption \ref{ass:projection},
\(
\beta(h)=O(h^{-\mathfrak d})\) for some \(
\mathfrak d>\frac{2\lambda}{\lambda-1}.
\)
\end{assumption}
Assumptions \ref{ass:projection} impose standard moment conditions; see, e.g., \citet{chiang2023standard} and \citet{chen2023fixed}. Assumption \ref{ass:mixing} imposes a \(\beta\)-mixing condition in order to invoke
the degenerate U-statistic result of \citet{yoshihara1976limiting}. This condition
can be weakened to \(\alpha\)-mixing at the cost of imposing additional smoothness
on $\xi_{ij}$, such as a Lipschitz-type continuity condition; see, for example,
\citet{jochmans2026two}. Define the long-run variance of the first-order node projection by
\begin{equation}
\Omega_\gamma
=E(\gamma_1\gamma_{1}')+
\sum_{h=1}^{\infty}
E(\gamma_1\gamma_{1+h}'+\gamma_{1+h}\gamma_1').
\label{eq:longrun}
\end{equation}
For the degenerate case, define
\(
v_h
=
E\!\left[
E\left(\zeta_{1,1+h}\zeta_{1,1+h}'\mid Z_1,Z_{1+h}\right)
\right],
\)
and
\begin{equation}
\Omega_\zeta
=
\lim_{n\to\infty}
\frac{4}{(n-1)^2}
\sum_{h=1}^{n-1}(n-h)v_h .
\label{eq:omega_zeta}
\end{equation}
\begin{assumption}[Variance]
\label{ass:variance}
One of the following two cases holds.
(i)
$\Omega_\gamma>0$; or
(ii) $\operatorname{Var}(\gamma_{i})=0$,
$\operatorname{Var}(\xi_{ij})=0$, and $\Omega_\zeta>0$.
\end{assumption}
Assumption~\ref{ass:variance} distinguishes two cases. In the
nondegenerate case, the first-order node projection $\gamma_i$ contributes
to the leading sampling variation. In the degenerate case imposed in Assumption \ref{ass:variance}(ii), the
first-order projection $\gamma_i$ and the
non-Gaussian component $\xi_{ij}$ are negligible, but the remaining dyad-level component based on $\zeta_{ij}$
has a nonzero limiting variance. The assumption therefore ensures that the
limiting distribution is Gaussian under the relevant normalization.\footnote{When the second-order component \(\xi_{ij}\) is not negligible, the limiting
distribution is generally non-Gaussian. For two-way clustering, max-type
statistics can deliver conservative inference because the two clustering
dimensions are distinct, so one can condition on one dimension and use the
other for one-way normalization; see \citet{mackinnon2024jackknife} and
\citet{davezies2025analytic}. This logic does not directly extend to dyadic
data, where both indices refer to the same node population and rows and
columns cannot be separated into two independent clustering directions.}
\begin{theorem}[Limit distribution]
\label{thm:clt}
Suppose Assumptions~\ref{ass:representation}-\ref{ass:variance} hold.
If Assumption~\ref{ass:variance}(i) holds, then
\begin{equation}
\sqrt n(\widehat\beta-\beta)
\Rightarrow
N(0,V_\gamma),
\qquad
V_\gamma
=
4Q^{-1}\Omega_\gamma Q^{-1}.
\label{eq:clt_beta_gamma}
\end{equation}
If Assumption~\ref{ass:variance}(ii) holds, then
\begin{equation}
n(\widehat\beta-\beta)
\Rightarrow
N(0,V_\zeta),
\qquad
V_\zeta
=
Q^{-1}\Omega_\zeta Q^{-1}.
\label{eq:clt_beta_zeta}
\end{equation}
\end{theorem}
The first result is the ordered-node analog of the standard
nondegenerate dyadic limit theory. The factor four in
\eqref{eq:clt_beta_gamma} comes from the fact that each node contributes
to approximately $n-1$ dyads. The second result covers the degenerate
case in which the first-order node projection is absent. In that case,
the rate becomes $n$ because the leading variation is generated by the
dyad-level residual component.
\begin{assumption}[Bandwidth]
\label{ass:bandwidth}
As $n\to\infty$, the bandwidth $L=L_n$ satisfies
$ L\to\infty$ and $
{L^2}/{n}=o(1)$.
\end{assumption}
\begin{theorem}
\label{thm:hac_jk}
Suppose Assumptions~\ref{ass:representation}-\ref{ass:bandwidth} hold.
If Assumption~\ref{ass:variance}(i) holds, then
\begin{equation}
n\widehat V_{\mathrm{DN}}
\to^P
V_\gamma,
\qquad
n\widehat V^{\mathrm{JK}}_{\mathrm{DN}}
\to^P
V_\gamma .
\label{eq:var_consistency_gamma}
\end{equation}
If Assumption~\ref{ass:variance}(ii) holds, then
\begin{equation}
n^2\widehat V_{\mathrm{DN}}
\to^P
V_\zeta,
\qquad
n^2\widehat V^{\mathrm{JK}}_{\mathrm{DN}}
\to^P
V_\zeta .
\label{eq:var_consistency_zeta}
\end{equation}
Consequently, for every fixed nonzero vector $a$,
\begin{equation}
\frac{a'(\widehat\beta-\beta)}
{\sqrt{a'\widehat V_{\mathrm{DN}}a}}
\Rightarrow
N(0,1),
\qquad
\frac{a'(\widehat\beta-\beta)}
{\sqrt{a'\widehat V^{\mathrm{JK}}_{\mathrm{DN}}a}}
\Rightarrow
N(0,1).
\label{eq:tstat_validity}
\end{equation}
\end{theorem}
Theorem~\ref{thm:hac_jk} shows that both proposed variance estimators
adapt to the relevant source of first-order variation. In the
nondegenerate case, both estimators consistently estimate the variance
of the $\sqrt n$ limit. In the degenerate case, both estimators
consistently estimate the variance of the $n$ limit. Therefore, the
studentized statistics in \eqref{eq:tstat_validity} are asymptotically
standard normal in both cases.
\section{Simulation evidence}
\label{sec:simulation}
This section studies the finite-sample performance of the proposed
dependent-node dyadic inference methods. The simulation uses the linear
dyadic regression model
\begin{equation}
y_{ij}=x_{ij}'\beta+u_{ij},
\qquad 1\le i<j\le n,
\label{eq:sim_model}
\end{equation}
where $\beta=(1,\ldots,1)'\in\mathbb{R}^{K}$. The null hypothesis concerns
the last component of $\beta$, and all tests are conducted at the nominal
$5\%$ significance level.
The data-generating process is designed to generate two forms of dependence.
First, two dyads that share a node are dependent through common latent node components. Second, because the latent ordered node components are dependent, two dyads that do not share a node may also be dependent when their endpoint nodes are close. Specifically, for each node $i$, let
$A_i^x\in\mathbb{R}^{K}$ and $A_i^u\in\mathbb{R}$ denote latent node shocks
generated by the stationary AR(1) processes
\begin{align}
A_i^x &= \rho A_{i-1}^x+\sqrt{1-\rho^2}\,\eta_i^x,
\qquad \eta_i^x\sim N(0,I_K), \label{eq:sim_Ax}\\
A_i^u &= \rho A_{i-1}^u+\sqrt{1-\rho^2}\,\eta_i^u,
\qquad \eta_i^u\sim N(0,1), \label{eq:sim_Au}
\end{align}
with innovations independent across $i$ and independent of all dyad-specific
shocks. The parameter $\rho\in[0,1)$ controls the strength of ordered-node
dependence. When $\rho=0$, the latent node shocks are independent over the
node index. When $\rho$ is large, nearby nodes are strongly dependent.
For each dyad $(i,j)$, the regressors and disturbance are generated as
\begin{align}
x_{ij} &= \omega(A_i^x+A_j^x)+e_{ij}^x, \label{eq:sim_x}\\
v_{ij} &= \omega(A_i^u+A_j^u)+e_{ij}^u, \label{eq:sim_v}\\
u_{ij} &= \{1+\gamma |x_{ij,K}|\}v_{ij}, \label{eq:sim_u}\\
y_{ij} &= x_{ij}'\beta+u_{ij}, \label{eq:sim_y}
\end{align}
where $e_{ij}^x\sim N(0,I_K)$ and $e_{ij}^u\sim N(0,1)$ are independent
dyad-specific shocks. The first component of $x_{ij}$ is then set equal to
one so that the regression includes an intercept. The parameter $\omega$
controls the strength of dyadic dependence generated by the latent node
components. When $\omega=0$, the common node components do not enter the DGP, and the dyadic dependence is weak. As $\omega$ increases, shared-node and ordered-node dependence become stronger. The parameter $\gamma$ controls the degree of conditional heteroskedasticity through the last regressor $x_{ij,K}$.\footnote{We also study different forms of heteroskedasticity through all regressors $\{x_{ij,k}\}_k$, and the result demonstrates a similar pattern.} The baseline design sets
\[
n=50,\qquad K=10,\qquad \omega=1,\qquad \gamma=0.5,
\]
and uses $5{,}000$ Monte Carlo replications.
\begin{figure}[t]
\centering
\includegraphics[width=0.55\linewidth]{vary_rho.png}
\caption{Rejection frequencies for dyadic inference methods, varying ordered-node dependence $\rho$.}
\label{fig:main}
\end{figure}
We compare five inference procedures:
\begin{enumerate}
\item \textbf{White}: the heteroskedasticity-robust estimator, which
treats all dyads as independent;
\item \textbf{TW}: the conventional two-way cluster-robust estimator
based on the two dyadic indices;
\item \textbf{Dyadic}: the conventional dyadic CRVE, which accounts for
arbitrary dependence between dyads that share a node, but does not
account for ordered-node dependence between distinct nodes;
\item \textbf{DN-Dyadic}: the dependent-node dyadic CRVE, which accounts
for shared-node dependence and ordered-node dependence;
\item \textbf{JK-DN-Dyadic}: the proposed row-column moving-block
jackknife, which deletes adjacent blocks of ordered nodes and removes
all dyads touching the deleted block.
\end{enumerate}
When the same bandwidth choice is used, the HAC and bootstrap procedures proposed by \citet{jochmans2026two} perform similarly to, and slightly better than, DN-Dyadic, but remain less accurate than JK-DN-Dyadic in the presence of node dependence. The difference arises because those procedures do not implement the double-counting correction. As a result, the estimated variance tends to be slightly larger, leading to somewhat more conservative tests. The trade-off is that these procedures become overly conservative when node dependence is absent.
This distinction is well known in the comparison between CRVEs without double-counting correction and CRVEs with double-counting correction in conventional two-way clustering; see \citet{cameron2011robust} and \citet{davezies2021empirical}. See also \citet{mackinnon2021wild} for theoretical results covering both approaches, and \citet{chiang2023standard} and \citet{chen2023fixed} for analogous methods in two-way clustering settings with a time dimension. For clarity of exposition, we report the additional simulation results in Appendix \ref{app:implementation}, including the naive iid homoskedastic variance estimator, one-way clustering CRVE, the \citet{jochmans2026two} method, and the jackknife procedure without double-counting correction.
Figure \ref{fig:main} varies the ordered-node dependence parameter $\rho$,
holding the other parameters at their baseline values. When $\rho$ is small,
ordered-node dependence is weak, and the conventional dyadic CRVE performs
reasonably well. The DN-Dyadic estimator is slightly more conservative in
this region, reflecting the finite-sample cost of allowing for additional
local dependence. As $\rho$ increases, however, all methods exhibit worse performance, and White, TW, and Dyadic exhibit more size distortion. This pattern is consistent with
their dependence restrictions: White ignores dependence, TW captures only
part of the dyadic dependence, and the conventional dyadic CRVE captures
shared-node dependence but not ordered-node dependence. The DN-Dyadic estimator improves size control over Dyadic, while
the JK-DN-Dyadic estimator is relatively robust to the varying level of ordered-node dependence compared to all other methods.
\begin{figure}[t!]
\centering
\begin{subfigure}{0.48\textwidth}
\caption{Varying dyadic dependence strength $\omega$}
\label{fig:rho70_omega}
\centering
\includegraphics[width=\textwidth]{vary_omega_50.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of nodes $n$}
\label{fig:rho70_n}
\centering
\includegraphics[width=\textwidth]{vary_n_50.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of regressors $K$}
\label{fig:rho70_K}
\centering
\includegraphics[width=\textwidth]{vary_K_50.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying heteroskedasticity $\gamma$}
\label{fig:rho70_gamma}
\centering
\includegraphics[width=\textwidth]{vary_gamma_50.png}
\end{subfigure}
\caption{Rejection frequencies for dyadic inference methods under moderate ordered-node dependence, $\rho=0.50$. The nominal significance level is $5\%$.}
\label{fig:main_rho50}
\end{figure}
Figure \ref{fig:main_rho50} fixes the ordered-node dependence parameter at the moderate level
$\rho=0.50$ and varies one design parameter at a time. Panel (a) varies $\omega$, which controls the strength of the latent node
component and hence the strength of dyadic dependence. When $\omega$ is close
to zero, the common node component is weak and the dyads are nearly
independent apart from the idiosyncratic shocks. In this case, the White method
is close to the nominal level, while DN-Dyadic can be conservative because it allows for additional local dependence. As $\omega$ increases, shared-node dependence becomes
stronger, and the methods that do not fully account for the dyadic dependence
begin to over-reject. White exhibits the largest size distortion because it ignores the dependence structure altogether. The two-way and conventional dyadic estimators improve upon White, reflecting their ability to account for part or all of the shared-node dependence. DN-Dyadic further improves slightly upon the conventional dyadic estimator.
The proposed jackknife estimator remains closest
to the nominal level, indicating that the row-column block deletion provides
additional finite-sample robustness.
Panel (b) varies the number of nodes $n$. White remains substantially
oversized as $n$ increases, whereas the other methods improve, reflecting
that they at least partially account for the dependence structure. The
performance of the Dyadic and DN-Dyadic estimators improves with $n$, but
they remain somewhat oversized in finite samples. Interestingly, the JK-DN-Dyadic estimator is already close to the nominal level when $n=10$. Although its rejection frequency increases slightly as the sample size becomes larger, it remains much more stable than the competing procedures. This suggests that the row-column block deletion delivers useful robustness even in very small samples.
Panels (c) and (d) vary the number of regressors $K$ and the
heteroskedasticity parameter $\gamma$, respectively. The rejection
frequencies are relatively stable across these variations, suggesting that
the main source of size distortion in this design is the dependence structure
rather than the number of regressors or the degree of heteroskedasticity.
Across all panels, the qualitative ranking of the methods is unchanged: White performs worst, the two-way and conventional Dyadic estimators improve upon White, DN-Dyadic performs slightly better than Dyadic under moderate ordered-node dependence, and JK-DN-Dyadic delivers the most reliable size control.
\begin{figure}[t]
\centering
\begin{subfigure}{0.48\textwidth}
\caption{Varying dyadic dependence strength $\omega$}
\label{fig:rho30_omega}
\centering
\includegraphics[width=\textwidth]{vary_omega.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of nodes $n$}
\label{fig:rho30_n}
\centering
\includegraphics[width=\textwidth]{vary_n.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of regressors $K$}
\label{fig:rho30_K}
\centering
\includegraphics[width=\textwidth]{vary_K.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying heteroskedasticity $\gamma$}
\label{fig:rho30_gamma}
\centering
\includegraphics[width=\textwidth]{vary_gamma.png}
\end{subfigure}
\caption{Rejection frequencies for dyadic inference methods under weak ordered-node dependence, $\rho=0.30$. The nominal significance level is $5\%$.}
\label{fig:main_rho30}
\end{figure}
\begin{figure}[h!]
\centering
\begin{subfigure}{0.48\textwidth}
\caption{Varying dyadic dependence strength $\omega$}
\label{fig:rho70_omega}
\centering
\includegraphics[width=\textwidth]{vary_omega_70.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of nodes $n$}
\label{fig:rho70_n}
\centering
\includegraphics[width=\textwidth]{vary_n_70.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying the number of regressors $K$}
\label{fig:rho70_K}
\centering
\includegraphics[width=\textwidth]{vary_K_70.png}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\caption{Varying heteroskedasticity $\gamma$}
\label{fig:rho70_gamma}
\centering
\includegraphics[width=\textwidth]{vary_gamma_70.png}
\end{subfigure}
\caption{Rejection frequencies for dyadic inference methods under strong ordered-node dependence, $\rho=0.70$. The nominal significance level is $5\%$.}
\label{fig:main_rho70}
\end{figure}
Figures \ref{fig:main_rho30} and \ref{fig:main_rho70} repeat the same experiments under weak and strong ordered-node dependence, with \(\rho=0.30\) and \(\rho=0.70\), respectively. The qualitative patterns are similar to that in Figure \ref{fig:main_rho50}, but the inference problem becomes more difficult at \(\rho=0.70\). White, TW, and Dyadic exhibit more severe over-rejection. The DN-Dyadic estimator substantially improves upon Dyadic, especially as \(n\) increases, because the ordered-node dependence is stronger and accumulates over a larger number of nodes. The JK-DN-Dyadic estimator delivers the most robust size control among the methods considered, although it can still over-reject when the dependence is strong. Furthermore, the results in Figure \ref{fig:main_bandwidth} demonstrate that the selected bandwidth is relatively robust and adapts well to different settings.
Overall, the simulation evidence supports the main message of the paper.
When the node index is ordered, and nearby nodes are dependent, conventional two-way clustering or
dyadic inference can be unreliable because dyads with no common node may
still be correlated. Accounting for ordered-node dependence improves size
control, and the row-column moving-block jackknife provides the most robust
finite-sample performance. We therefore recommend the JK-DN-Dyadic estimator as the default procedure for applications in which ordered-node dependence may be empirically meaningful.
\section{Empirical Illustration: Impact of Free Trade Agreements on Trade}
\label{sec:empirical_trade}
We use the proposed inference procedures to revisit a central question in international trade: do free trade agreements (FTAs) significantly increase bilateral trade? This question is both empirically important and policy relevant. FTAs are among the most widely used policy instruments for reducing trade barriers, strengthening economic integration, and reshaping global trade patterns. At the same time, their empirical effects remain actively debated, because countries do not enter FTAs randomly and because bilateral trade flows are subject to rich cross-country dependence; see, for example, \citet{baier2007free}, \citet{magee2008new}, and \citet{egger2008interdependent}. Gravity regressions provide the standard empirical framework for studying this question because they relate bilateral trade flows to trade costs, country-pair characteristics, and trade-policy variables.\footnote{The gravity specification follows the extensive empirical trade literature initiated by \citet{tinbergen1962shaping} and further developed by \citet{anderson1979theoretical}, \citet{anderson2003gravity}, and \citet{silva2006log}. Similar dyadic regression frameworks are widely used to study the determinants of bilateral trade flows and international economic integration.} We estimate the gravity model using the CEPII Gravity Database of \citet{conte2022cepii}. The sample consists of ($n=156$) countries observed from 1996 to 2000. To obtain a cross-sectional dyadic dataset, we average the variables over this period for each country pair.
We order countries by their average GDP per capita. {The ordering is constructed from predetermined average GDP-per-capita measures rather than estimated from the regression residuals.} Countries at similar levels of development may be exposed to similar global demand shocks, financial conditions, supply-chain disruptions, institutional constraints, and trade-policy environments. These common forces may induce dependence not only between dyads sharing a country, but also between dyads whose endpoint countries are close in the economic ordering. For example, trade flows among high-income economies may respond similarly to global financial conditions or supply-chain disturbances, even when the corresponding country pairs do not overlap. {In the application, the data-driven bandwidth selector, implemented as described in Appendix \ref{app:implementation}, chooses \(L=7\), suggesting that the relevant dependence extends beyond exact country overlap.}
The dependent variable is \(y_{ij}=\log(1+\text{Manufacturing Trade}_{ij})\), where \(\text{Manufacturing Trade}_{ij}\) denotes the undirected BACI manufacturing trade flow between countries \(i\) and \(j\). We estimate the following gravity specification:
\[
y_{ij}
=
\alpha_i+\alpha_j
+
\beta_1\text{FTA}_{ij}
+
\beta_2\text{Language}_{ij}
+
\beta_3\log(\text{Distance}_{ij})
+
\beta_4\text{Border}_{ij}
+
\beta_5\text{Sibling}_{ij}
+
u_{ij},
\ i<j.
\]
Here, \(\alpha_i\) and \(\alpha_j\) are country fixed effects. The bilateral controls include a common official language indicator, log distance, a common-border indicator, and an indicator for whether the country pair ever shared the same colonizer.\footnote{Our objective is not to identify a causal effect of FTAs, but rather to illustrate how alternative dyadic
inference procedures affect statistical conclusions in a standard gravity framework.} Our primary parameter of interest is \(\beta_1\), which measures the association between FTA coverage and bilateral manufacturing trade after controlling for country fixed effects and standard gravity covariates.\footnote{Country fixed effects are included to absorb country-level heterogeneity. Although the theoretical results are stated without explicitly modeling fixed effects, the empirical exercise applies the proposed inference procedure to the corresponding fixed-effect transformed estimating equation.
}
\begin{table}[!t]
\centering
\caption{Inference for the FTA coefficient in the manufacturing-trade gravity regression}
\label{tab:empirical_fta}
\begin{tabular}{lccccc}
\hline\hline
Estimate & White & TW & Dyadic & DN-Dyadic & JK-DN-Dyadic \\
\hline
0.1680 & 0.0125 & 0.0582 & 0.0767 & 0.1010 & 0.1198 \\
\hline\hline
\end{tabular}
\begin{flushleft}
\footnotesize
Notes: The table reports the estimated FTA coefficient and the corresponding \(p\)-values under different inference procedures. The dependent variable is \(\log(1+\texttt{manuf\_tradeflow\_baci})\). The regression includes country fixed effects, common language, log distance, common border, sibling-pair status, and the FTA indicator. The node ordering is based on countries' average GDP per capita. The selected bandwidth is \(L=7\), and the sample contains \(n=156\) countries.
\end{flushleft}
\end{table}
Table \ref{tab:empirical_fta} reports the estimated FTA coefficient and the corresponding \(p\)-values. The point estimate is positive, equal to 0.1680, which is consistent with the view that FTAs are associated with higher bilateral manufacturing trade. However, the statistical conclusion depends substantially on how cross-dyad dependence is handled. Under White standard errors, the \(p\)-value is 0.0125, suggesting a statistically significant FTA effect. Once dyadic dependence is taken into account, the evidence becomes weaker: the two-way \(p\)-value increases to 0.0582, and the conventional dyadic \(p\)-value increases to 0.0767. {The increase in estimated uncertainty becomes even more pronounced once ordered-node dependence across economically similar countries is incorporated.} The \(p\)-value rises to 0.1010 under DN-Dyadic and to 0.1198 under JK-DN-Dyadic.
{The progressive increase in \(p\)-values across inference procedures indicates that accounting for richer dependence structures leads to substantially larger estimated standard errors.}
These results highlight the empirical relevance of node dependence in gravity applications. If shared-node or ordered-node dependence across economically similar countries is ignored, the evidence in favor of a statistically significant FTA effect appears stronger. {Overall, once both shared-node and ordered-node dependencies are accounted for, the statistical evidence in favor of a significant FTA effect becomes substantially weaker. Under the proposed JK-DN-Dyadic procedure, we claim that the estimated effect of FTAs on bilateral manufacturing trade flows is not statistically significant at the 10\% level.
}
\section{Conclusion}
\label{sec:conclusion}
This paper studies inference for dyadic regressions when the nodes are ordered
and the latent node shocks are weakly dependent along the node index. In this
setting, conventional dyadic clustering can be insufficient because two dyads
may remain correlated even when they do not share a node, provided that their
endpoint nodes are sufficiently close. The key observation is that the leading
component of the dyadic score behaves like a weakly dependent sequence indexed
by nodes. {Consequently, when such ordered-node dependence is present, valid inference must account not only for shared-node dependence, but also for local dependence along the ordered node index.}
{We propose two variance estimators. The first is a dependent-node dyadic CRVE
that retains covariance terms between dyads with nearby endpoints. The second
is a row-column moving-block jackknife that deletes adjacent blocks of nodes
together with all dyads touching the deleted block. This deletion scheme
preserves both shared-node dependence and ordered-node dependence. Under
standard moment and weak-dependence conditions, we show that both estimators
consistently estimate the asymptotic variance and deliver valid studentized
inference.}
{The Monte Carlo evidence supports the theory and suggests that the proposed row-column moving-block jackknife provides a reliable default procedure for dyadic applications with dependent ordered nodes. The
empirical illustration based on international trade gravity regressions further
shows that accounting jointly for shared-node dependence and ordered-node
dependence can substantially weaken the statistical evidence in favor of free
trade agreement effects on bilateral manufacturing trade flows.}
\clearpage