Dyadic Regression with Sample Selection

Dyadic Regression with Sample Selection

\address{University of Wisconsin-Madison} \email{[email removed]}

Keywords: Dyadic Data, Sample Selection, Fixed Effects, Network Formation, Bias Correction.

Introduction

Dyadic data describe pairwise outcomes, such as trade volume between countries. Numerous applications have analyzed such data using the regression model, referred to as dyadic regression. Examples include gravity equations in trade, migration, and urban economics Helpman2008, Moretti2017, Monte2018, and risk-sharing networks in development economics Fafchamps2007. One of the prominent features of dyadic data is the non-negligible number of zeros in the outcomes of interest, \footnote{Helpman2008 document that there was no trade among roughly 50% of country pairs from 1970 to 1997. In 2017, there was no migration among about 60% of country pairs (the author calculated using the data available from the World Bank (\url{https://www.worldbank.org/en/topic/ migrationremittancesdiasporaissues/brief/migration-remittances-data}).} possibly due to economic mechanisms such as prohibitive fixed costs. This paper deals with panel dyadic data, where zeros are prevalent both across cross-sections and over time.

How should we treat zeros in dyadic regression? In applications, zeros are often discarded due to the log-linear specification Moretti2017. The Poisson pseudo-maximum-likelihood (PPML) estimator is also frequently used to avoid discarding zeros and address issues related to log-linearization SantosSilva2006. These approaches implicitly assume that zeros occur exogenously. Since a zero in a pairwise outcome results from no link between two units, we can associate zeros with the underlying network formation mechanism that determines which pairs appear in a sample. If the network is formed endogenously as a result of an interaction between two agents, the empirical practices mentioned above can be subject to sample selection bias, as in J.Heckman1979.

This paper has two primary objectives. First, we aim to jointly model network formation and the outcome generation on such networks. This joint modeling allows identification of the effects of changes in pair-level or individual-level characteristics, separating them from the effects caused by changes in networks. In contrast, the dyadic regression literature has primarily focused on regression with fixed or exogenous networks. Second, we develop a robust inference method that accounts for the dyadic dependence structure. Pairwise outcomes are likely to be dependent on each other through common shocks to individuals. This dyadic dependence can be especially important in the presence of zeros and the network formation because a few individuals can have significantly more links than others, \footnote{For example, in Moretti2017's migration flow data, star scientists' migration from or to California constituted approximately 14% of the links in the sample on average. This percentage is much higher than the expected 2% when considering all potential links in the sample.} which strengthens the influence of shocks to those individuals on the dyadic dependence. At the same time, it is known that with dyadic data, we can have different asymptotic regimes depending on the nature of those individual-level shocks Menzel2021. To be practitioner-friendly, our inference method needs to consider the dyadic dependence and ensure adaptivity to different resulting asymptotic regimes.

Our setup will be a linear panel dyadic regression model, featuring the network formation process as a sample selection mechanism that generates both zeros and unobservable outcomes. To capture the dyadic dependence structure, we incorporate two types of unobservable individual heterogeneity into the model: time-invariant fixed effects and time-varying random effects, which is a new modeling strategy in the literature. We extend Kyriazidou1997's identification argument, originally designed for individualistic data, to dyadic data, and correspondingly propose a semiparametric, kernel-based estimator that assigns weights to pairs whose selection index remains stable over time. A significant challenge we face when analyzing our estimator is the need to address the dependence structure caused by node-level shocks, which is absent in individualistic data models analyzed in Kyriazidou1997. To control for this type of dependence, we utilize the U-statistic-like structure of our estimator, which gives us a mutually uncorrelated decomposition into the node-level Hájek projection part and the dyad-level projection error part.

We show that our estimator is asymptotically normal with two different convergence rates depending on the nature of errors. If the H\'{a}jek projection is non-degenerate (i.e., each summand has positive variance), our estimator achieves $\sqrt{n}$-asymptotic normality, where $n$ is the number of nodes. In this case, we not only have zero asymptotic bias but also share the same convergence rates as the usual fixed effect estimator and PPML estimator when its leading term is also non-degenerate. The latter point implies that there is no loss in effective sample sizes with our estimator for using a kernel-based local method compared with the usual non-weighted estimator. If the H\'{a}jek projection is degenerate, our estimator achieves $\sqrt{Nh_{n}}$-asymptotic normality, where $N\sim n^{2}$ is the number of dyads and $h_{n}$ is a bandwidth. While the usual fixed effect estimator and the PPML estimator can be non-Gaussian in the limit Menzel2021, our estimator is guaranteed to be asymptotically normal regardless of degeneracy. This result is analogous to Hall1984's central limit theorem for degenerate U-statistics, allowing common statistics of interest, such as confidence intervals, to be constructed in a standard manner. In the degenerate case, our estimator exhibits asymptotic bias, which motivates us to introduce a bias correction.

We propose a variance estimator and bias-corrected confidence intervals that adapt to the degeneracy. Our variance estimator is similar to the one proposed by Graham2019 for nonparametric dyadic density estimation. We show that our estimator is consistent for the asymptotic variances in both non-degenerate and degenerate cases, after being rescaled by $\sqrt{n}$ or $\sqrt{N h_{n}}$, respectively. For the bias correction, we use a consistent estimator for the asymptotic bias in the degenerate case. We show that the correction term is negligible in the non-degenerate case after being rescaled by $\sqrt{n}$. Combining both bias-corrected estimator and variance estimator, we can construct bias-corrected confidence intervals for our estimator. These intervals have asymptotically correct sizes regardless of the (non-)degeneracy of the leading term in our estimator.

We conduct a simple simulation exercise to demonstrate the performance of our estimators compared to the usual fixed effect estimator and PPML estimator, as we vary the fraction of selected dyads from 10% to 90%. Our proposed estimator exhibits better finite sample properties than the other two estimators. Our bias-corrected confidence intervals also outperform the alternatives in coverage probabilities, regardless of degeneracy. This result underscores the importance of bias correction in finite samples, even though the asymptotic bias is zero in the non-degenerate case, which is a new finding in the literature.

We apply our estimator to the regression specification proposed by Moretti2017, which estimates the effects of state tax differences on the internal migration flows within the U.S. Comparing our proposed estimator with Moretti2017's, we find that their conclusion, which suggests that state tax differences have a significant impact on internal migration, may not be robust in the presence of a dyadic dependence structure and sample selection biases.

This paper is closely related to the growing literature on dyadic regression Cameron2014, Tabord-Meehan2019, Bonhomme2020, Zeleneev2020, Graham2020, graham2021minimax,sassi2023. With the exception of Bonhomme2020 and Zeleneev2020, most of these papers do not address non-random sample selection, but instead focus on the consequences of dyadic dependence. Bonhomme2020 primarily studies cases where selection is conditionally random with random effects, and briefly discusses conditionally non-random selection without providing a theoretical analysis. Zeleneev2020 investigates identification and estimation in cross-sectional dyadic regression models with more flexible combinations of node-level fixed effects, including fixed selection effects as a special case. However, this flexibility comes at the cost of more complex inference, which is not covered in their paper. In contrast, our paper focuses on models with an additional time dimension, enabling us to develop asymptotic distribution theory and a practical inference method that adapts to degeneracy.

This paper also contributes to the literature on econometric analysis of models with endogenous network formation. Examples include Moon2021, Auerbach2022, and Jochmans2023. While these papers study social interaction/peer effects type models where outcomes of interest are individualistic, our paper studies the direct consequence of network formation on dyadic outcomes.

Model

Setup

\\ There are $n$ nodes in the data (e.g., states, countries), indexed by $i=1,...,n$. Let $\{(X_{it},Z_{it})_{t=1,...,T})\}_{i=1}^{n}$ be a node-level observation, where $X_{it}\in\mathbb{R}^{q_{x}}$ and $Z_{it}\in\mathbb{R}^{q_{z}}$. For each dyad $ij$ and time $t$, $Y_{ijt}\in\mathbb{R}$ is a main outcome, and we observe a binary variable $d_{ijt}\in\{0,1\}$, which indicates that $Y_{ijt}$ is observable \footnote{Since we focus on a linear model, we can interchange unobservability with zero. Alternatively, we can interpret $Y_{ijt}$ as the logarithm of $\tilde{Y}_{ijt}\geq 0$.} only if $d_{ijt}=1$. We can interpret the adjacency matrix $D_{t}\equiv [d_{ijt}]_{i,j=1,...,n}$ as a network that summarizes the existence of interactions between nodes. In this paper, we restrict our attention to a model with $T=2$ and an undirected graph where $Y_{ijt}=Y_{jit}$, $d_{ijt}=d_{jit}$ for all $i,j,t$. We also rule out self-loops by convention: $Y_{iit}=d_{iit}=0$ for all $i,t$. An extension to $T>2$ and a directed graph is discussed in Section 4.1.

The data is generated according to the following model:

The regressors $W_{ijt}\in\mathbb{R}^{q_{w}}$ and $R_{ijt}\in\mathbb{R}^{q_{r}}$ are constructed from some user-specified symmetric functions $w:\mathbb{R}^{q_{x}}\times\mathbb{R}^{q_{x}}\to\mathbb{R}^{q_{w}}$ and $r:\mathbb{R}^{q_{z}}\times \mathbb{R}^{q_{z}}\to\mathbb{R}^{q_{r}}$ such that $w(x,y)=w(y,x)$ and $r(x',y')=r(y',x')$ for any $x,y\in\mathbb{R}^{q_{x}}$ and $x',y'\in\mathbb{R}^{q_{z}}$. For example, we can specify $w$ to be a pairwise summation $w(x,y)=x+y$. The symmetry in these functions is needed as our graphs are undirected; we can relax this requirement with directed graphs, as discussed in Section 4.1. The node-level fixed effects $A_{i},B_{i}\in\mathbb{R}$ are unobservable, and we allow them to correlate with the regressors, as in the usual fixed effect model. The functions $\psi:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ and $\varphi:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ are unknown symmetric functions that capture the interaction between two nodes through their fixed effects.

We specify the structure of errors $\epsilon_{ijt},\eta_{ijt}$ as follows: For $1\leq i<j\leq n$,

where $U_{i}\equiv (U_{i1},U_{i2})$ and $U_{ij}\equiv (U_{ij1},U_{ij2})$ are node-level and dyad-level random vectors, respectively, and $\tau$ is an unknown multivariate function.\footnote{Here, we need not specify the dimensions of those vectors and the function since the following results do not depend on them as long as those dimensions are fixed.}

Let $\xi_{i}\equiv (X_{i1},X_{i2},Z_{i1},Z_{i2},A_{i},B_{i})$ be a vector that contains observed and unobserved information in the two periods with respect to node $i$. We impose the following distributional assumption:

Part (1) imposes homogeneity on the node-level data-generating process. Parts (2) and (3) are new to the literature on dyadic regression with fixed effects. While the previous literature assumes conditional independence of dyadic-level errors Graham2017,Zeleneev2020,Candelaria2020, our error structure ((ref)) allows for the conditional dependence between errors with a common node (e.g., $\epsilon_{ij1}$ and $\epsilon_{ik1}$) through $U_{i}$, but also includes conditional independence as a special case where node-level random vectors $U_{it},U_{jt}$ are degenerate given $\{\xi_{i}\}_{i=1}^{n}$. Part (4) is the standard assumption in the literature and excludes "externalities," where dyad $ij$ can be affected by nodes other than $i$ or $j$. Part (5) ensures the conditional exchangeability of $(\epsilon_{ijt},\eta_{ijt})_{t=1,2}$ across dyads.

Identification

\\ The following two assumptions are crucial for the identification of $\beta$:

Assumption (ref) excludes cases where, for example, the conditional variance of $\epsilon_{ijt}$ depends only on period $t$'s information: $Var(\epsilon_{ijt}|\xi_{i},\xi_{j})=\sigma^{2}\times W_{ijt}'\beta$. However, it allows time invariant heteroskedasticity such as $Var(\epsilon_{ijt}|\xi_{i},\xi_{j})=\sigma^{2}(W_{ij1}+W_{ij2})'\beta\times A_{i}\times A_{j}$. From ((ref)), this assumption is implied by the conditional exchangeability of $U_{it}$ and $U_{ijt}$ with respect to time and the symmetry of function $\tau$ in the sense that $\tau(u_{1},u_{2},v_{1},v_{2},w_{1},w_{2})=\tau(u_{2},u_{1},v_{2},v_{1},w_{2},w_{1})$ for any $u_{1},u_{2},v_{1},v_{2},w_{1},w_{2}$. Assumption (ref) excludes cases where $W_{ijt}$ is exactly the same as $R_{ijt}$ and implies that some variables in $R_{ijt}$ must be excluded from $W_{ijt}$. Since

this assumption also implies that the networks $D_{1},D_{2}$ are locally dense across time in the sense that $Pr(d_{ij1}d_{ij2}=1|\Delta R_{ij}'\gamma=0)>0$.

Our identification argument is summarized in the following two steps, similarly to Kyriazidou1997. First, take the time-difference on observed outcomes (dyads with $d_{ij1}=d_{ij2}=1$) to eliminate the fixed effects:

If we take expectation of both sides conditionally on $d_{ij1}=d_{ij2}=1$ and $\xi_{i},\xi_{j}$,

Note that, in general, the sample selection effect is not $0$.

Second, we seek to find conditions to eliminate the selection effect. Assumption (ref) is equivalent to

where $F$ is the conditional distribution of the errors given $\xi_{i},\xi_{j}$. Then, for dyad $ij$ with $\Delta R_{ij}'\gamma=R_{ij1}'\gamma-R_{ij2}'\gamma=0$,

Hence, the conditional expectation of $\Delta Y_{ij}$ given $d_{ij1}d_{ij2}=1$, $\xi_{i},\xi_{j}$, and $\Delta R_{ij}'\gamma=0$ is

Multiplying the both sides by $\Delta W_{ij}$ and aggregating $\xi_{i},\xi_{j}$, we get

Then, under Assumption (ref), $\beta$ is uniquely written as

Estimation

\\ Estimation is done in two steps. In the first step, we estimate $\gamma$ with a consistent estimator $\hat{\gamma}_{n}$, and in the second step we estimate $\beta$ with $\hat{\beta}_{n}$, a sample analogue of the identified $\beta$ with $\gamma$ replaced by $\hat{\gamma}_{n}$.

In the following, we focus on the second step. The sample-analogue of ((ref)) is given by

where $\sum_{i<j}=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}$, $K_{h_{n}}(v)=h_{n}^{-1}K(v/h_{n})$ is a kernel, and $h_{n} $ is a bandwidth. The weight function is used to smooth the condition $\Delta R_{ij}'\gamma=0$ and puts larger weight on observations with small $\Delta R_{ij}'\hat{\gamma}_{n}$.

To evaluate $\hat{\beta}_{n}$ in terms of $\beta$, rewrite the time-differenced model as

where

Note that $E[\nu_{ij}|d_{ij1}d_{ij2}=1,\xi_{i},\xi_{j}]=0$ by construction. Define

Substituting $\Delta Y_{ij}$ into $\hat{\beta}_{n}$ yields

The terms $\hat{S}_{WW}^{-1}\hat{S}_{W\lambda}$ and $\hat{S}_{WW}^{-1}\hat{S}_{W\nu}$ can be understood as the selection bias term and the stochastic error term of the estimator, respectively.

Asymptotic Analysis

Regularity Conditions

\\ For ease of notation, we write the following conditions in terms of dyads $12$ and $13$, which entails no loss of generality under the undirected graph and Assumption (ref).

Let $f_{R\gamma,2}$ be the joint density of $\Delta R_{12}'\gamma$ and $\Delta R_{13}'\gamma$ when it exists and $f_{R\gamma,2|\xi_{2},U_{2},\xi_{3},U_{3}}$ be the conditional density given $\xi_{2},U_{2},\xi_{3},U_{3}$. Let $f_{R\gamma}$ be the marginal density and $f_{R\gamma|\xi_{1},U_{1}}$ be the conditional density given $\xi_{1},U_{1}$.

Part (1) is a smoothness assumption on the density as in the nonparametric regression literature. Part (3) ensures that we observe $\Delta R_{12}'\gamma$ around $0$, which is crucial for identification. Parts (2) and (4) essentially requires well-behaved $r(\cdot,\cdot)$ in ((ref)).

Define $(w_{1},w_{2})\mapsto \Lambda(w_{1},w_{2},\xi_{1},\xi_{2})$ as

with $t,s=1,2,t\neq s$. This $\Lambda$ is the sample selection effect caused by the correlation between errors $\epsilon_{12t},\epsilon_{12s}$ and $\eta_{12t},\eta_{12s}$. Note that the function $\Lambda$ does not depend on time $t$ or $s$ because of Assumption (ref).

This assumption is essential for controlling the sample selection effect and characterizing the asymptotic bias in some cases. An implication of this assumption is that for some $\Lambda_{12}\equiv \tilde{\Lambda}(w_{1},w_{2},\xi_{1},\xi_{2})$,

by the multivariate mean-value theorem. Note that the function $\Lambda$ does not depend on time $t$ or $s$ because of Assumption (ref). This assumption is strong because the difference in $\Lambda$ must be exactly linear in the first and second elements. If we focus on the degenerate case discussed below, since the asymptotic bias is $0$ in that case, we can relax the differentiability to Lipschitz-like continuity on $\Lambda$: $|\Lambda(w_{1},w_{2},\xi_{1},\xi_{2})-\Lambda(w_{2},w_{1},\xi_{1},\xi_{2})|\leq |\Lambda_{12}|\times |w_{1}-w_{2}|$.

Let $\|\cdot\|$ denote a Euclidian norm of vectors.

Part (1) assumes the existence of conditional moments for the relevant variables. The conditioning on $\xi_{1}$ and $U_{1}$ is needed for controlling the dyadic dependence structure. Part (2) is crucial for obtaining the convergence results used below, and the positive definiteness is needed for ensuring the non-degeneracy of our estimator in the limit. Part (3) is used for characterizing the asymptotic bias provided below. Part (4) is essential for the negligibility of the approximation error of our variance estimator.

Additionally to Assumption (ref), which restricts the moments locally around $(0,0)$ or $(0)$, we use the existence of these unconditional moments when bounding error terms coming from the usage of $\hat{\gamma}_{n}$.

For example, a fourth-order biweight kernel $K(x)=106/64(1-3u^{2})(1-x^{2})^{2}\boldsymbol{1}\{|x|< 1\}$ satisfies this assumption with $\kappa=1$ and $k=3$.

This assumption is standard in the nonparametric regression literature. We impose further conditions on $\{h_{n}\}$ in each statement below.

This assumption requires the first-step estimator to be consistent and converge faster than our estimator. For example, if $\eta_{ijt}\sim Logistic(0,1)$ independently across $ij$ and $t$, we can show that Chamberlain1980's conditional logit estimator satisfies $\hat{\gamma}_{n}-\gamma=O_{p}(1/\sqrt{N})$ so that $\sqrt{Nh_{n}}(\hat{\gamma}_{n}-\gamma)=O_{p}(\sqrt{h_{n}})=o_{p}(1)$. In Section 3.6, we discuss the availability of alternative estimators for $\gamma$. We leave the case where $\hat{\gamma}_{n}$ converges slower than required in this assumption for future research.

Asymptotic Normality

\\ Define the following components that will appear in the asymptotic bias and variance expression:

We have the following result:

Parts (1) and (2) of Theorem (ref) show that our estimator is asymptotically normal, with different convergence rates depending on $\Sigma_{W\nu,1}$. Part (1) differs from Kyriazidou1997 in that the convergence rate is parametric and based on the number of nodes $n$, rather than the number of dyads $N$. When $c_{W}'\Sigma_{W\nu,1}c_{W}>0$ in part (1), the covariance between summands sharing a common node (e.g., dyads $ij$ and $ik$) does not vanish asymptotically, reducing the effective sample size to $n$. The leading term is an average of conditional means given $\xi_{i}$ and $U_{i}$, which averages out and eliminates $h_{n}$ from the convergence rate. This $\sqrt{n}$-asymptotic normality matches results in the dyadic nonparametric density estimation literature Graham2019. When $\Sigma_{W\nu,1}=0$, as in part (2), our result aligns with Kyriazidou1997, with nonparametric convergence rates based on the number of dyads. This corresponds to the degenerate case for dyadic dependence, as discussed in the literature Graham2019,cattaneo2024. Part (3) of Theorem (ref) shows that, with suitable normalization, our estimator converges to the asymptotic bias term regardless of degeneracy. This property is used to construct the bias-corrected estimator in the following section.

We can compare our estimator with the usual fixed effect estimator:

which is biased because of the selection effect $\lambda_{ij}$. First, in the case of non-degeneracy, our estimator and the re-centered (infeasible) fixed effect estimator share the same convergence rates of $\sqrt{n}$ Daveziez2021. This implies that there is no reduction in the effective sample size for using our kernel-based local estimator, which amends the need for fairly large samples as discussed in Kyriazidou1997. Second, in the case of degeneracy, the fixed effect estimator applied to our model can exhibit a non-Gaussian distribution in the limit Menzel2021, while our estimator is asymptotically normal regardless of the degeneracy. This guaranteed asymptotic normality is analogous to Hall1984's central limit theorem for degenerate U-statistics, and thus the common statistics of interest, such as confidence intervals, can be constructed in a standard manner.

If we interpret the structural equation (ref) as the log-linearized version of the canonical gravity model SantosSilva2006,Head2014 with additive fixed effects,

the Poisson pseudo-maximum-likelihood estimator (PPML) for $\beta$ can be compared with our estimator. The PPML estimator $\hat{\beta}_{PPML}$ with two-way fixed effects is defined as the solution to:

where $a_{1},...,a_{n}$ satisfy

We can make a similar comparison as in the fixed effect estimator based on the results by Daveziez2021 and Menzel2021: $\hat{\beta}_{PPML}$ will be biased because of the misspecified errors, and the re-centered $\hat{\beta}_{PPML}$ is asymptotically normal at the rate of $\sqrt{n}$ in the non-degenerate case and can be non-Gaussian in the degenerate case.

Variance Estimation

\\ Since our estimator exhibits different asymptotic distributions depending on $\Sigma_{W\nu,1}$, it is desirable to have a variance estimator that adapts to the degeneracy.

First, we estimate $\Sigma_{W\nu,1}$. Define

where $\Delta\hat{\epsilon}_{ij}$ is a residual $\Delta Y_{ij}-\Delta W_{ij}'\hat{\beta}_{n}$. Then, we propose an estimator for $\Sigma_{W\nu,1}$ as

Next, we estimate $\Sigma_{W\nu,2}$ by

The following result shows consistency of these estimators and their usefulness in adaptive variance estimation.

We now propose our variance estimator as follows:

We can see that this estimator is adaptive to the degeneracy: When $\Sigma_{W\nu,1}$ is positive definite, since $n/(Nh_{n})=o(1)$,

as $n\to\infty$ by Proposition (ref) and Lemma (ref) in Appendix A. When $c_{W}'\Sigma_{W\nu,1}c_{W}=0$, since $nh_{n}c'\hat{S}_{WW}^{-1}\hat{\Sigma}_{W\nu,1}\hat{S}_{WW}^{-1}c=o_{p}(1)$ by Proposition (ref),

as $n\to\infty$.

Our variance estimator is adapted from the one provided in Graham2019 for a dyadic nonparametric density estimator. They show that this type of estimator can be adaptive to the "knife edge" case, where $nh_{n}$ is bounded from above and below asymptotically so that $Nh_{n}\sim n$. Here, we additionally show that the estimator is adaptive to the degeneracy by showing that the term involving $\hat{\Sigma}_{W\nu,1}$ decays fast enough to be negligible when the convergence rate is $\sqrt{Nh_{n}}$.

Bandwidth Selection

\\ From the asymptotic distributional approximation result in Theorem (ref), we can write down the mean squared error of our estimator (without negligible parts)

The optimal solution for minimizing this mean squared error with respect to $h_{n}$ is given by

We can estimate $h^{*}$ by the plug-in method. By Proposition (ref), we have a consistent estimator for the variance part. For the bias part, we use a pilot bandwidth given by

for some $\delta\in(0,\frac{2k+3}{4k+4})$ and $h>0$. Let $\hat{\beta}_{n,\delta}$ be our estimator calculated with $h_{n,\delta}$. We can check that this bandwidth satisfies $Nh_{n,\delta}^{2k+3}\to\infty$ and $nh_{n,\delta}^{2k+2}\to\infty$. Thus, by Theorem (ref),

as $n\to\infty$. By replacing $\beta$ by $\hat{\beta}_{n}$, calculated with $h_{n}=hN^{-\frac{1}{2k+3}}$, we have the following result:

Thus,

is a consistent estimator for $h^{*}$ by Propositions (ref) and (ref).

Bias Correction

\\ Notice that our estimator has the asymptotic bias of $\sqrt{h}\Sigma_{WW}^{-1}\Sigma_{W\lambda}$ in the case of degeneracy, $\Sigma_{W\nu,1}=0$ from Theorem (ref). If the bias is non-negligible, it distorts the coverage probability of the confidence interval. Correcting the bias part is desirable as it is generally unknown whether the degeneracy occurs. Fortunately, given the similar asymptotic distributional result as Kyriazidou1997 in the degenerate case, we can use her bias correction strategy as follows.

Note that $h_{n,\delta}^{-(k+1)}(\hat{\beta}_{n,\delta}-\beta)$ directly estimates the asymptotic bias from Theorem (ref). We can construct a bias-corrected estimator $\hat{\beta}_{n,bc}(\beta)$ by subtracting this bias estimator from the original estimator with suitable normalization: Let $r_{n,\delta}=N^{(1-\delta)/(2k+3)}$. The bias-corrected estimator is given by

We can check that this estimator is asymptotically unbiased regardless of the degeneracy: When $c_{W}'\Sigma_{W\nu,1}c_{W}>0$,

as $n\to\infty$. When $c_{W}'\Sigma_{W\nu,1}c_{W}=0$,

as $n\to\infty$. Thus, given the adaptivity of $\hat{\Sigma}$ to the degeneracy, we have

as $n\to\infty$ for an arbitrary non-zero vector $c\in\mathbb{R}^{q_{w}}$.

Then, we can construct the bias-corrected confidence interval as follows: Letting $\Phi_{1-\alpha/2}^{-1}$ be $1-\alpha/2$ quantile of the standard normal distribution, we have

where

The full inference procedure is summarized as follows:

First-step Estimator

\\ Remember that we want to estimate $\gamma$ from the selection equation or network formation process ((ref)):

This DGP can be interpreted as a panel discrete choice model as well as a network formation model. Estimators for discrete choice models such as Chamberlain1980, Manski1987, or Horowitz1992 can be candidates for estimating $\gamma$. Also, estimators for network formation models such as Graham2017 or Candelaria2020 can be applicable under additional conditions.

Whether those estimators can be used as our first-step estimator $\hat{\gamma}_{n}$ boils down to their convergence rates: Recall that Assumption (ref) requires that $\sqrt{Nh_{n}}(\hat{\gamma}_{n}-\gamma)=o_{p}(1)$, which implies that the first-step estimator needs to converge faster than $\hat{\beta}_{n}$. We can conjecture that, without additional conditions on $\eta_{ijt}$, the convergence rates of those estimators are $\sqrt{n}$ in worst cases due to the conditional dependence across dyads. Obviously, $\sqrt{n}$-rate is incompatible with Assumption (ref). In the following, we discuss what kind of additional conditions are needed to ensure Assumption (ref).

We may assume additive separability for $\eta_{ijt}$: $\eta_{ijt}=V_{it}+V_{jt}+V_{ijt}$, where conditionally on $\{\xi_{i}\}_{i=1}^{n}$, $(V_{i1},V_{i2}),i=1,...,n$ is independent, $(V_{ij1},V_{ij2}),1\leq i<j\leq n$ is independent, and both are mutually independent. This assumption is weaker than assuming $(\eta_{ij1},\eta_{ij2}),1\leq i<j\leq n$ is conditionally independent given $\{\xi_{i}\}_{i=1}^{n}$, where $V_{it}$ is treated as degenerate. With additional conditions, we can directly apply Graham2017's joint maximum likelihood estimator or Candelaria2020's semiparametric estimator, both of which leverage the cross-sectional variation in $d_{ijt}$ and $R_{ijt}$. We can show that in our setting (especially Assumptions (ref) and $\ref{asm:moments}$), the limiting networks are dense, which implies that both Graham2017 and Candelaria2020's estimators satisfy $\sqrt{N}(\hat{\gamma}_{n}-\gamma)=O_{p}(1)$ and Assumption (ref).

Alternatively, we may assume that $(\eta_{ij1},\eta_{ij2}),1\leq i<j\leq n$ is conditionally independent given $\{\xi_{i}\}_{i=1}^{n}$ and $\varphi(B_{i},B_{j})=B_{i}+B_{j}$. Graham2017 and Candelaria2020's estimators still satisfy Assumption (ref), but we can also show that Chamberlain1980's conditional logit estimator and Horowitz1992's smoothed maximum score estimator can satisfy Assumption (ref). Under the conditional independence assumption, the latter two estimators can be written in an asymptotically locally linear form where the corresponding influence function is indexed by $ij$ with $0$ covariances. Thus, the convergence rates are based on $N$ and Assumption (ref) can be satisfied depending on the tuning parameters.

Extension

Directed Graph with Multiple Periods

\\ In the above analysis, we restricted our attention to an undirected graph; the variables are all symmetric with respect to nodes (e.g., $Y_{ijt}=Y_{jit}$). Also, there were only two time periods, $t=1,2$. The extension to a directed graph case with $t=1,...,T\, (T\geq 2)$ is straightforward; Letting $\Delta_{st}A\equiv A_{s}-A_{t}$ denote the time difference between $s$ and $t$, we propose the following estimator:

All the results and their proofs are valid with some modification because we can always rewrite the double sum $\sum_{i=1}^{n}\sum_{j\neq i}A_{ij}$ as $\sum_{i<j}(A_{ij}+A_{ji})$ for any variables $\{A_{ij}\}$. We will use this version of the estimator in our empirical application.

Pairwise Fixed Effects

\\ In the model ((ref)) and ((ref)), all the fixed effects are node-wise. Since we are interested in coefficients on time-varying dyadic variables, it is possible to include pairwise fixed effects $A_{ij}$ and $B_{ij}$ in each equation, additionally to $A_{i},A_{j}$ and $B_{i},B_{j}$. Clearly, with pairwise fixed effects, the identification and estimator will be the same as with node-wise fixed effects since we are leveraging the time variation. Thus, a similar asymptotic analysis will also hold as long as $(A_{ij},B_{ij}),1\leq i<j\leq n$ are independently distributed conditionally on $\{\xi_{i}\}_{i=1}^{n}$.

Alternatively, we can also do away with the additive separability by incorporating node-wise fixed effects into pairwise ones:

where $\tilde{\tau}$ is some unknown function, $\tilde{A}_{i}$ is a node-wise fixed effect, and $\tilde{A}_{ij}$ is a pairwise fixed effect. We can impose a similar structure for $B_{ij}$. Again, the asymptotic analysis will hold as long as $(\tilde{A}_{ij},\tilde{B}_{ij}),1\leq i<j\leq n$ are conditionally independent. With a more general dependence structure, we could show a similar asymptotic result using Kojevnikov2020's central limit theorem for $\psi$-dependent data.

Sparsity

\\ Above, we argue that our model and assumptions imply that the limiting networks $D_{1}$ and $D_{2}$ are locally dense around $\Delta R_{ij}'\gamma\sim 0$. Thus, we limit our attention to cases where the number of dyads in the sample must be proportional to $N$. Our modeling is appropriate in some applications, such as trade or migration, where the number of dyads is rather dense. However, ours can be inappropriate for some applications where the networks are sparse such as employee-employer, bank-firm matched data (e.g., Abowd1999, jimenez2014).

We can accommodate sparse networks by the following modification; let us modify Assumption (ref) so that $\xi_{i},i=1,...,n$ are drawn from some distribution that is allowed to depend on $n$. For example, as argued in Graham2017, we can consider a distribution where the fixed effects are such that $\liminf_{1\leq i\leq n}B_{i}=-\infty$. Then, we can discuss identification and estimation with fixed $n$, and the moments of interest are all dependent on $n$. Especially, we can consider the sequence of networks such that $Pr(d_{121}d_{122}=1|\Delta R_{12}'\gamma=0)\to 0$ and $r_{n}Pr(d_{121}d_{122}=1|\Delta R_{12}'\gamma=0)=\Omega(1)$ for some $r_{n}\to\infty$ to incorporate sparsity. We do not pursue sparsity in this paper and leave it for future projects.

Simulation

To see the performance of the estimator, we conduct some simulation exercises. Consider the following data-generating process:

where

Note that $\beta=1$ and $\gamma=(1,1)'$. We have $\theta\in\{-0.3,-2.0,-3.0\}$ inside of $d_{ijt}$ to control for the fraction of zeros in the simulated data set:

We also change $\sigma\in\{0.0,1.0\}$ for $U_{it}$ so that $\sigma=0.0$ ($\sigma=1.0$) corresponds to the degenerate (non-degenerate) case.

As described above, we can interpret this data-generating process as a log-linearized version of the canonical gravity model (Head2014); by writing $\tilde{Y}_{ijt}$ as an observable outcome, we redefine the main equation as

We can take a log and recover the original model for a unit with $d_{ijt}=1$. This modeling allows a mass at $\tilde{Y}_{ijt}=0$, one important feature of dyadic data.

We conduct experiments for $n\in\{50,100,150,200\}$, $\theta\in\{-0.3,-2.0,-3.0\}$, and $\sigma\in\{0.0,1.0\}$, and iterate $2000$ times for each one. We calculate $\hat{\gamma}_{n}$ by Chamberlain1980's conditional logit estimator:

where $\mathcal{G}$ is a compact subset of $\mathbb{R}^{q_{r}}$ and

For $\hat{\beta}_{n}$, we use a biweight kernel for $K(\cdot)$, given by $K(x)=15/16(1-x^{2})^{2}\boldsymbol{1}\{|x|\leq 1\}$. This choice implies that we assume that the smoothness of the model is given by $k=2$. We set $\delta=0.4$ and $h=3.0$ and calculate each estimator and confidence interval according to the inference procedure discussed above.

For comparison, we calculate the fixed effect estimator $\hat{\beta}_{FE}$ given by ((ref)). The standard error is calculated by $\hat{\Sigma}$, with $K_{h_{n}}(\cdot)$ replaced by $1$. We also calculate the Poisson pseudo-maximum-likelihood (PPML) estimator $\hat{\beta}_{PPML}$ given by ((ref)). We compute $\hat{\beta}_{PPML}$ and its standard error by the penppml package in R penppml. The standard error is clustered at the node level, which is close to $\hat{\Sigma}_{WW}^{-2}\hat{\Sigma}_{W\nu,1}$ in our setting Graham2020.

The result is summarized in the following TABLE 1 and 2. In TABLE 1, we evaluate the three estimators by mean and median biases (MeanBias), root mean square error (RMSE) for $\sigma=0,1$. In TABLE 2, we compute 95% coverage probabilities (Coverage) of four different confidence intervals: $CI_{conv}$ (conventional CI from $\hat{\hat{\beta}}_{n}$ and $\hat{\Sigma}$), $CI_{bc}$ (bias-corrected CI given by $CI_{L,0.05}$ and $CI_{U,0.05}$), $CI_{FE}$ (conventional CI from $\hat{\beta}_{FE}$ and $\hat{\Sigma}$ with a flat kernel.), and $CI_{PPML}$ (conventional CI from $\hat{\beta}_{PPML}$ and its node-level clustered standard error).

From TABLE (ref) and (ref), we can see that our estimator performs better than the fixed effect estimator and the PPML estimator in terms of bias, which shows that the weights given by the first step estimator work well in eliminating the bias. Our estimator also outperforms the competitors regarding RMSE, which implies that the loss in precision is not severe. Our estimator also performs well even when there is a large fraction of zeros in $Y$ ($Pr(D_{ij1}\times D_{ij2})\sim 90\%$ when $\theta=-3.0$). There is little difference between $\sigma=0$ and $\sigma=1$ other than added variances in the estimators.

From TABLE (ref) and (ref), we can see that $CI_{bc}$ is close to 95% regardless of the degeneracy ($\Sigma_{W\nu,1}=0$ or $>0$) while the others are off from the targeted nominal coverage. This result confirms the effectiveness of the bias correction strategy as well as the adaptivity of our variance estimator, as claimed in Section 3.3. Also, it is notable to see that the bias correction is important for obtaining correct coverage probabilities even though the asymptotic bias is $0$ in the case of $\sigma=1.0$ so that $\Sigma_{W\nu,1}>0$ (Theorem (ref)) and $CI_{conv}$ would return an asymptotically correct coverage.

Empirical example

Background

\\ As a leading application of our model, consider Moretti2017. They study how state-level tax differences affect migration by top scientists in the U.S. Specifically, they estimate the following model implied by their economic theory:

where $P_{ijt}$ is the number of scientists migrating to state $j$ from state $i$ at year $t$, $\tau_{it}$ and $\tau_{it}'$ are personal and corporate taxes imposed in state $i$ at year $t$, $\gamma_{i}$ is a state fixed effect, and $u_{ijt}$ is an error term.

Note that if there is no migration from $i$ to $j$ at year $t$, $P_{ijt}=0$ and $log(P_{ijt}/P_{iit})$ is undefined. In Moretti2017's dataset, more than 70% of state-pairs exhibit no migration flow:

When running a regression, they are concerned with a potential sample selection bias stemming from these undefined outcomes. They argue that if the main regressors are not systematically associated with the probability of positive migration flows, the selection bias should be minimal. Running OLS on the linear probability model, they find little correlation between the main regressors and no flow. Re-estimating their model with our method provides a check on the validity of their argument and the appropriateness of using the linear probability model.

When applying our model to their context, we must consider what $R_{ijt}$ should be. Since Moretti2017's underlying theory is based on scientists' and firms' discrete choice, one consistent way to generate zero migration flows between some states is to consider endogenous choice sets as in Dube2021. Formally, we can write the choice set of representative scientists in state $i$ as $C_{it}=\{j\in\{1,...,51\}:d_{ijt}=1\}$. Here, $\{d_{ijt}\}$ represents the job-market network; if $d_{ijt}=1$, it is possible to move from $i$ to $j$, and vice versa. We can attribute the determinants of the network to the utilities and profits of scientists and firms, as well as the matching costs between the two parties. Such costs are not present in the structural equation if those costs are not compensated through wages; The structural equation consists of the determinants of log wage differences between two states. Thus in the selection equation ((ref)), in addition to $W_{ijt}$, we can include variables in $R_{ijt}$ that capture non-monetary matching costs between two states $i$ and $j$, which does not violate Assumption (ref) as $R_{ijt}$ satisfies the exclusion restriction.

Implementation

\\ For $W_{ijt}$, as in Moretti2017, we include the state-to-state differences in (i) an individual income average income tax rate (ATR) faced by a hypothetical taxpayer at 99% quantile of the national income distribution, (ii) the corporate tax rate (CIT), (iii) the investment tax credit (ITC), and (iv) the R&D tax credit (R&D credit). This is the same set of regressors as Moretti2017's baseline regression. For $R_{ijt}$, we use $W_{ijt}$ plus state-to-state difference in the logarithm of population (POP) and a dummy variable that indicates whether $i$ and $j$ share their governors' political parties (GOV). The additional variables in $R_{ijt}$ arguably measure non-monetary costs of connecting firms and workers in two states.

We implement the first step estimation as follows. We use the conditional logit estimator extended to a directed graph with multi-periods case, which is given by

where $\mathcal{G}$ is a compact subset of $\mathbb{R}^{q_{r}}$ and

TABLE (ref) reports the first step estimation result. We can see that the coefficients on the newly added variables $GOV$ and $POP$ deviate from zero, which implies that a part of the identification assumptions (Assumption (ref)) is satisfied. Also, for these two variables, the estimated coefficients imply that the job market network exhibits homophily; similar states are more likely to be connected.

For the second step estimator, we use our $\hat{\beta}_{n}$ defined above and extend it to the directed graph with multiple period cases, as discussed in Section 4. We use a biweight kernel for $K(\cdot)$ (so $k=2$), choose $h=3.0$ as an initial constant for pilot bandwidths, and use $\delta=0.4$ for calculating $h_{n,\delta}$. We extend and use $\hat{\Sigma}$ to calculate the standard error while taking into account the correlation across time (case 1). We also calculate the bias-corrected $95\%$-confidence interval by computing $CI_{L,0.05}$ and $CI_{U,0.05}$ as defined above. Also, we list $\hat{\beta}_{MW}$, an estimate from Moretti2017 (page 1883, TABLE 2A, specification (3)) and calculate the conventional $95\%$-confidence interval based on their standard errors. Note that $\hat{\beta}_{MW}$ is the fixed estimator, where its standard error is calculated by clustering across time and the origin, destination, and origin-destination pairs.

We summarize the result in TABLE (ref). We can see that $\hat{\beta}_{n}$ returns similar values as $\hat{\beta}_{WM}$, which claims the robust positive effect of income and corporate-related tax differences on migration. Thus, Moretti2017's estimates are not likely to be qualitatively affected by the sample selection effects. However, while Moretti2017's estimates are statistically significant at $5\%$ level, our confidence intervals show that all of the estimates are no longer statistically significant at that level except for ITC. Our insignificance result is driven by both the increase in standard errors \footnote{Our standard error is from $\hat{\Sigma}$, which takes fully into account the dependence among pairs that share origin and destination, such as California$\to$Wisconsin and New York$\to$California. Moretti2017's standard error calculation ignores such dependence structure.} and the asymptotic bias correction. Thus, our exercise shows that some of the results in Moretti2017 may not be robust to the presence of sample selection due to the endogeneity of the job market network.

Conclusion

This paper studies identification and inference of a panel dyadic data sample selection model. We show that Kyriazidou1997's identification strategy can be extended to our dyadic data setting, and we prove asymptotic normality of the proposed estimator.

Our estimator has some appealing properties. The distributional result implies that our estimator has the same convergence rates as the usual estimators used in practice in the non-degenerate case, and there is no loss of effective sample size for using our nonparametric type estimator. Also, our estimator is guaranteed to be asymptotically normal, while others can be non-Gaussian in the limit.

We also provide consistent estimators for asymptotic bias and variance that adapts to the degeneracy. Specifically, the bias-corrected confidence interval has an asymptotically correct size. Our simple simulation exercise confirms the validity of these estimators and highlights the importance of bias correction in both degenerate and non-degenerate cases.