Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,006 characters · 11 sections · 42 citation commands
Robust Inference for Dyadic Data with Dependent Ordered Nodes
\noindentKeywords: dyadic data, locally dependent, cluster-robust variance estimation, jackknife.\\ JEL Classification: C12, C15, C21, C31.
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., cameron2011robust, thompson2011simple, aronow2015cluster, tabord2019inference, menzel2021bootstrap, davezies2021empirical, and davezies2025analytic. Related work on network formation and sparse network asymptotics includes fafchamps2007formation and 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 jochmans2026two, who studies non-exchangeable dyadic data with dependence that decays over an ordered index distance. The distinction is useful to make explicit. 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. chiang2023standard, chen2023fixed, and 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 hansen2022jackknife, mackinnon2023leverage, mackinnon2024jackknife, and 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) introduces the dyadic regression model and ordered-node dependence. Section (ref) defines the DN-Dyadic CRVE and the JK-DN-Dyadic CRVE. Section (ref) presents the asymptotic validity results. Section (ref) reports the simulation evidence. Section (ref) illustrates the practical relevance of the proposed approach through an empirical application. Section (ref) concludes. Proofs are collected in the Appendix.
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
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
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
Thus, the dependence structure relevant for inference is the dependence structure of the dyadic score array $\{s_{ij}:(i,j)\in\mathcal{D}_{n}\}$.
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, aldous1981representations; hoover1979relations; kallenberg1989representation) representation, but adapted to ordered weakly dependent node variables.
Assumption (ref) 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
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
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
and the dyad-level residual component
Then, for each $(i,j)\in\mathcal D_n$,
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 (ref) over all dyads gives
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 (ref) 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.
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.
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 (ref) extends the usual dyadic-neighborhood relation to the ordered-node setting.
Let
be the Bartlett kernel, where $L$ denotes the bandwidth or block length. The DN-Dyadic CRVE meat is
The corresponding variance estimator for $\widehat{\beta}$ is
Although (ref) 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
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 (ref). The estimator in (ref) assigns a separate kernel weight to each close endpoint pairing, whereas (ref) 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 (ref) 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 (ref) 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).
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.}
The set of dyads touching $B_\ell$ is
The delete-block dyadic sample is \( \mathcal{D}_{n,-\ell} = \mathcal{D}_{n}\setminus\mathcal{A}_\ell. \) The corresponding delete-block estimator is
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
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
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.
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
We use the projection notation from Section (ref). 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.
Assumptions (ref) impose standard moment conditions; see, e.g., chiang2023standard and chen2023fixed. Assumption (ref) imposes a \(\beta\)-mixing condition in order to invoke the degenerate U-statistic result of 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, jochmans2026two. Define the long-run variance of the first-order node projection by
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
Assumption (ref) 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)(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 mackinnon2024jackknife and 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.}
The first result is the ordered-node analog of the standard nondegenerate dyadic limit theory. The factor four in (ref) 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.
Theorem (ref) 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 (ref) are asymptotically standard normal in both cases.
This section studies the finite-sample performance of the proposed dependent-node dyadic inference methods. The simulation uses the linear dyadic regression model
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
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
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.
We compare five inference procedures:
When the same bandwidth choice is used, the HAC and bootstrap procedures proposed by 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 cameron2011robust and davezies2021empirical. See also mackinnon2021wild for theoretical results covering both approaches, and chiang2023standard and 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), including the naive iid homoskedastic variance estimator, one-way clustering CRVE, the jochmans2026two method, and the jackknife procedure without double-counting correction.
Figure (ref) 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.
Figure (ref) 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.
Figures (ref) and (ref) 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), 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) 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.
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, baier2007free, magee2008new, and 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 tinbergen1962shaping and further developed by anderson1979theoretical, anderson2003gravity, and 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 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), 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. }
Table (ref) 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. }
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.}