Penalized Likelihood for Dyadic Network Formation Models with Degree Heterogeneity

Penalized Likelihood for Dyadic Network Formation Models with Degree Heterogeneity

\setcounter{footnote}{0}

Introduction

Empirical analysis of network formation routinely relies on dyadic models to relate the presence of a link to observable characteristics of the pair. Because link formation also depends on persistent agent specific attributes that are often unobserved by the researcher, it is standard practice to control for degree heterogeneity through sender and receiver fixed effects. In many empirical networks, however, this creates an estimation trilemma.

First, the fixed effects maximum likelihood estimator (MLE) frequently fails to exist in finite samples. Crucially, this existence failure is not confined to globally sparse networks, it routinely occurs in moderately dense data whenever local degree sparsity is present, such as agents sending no links, receiving no links, or connecting to everyone. In practice, researchers attempt to restore existence by trimming these boundary-degree nodes. However, in directed networks, this mechanical trimming is iterative and cascading. The resulting MLE is therefore computed on an endogenously selected subnetwork rather than the full observed network, inducing severe sample selection bias. Second, even when the standard MLE exists, the incidental parameter problem induces first order bias in the common parameters neyman1948consistent. Third, while analytical bias corrections exist, extending them to structural functions such as average partial effects (APEs) typically requires cumbersome higher-order derivative calculations.

We address these problems jointly with a penalized likelihood (PL) approach to dyadic network formation models with degree heterogeneity. The PL estimator remains feasible in network samples with boundary-degree nodes, and yields bias corrected estimates of both coefficients and average structural functions. Our leading specification is a directed network formation model with both unilateral incentives and bilateral reciprocity. Let $g_{ij}=1$ if agent $i$ forms a link to agent $j$, and let $g_{ij}=0$ otherwise. The payoff to agent $i$ takes the form

where $\alpha_i$ and $\gamma_j$ denote sender and receiver heterogeneity, and the mutual utility term captures reciprocity via bilateral dependence in link formation. Following the dynamic game approach of mele2017structural, we show that, under an ergodicity condition, a stochastic best response dynamic yields a tractable stationary likelihood and a unique equilibrium distribution. This foundation motivates the likelihood, and our PL estimator is built on it.

The reciprocal directed specification serves as our leading case because it structurally nests the standard dyadic models in the literature. Setting the mutual utility parameter to zero ($\rho=0$) yields the directed model without reciprocity jochmans2018semiparametric, dzemski2019empirical, yan2019statistical, hughes2025estimating, while equating the sender and receiver effects ($\alpha_i=\gamma_i$) under symmetric interactions recovers the undirected model graham2017econometric. We develop the estimator and asymptotic theory for the general reciprocal case, and the same penalization logic applies to these benchmark models as special cases.

This paper makes four econometric contributions. First, the PL approach guarantees the estimator exists in finite sample without trimming isolated or boundary-degree nodes. The penalty keeps the objective function well behaved when fixed effects drift toward the parameter boundary. This avoids the sample selection bias induced by repeated trimming the network data.

Second, the same likelihood-based penalty corrects the leading incidental parameter bias in the common parameters. The penalty is constructed from the block-diagonal part of the fixed-effects Hessian, with a \(2\times2\) sender-receiver Hessian block for each node in the reciprocal directed model. Our penalty is non-data-dependent and does not require the estimation of additional expectation terms. By contrast, an estimator-level analytical correction can also be derived in the reciprocal setting. However, it still requires existence of the unpenalized estimator and therefore does not resolve the finite-sample nonexistence problem. Conditional likelihood methods provide useful benchmarks in nested settings. In the reciprocal directed model, however, they rely on six-node hexiad cycles rather than tetrads or quadruples, which is computationally more demanding in large networks. In addition, because fixed effects are conditioned out, recovery of APEs is less direct.\footnote{For our reciprocal setting, we also develop a hexiad conditional likelihood approach, which generalizes the tetrad logit of graham2017econometric and the quadruple logit of jochmans2018semiparametric. See Appendix (ref) for details.}

Third, we extend the correction to APEs. This matters because empirical work typically reports partial effects, average link probabilities, and related structural functions rather than coefficients alone. After correcting the likelihood, we show the remaining leading bias in APEs comes from fixed effects uncertainty and can be removed with a simple adjustment based on the second derivatives of the APE function and the fixed-effects Hessian.

Beyond the finite-sample existence and bias correction, we establish asymptotic existence, consistency, and asymptotic normality without standard compactness assumptions. In the existing literature, assuming a compact parameter space for fixed effects implicitly forces expected node degrees to grow proportionally with the network size, leading to a dense limiting network. We break this limitation by establishing asymptotic results while allowing the true fixed effects to diverge at a logarithmic rate. This growth condition accommodates the degree sparsity often observed in large empirical networks.

Fourth, we approximate the inverse fixed-effects Hessian needed for sparse network asymptotics in the reciprocal directed model. In benchmark undirected, non-reciprocal directed, and two-way fixed effects panel models, a diagonal approximation to the fixed-effects inverse Hessian is often sufficient fernandez2016individual,hughes2025estimating. Reciprocity breaks this simplification because sender and receiver effects are coupled within each node and across reciprocal dyads. We derive a closed form hybrid approximation that combines \(2\times2\) node level Hessian blocks with a low-rank correction for two aggregate weak directions. The approximation remains valid when the network grows and some link probabilities shrink. It provides the uniform bounds used for consistency, asymptotic normality, and the APE correction without compactness assumptions on the fixed effects.

We then apply the method to the global textiles trade network. This network is not globally sparse. Yet the standard fixed-effects MLE is still not available on the observed network. 11 countries have zero in-degree, and trimming them induces 5 major economies (France, Germany, Italy, Japan, and the UK) to become full out-degree countries in the remaining network. The application shows that PL estimation remains feasible on the full observed network rather than on a selected subnetwork.\footnote{To facilitate empirical use, we provide a Python package, NetworkFm, that implements the PL estimator and the nested models used in the paper: \url{github.com/zizhongyan/NetworkFm}.}

Related literature. This paper connects to the literatures on fixed-effects estimation in nonlinear panel and network formation models. In nonlinear panel models with two-way fixed effects, bias correction for parameters and average partial effects is extensively studied fernandez2016individual, jochmans2019likelihood, chen2021nonlinear, leng2025debiased, yan2026robust.\footnote{Because the PL correction acts on the likelihood itself, it preserves the extremum structure of the objective function. This makes likelihood-based inference a natural extension, in line with recent results for nonlinear panel models in leng2025debiased and yan2026robust. We do not develop that theory here.} In the network econometric literature, graham2017econometric, jochmans2018semiparametric,dzemski2019empirical,hughes2025estimating,li2025bagging develop fixed-effects estimators. From a statistical perspective, yan2016asymptotics,yan2019statistical establish asymptotic properties for dyadic network models under diverging parameter spaces. ma2021determining identify latent community structures via random block models rather than individual fixed effects. We build on this literature by simultaneously accommodating reciprocity, ensuring finite-sample existence without network trimming, and providing a unified likelihood-based correction for both coefficients and APEs.

A related literature studies network formation under broader strategic interactions and network externalities, including miyauchi2016structural, mele2017structural, boucher2017my, de2018identifying, leung2019inference, sheng2020structural, hoshino2022pairwise, ridder2025two, leung2026normal, menzel2026strategic, and related work. We take a narrower focus on dyadic dependence generated by reciprocity, which yields a tractable estimator for sparse directed networks with unrestricted degree heterogeneity.

Outline of the paper. Section (ref) introduces the reciprocal model and the likelihood foundation. Section (ref) presents the penalized likelihood estimator and the APE correction. Section (ref) establishes the existence results and the large sample theory. Section (ref) reports the Monte Carlo evidence. Section (ref) presents the trade application. The final section concludes. Proofs and additional computational details are collected in the appendices.

Dyadic Network Formation Model and Likelihood Foundation

Payoff Specification and Nested Models

Consider a network of $n$ agents (e.g., individuals, firms, or countries). The observed network is an $n\times n$ adjacency matrix $G=(g_{ij})$, where $g_{ij}\in\{0,1\}$ indicates whether agent $i$ forms a directed link to agent $j$. In the trade network, for example, a directed link indicates that one country exports to another. We rule out self-links, so $g_{ii}=0$ for all $i$. Let $\mathcal{U}_{ij}$ denote the payoff to agent $i$ from establishing a directed link to agent $j$. We propose a transferable utility function that incorporates both directed utility and mutual reciprocity:

where the first directed utility term captures unilateral incentives to form a link from $i$ to $j$, while the second mutual utility term captures reciprocity through bilateral interdependence. $\alpha_i$ and $\gamma_j$ represent sender and receiver effects, respectively; $X_{ij}$ are covariates that enter the directed utility, and $Z_{ij}$ are covariates that enter the mutual component. The micro-foundation requires the covariates $Z_{ij}$ driving the reciprocal effect to be symmetric across dyads ($Z_{ij} = Z_{ji}$), whereas the directed covariates $X_{ij}$ can be asymmetric. We assume both $X_{ij}$ and $Z_{ij}$ include a constant term to capture baseline linking utilities.

The payoff specification in (ref) provides a unified framework that encompasses two benchmark models that are standard in the literature:

1) Directed networks without reciprocity. If we assume that link formation relies purely on unilateral incentives, we shut down the mutual utility component by restricting $\rho = 0$. The payoff collapses to

which yields the standard directed model with degree heterogeneity studied by jochmans2018semiparametric,yan2019statistical,dzemski2019empirical, among others. Under this restriction, agent $i$'s decision to link with $j$ is entirely independent of agent $j$'s linking decision.

2) Undirected networks. An undirected network requires that connections are symmetric by construction, so that $g_{ij} = g_{ji} \equiv g_{ij}^*$. The same specification also arises as a special case of the directed model when link payoffs depend only on reciprocity. In this case, the dyadic outcome is whether two agents are mutually linked, which can be written as $g_{ij}g_{ji}\equiv g_{ij}^*$. Under symmetry, sender and receiver heterogeneity are no longer separately distinguished and collapse to a single node-specific effect, so that $\alpha_i=\gamma_i$. The payoff reduces to

which is the undirected specification studied by graham2017econometric.

In Section (ref), we develop the penalized likelihood estimator for the generalized reciprocal model in (ref) as the leading case, and our estimation and bias-correction methods automatically accommodates these benchmark cases.

Stationary Equilibrium and Likelihood

In the directed model without reciprocity (ref), each link decision is unilateral, so the unique equilibrium is immediate. In the setting (ref), the dyad is summarized by a single symmetric link indicator, so the unique equilibrium is again straightforward.

However, once reciprocity is introduced to the directed model (ref), the pair $(g_{ij},g_{ji})$ becomes jointly determined. The presence of the mutual utility term $g_{ij}g_{ji}Z_{ij}'\rho$ introduces strategic interdependence between agent $i$ and agent $j$. As graham2020network notes, combining both unilateral incentives and bilateral reciprocity in cross-sectional network data requires additional structure. It is also well documented that in static, simultaneous-move games with such bilateral interactions, the network formation process generally suffers from multiple equilibria lewbel2007coherency, ciliberto2009market. hoshino2022pairwise highlights that such multiple equilibria render an incomplete model.

To bypass this problem and obtain a well-defined likelihood, we follow the dynamic game approach of mele2017structural. We model the observed cross-sectional network as the long-run realization of a stochastic best-response dynamic process blume1993statistical. Suppose that in each period $t$, agent $i$ meets $j$ and decides whether to update the connection $g_{ij}^t$ to maximize utility, taking $j$'s previous connection status $g_{ji}^{t-1}$ as given. Before updating, agent $i$ receives an idiosyncratic preference shock $\epsilon^t$. Agent $i$ sets $g_{ij}^t = 1$ if and only if $\mathcal{U}_{i}(g_{ij}^t = 1, g_{ji}^{t-1}, X_{ij}, Z_{ij}) + \epsilon_1^t \ge \mathcal{U}_{i}(g_{ij}^t = 0, g_{ji}^{t-1}, X_{ij}, Z_{ij}) + \epsilon_0^t$. We impose the following standard conditions on this dynamic process:

Under Assumption (ref), the network formation process operates as a potential game monderer1996potential. With strict positivity of meeting probabilities, the resulting Markov chain of network sequences is ergodic. Ergodicity ensures that the dynamic process converges to a unique stationary distribution, independent of the initial network configuration. This allows us to characterize the equilibrium distribution of the network as a standard exponential family model:

Proposition (ref) provides the equilibrium distribution that is essential for our econometric strategy built in Section (ref). Because the model rules out external dependence beyond the dyad $(i,j)$, the unique stationary distribution for the entire network of $n$ agents is simply the product of the dyadic distributions: $\prod_{i=1}^n \prod_{j>i}^n \pi(g_{ij}, g_{ji} | \cdot)$.

This dyadic independence yields a likelihood with a tractable denominator in (ref) and keeps estimation within the frequentist MLE framework. By contrast, models with broader link interdependence typically lead to an intractable denominator and require Bayesian MCMC methods.\footnote{To overcome the estimation difficulty, mele2017structural, mele2020does applies the Bayesian double Metropolis-Hastings algorithm based on murray2012mcmc to avoid the calculation of the denominator. mele2023approximate developed a variation inference algorithm to estimate a simplified version of the model.}

Penalized Likelihood Estimators

This section details the estimation methodology. We first construct the log-likelihood function to show why standard MLE fails in sparse networks. We then introduce a PL estimator that restores finite-sample existence without data trimming and simultaneously corrects incidental-parameter bias.\footnote{In the context of MLE, the existence of an estimator means that the supremum of the (log-)likelihood function can be achieved using a finite set of parameter values.} Finally, we provide a computationally tractable procedure to estimate bias-corrected APEs.

Likelihood and Existence

We adopt a fixed effects approach for estimation, treating $\alpha$ and $\gamma$ as parameters to be jointly estimated with the common parameters $\theta = (\beta', \rho')'$. We collect $2n-2$ fixed effects in the vector ordered by node: $\lambda=(\alpha_1,\gamma_1,\alpha_2,\gamma_2,\ldots,\alpha_{n-1},\gamma_{n-1})'$, subject to the location normalizations $\alpha_n = \gamma_n = 0$.\footnote{Equivalent normalizations are possible. Our choice keeps constant terms in both $X_{ij}$ and $Z_{ij}$, to capture baseline directed and reciprocal utility. Section (ref) returns to identification and population maximization problem.} Using the stationary distribution in Proposition (ref), the log-likelihood for the observed network $G$ is

where $\ell_{ij}\left(\theta,\alpha_i,\gamma_j,\alpha_j,\gamma_i \right):= \ln\left[\pi(g_{ij},g_{ji}\mid \theta,\alpha_i,\gamma_j,\alpha_j,\gamma_i)\right]$.

It is useful to rewrite (ref) in exponential family form. Let the out-degree and in-degree of agent $i$ be $d_i:=\sum_{j\neq i}^n g_{ij}$ and $b_i:=\sum_{j\neq i}^n g_{ji}$, respectively. Define the corresponding degree sequences $d=(d_1,\ldots,d_{n-1})'$ and $b=(b_1,\ldots,b_{n-1})'$. The log-likelihood can then be expressed as $$ \ell(\theta,\lambda)=\textstyle\sum_{i=1}^n\sum_{j\neq i}^n g_{ij}X_{ij}'\beta + \sum_{i=1}^n\sum_{j> i}^n g_{ij} g_{ji}Z_{ij}'\rho + d'\alpha + b'\gamma -\sum_{i=1}^n\sum_{j\neq i}^n \mathcal{C}_{ij}(\theta,\lambda), $$ where $\mathcal{C}_{ij}(\theta,\lambda)$ is the log partition function.\footnote{The log partition function is $ \mathcal{C}_{ij}(\theta,\lambda)=\frac{1}{2}\ln\big(1+\exp(B_{ij})+\exp(B_{ji})+\exp(B_{ij}+B_{ji}+C_{ij})\big), $ where $B_{ij}$ and $C_{ij}$ defined in Proposition (ref). Note that $\mathcal{C}_{ij}(\theta,\lambda)$ is symmetric in $i$ and $j$. }

This exponential family representation explicitly isolates the sufficient statistics. For the common parameters $(\beta,\rho)$, they are $\sum_{i\neq j} g_{ij}X_{ij}$ and $\sum_{i<j} g_{ij}g_{ji}Z_{ij}$. For the fixed effects $\lambda$, the sufficient statistics are the out-degree and in-degree sequences $(d', b')'$. This representation is convenient because the existence problem can be stated directly in terms of degree sequences.

Let $(\widehat\theta, \widehat\lambda) = \mathop{\rm arg\,max}_{\theta,\lambda} \ell(\theta,\lambda)$ denote the fixed effects MLE, provided it exists. We assume that the common parameter space is compact and the covariates have bounded support. Under these standard regularity conditions, the finite-sample existence of the MLE dependent entirely on the fixed effects $\lambda$. More importantly, we do not impose compactness on the parameter space of $\lambda$, because doing so would rule out the slowly diverging heterogeneity that is central to the sparse network asymptotics studied later in the paper.

For exponential family models, existence of the fixed effects MLE is equivalent to an interior condition on the sufficient statistics brown1986fundamentals. In the present model, this means that the observed bi-degree sequence $(d', b')$ must lie in the interior of the convex hull generated by all possible bi-degree sequences:

Empirical dataset routinely violate this convex-hull condition. The most obvious violations occur when an agent is an “isolated node”---sending no links ($d_i = 0$) or receiving no links ($b_i = 0$)---or a “fully connected node” ($d_i = n-1$ or $b_i = n-1$). To circumvent this issue, researchers often trim such nodes until the remaining sample becomes estimable. In directed networks, however, this procedure may need to be repeated, because deleting one node can create new null in-degrees or out-degrees in the remaining graph. The standard MLE is then based on an endogenously selected subnetwork rather than on the observed one. In the trade application of Section (ref), for example, MLE is only computable after trimming countries with zero in-degree. This is precisely the problem that the penalized likelihood estimator is designed to avoid.

Nevertheless, bounding the nodal degrees strictly away from $0$ and $n-1$ is merely a necessary condition for existence, not a sufficient one. As rinaldo2013maximum demonstrate, there are still numerous other cases where the MLE is undefined, and these cases can be difficult to detect in practice.\footnote{rinaldo2013maximum focus on the undirected $\beta$-model and systematically enumerate graph configurations on 4--6 nodes where MLE is undefined, despite all node degrees being strictly bounded away from zero and $n-1$.}

Even when the MLE exists in a sufficiently dense network, it still suffers from the incidental parameter problem. Because the dimension of $\lambda$ grows with the sample size $n$, estimation noise in $\widehat{\lambda}$ contaminates the concentrated likelihood, inducing a first-order bias of $\mathcal{O}(n^{-1})$ in $\widehat{\theta}$, which subsequently distorts structural functions such as APEs.

Section (ref) establishes the asymptotic existence properties, and shows that the estimator can exist asymptotically if the true fixed effects diverge sufficiently slowly. In finite samples, however, non-existence remains a severe barrier in empirical practice. The penalized likelihood estimator introduced next is designed specifically to resolve this finite-sample failure and the incidental parameter bias simultaneously, without data trimming.

Penalized Likelihood and Bias Reduction

We construct the penalty term using the block-diagonal extraction of the negative fixed-effects Hessian. Let ${D}(\theta,\lambda) = \operatorname{blkdiag}(D_1(\theta,\lambda), \dots, D_{n-1}(\theta,\lambda))$ denote this extraction from $-H_{\lambda\lambda}(\theta, \lambda) = - \partial_{\lambda\lambda'} \ell(\theta, \lambda)$. For our leading reciprocal directed model, each block $D_i(\theta, \lambda)$ is a $2 \times 2$ matrix corresponding to the sender and receiver effects $(\alpha_i, \gamma_i)$ of node $i$. We define the penalty via the log-determinant:

yielding the penalized likelihood estimator $(\widehat{\theta}_{\mathrm{PL}}, \widehat{\lambda}_{\mathrm{PL}}) = \mathop{\rm arg\,max}_{\theta,\lambda} \left[ \ell(\theta, \lambda) + \eta(\theta, \lambda) \right]$.

A compelling theoretical advantage of this formulation is that it adapts to the nested benchmark models discussed in Section (ref). In a directed model without reciprocity ($\rho=0$), link formation depends purely on unilateral incentives, rendering the cross-partial derivatives $\partial_{\alpha_i \gamma_i}\ell$ identically zero. The block $D_i(\theta, \lambda)$ then naturally collapses to a strictly diagonal matrix, reducing the penalty to the standard sum of univariate log-Hessians: $\frac{1}{2} \sum_{i=1}^{n-1} [\ln(-\partial_{\alpha_i \alpha_i} \ell) + \ln(-\partial_{\gamma_i \gamma_i} \ell)]$.

Similarly, in an undirected network where degree heterogeneity is captured by a single node-specific parameter ($\alpha_i = \gamma_i$), the block $D_i$ degenerates to a scalar. The penalty then mathematically simplifies to the exact adjustment for one-way fixed effects: $\frac{1}{2} \sum_{i=1}^{n-1} \ln(-\partial_{\alpha_i \alpha_i} \ell)$. Equation (ref) thus provides a completely unified bias-reduction framework across the entire family of dyadic models.

The mechanism of $\eta(\theta,\lambda)$ is twofold. First, it restores finite-sample existence. When any sender or receiver effect drifts to $\pm\infty$ (e.g., due to an isolated node), the corresponding diagonal Hessian term collapses toward zero, so $\eta(\theta,\lambda)\to -\infty$. The penalized log-likelihood therefore cannot be maximized at the boundary. Theorem (ref) in the next section formalizes this argument. Second, the same penalty removes the leading bias in the common parameters. Intuitively, the score contribution from $\eta(\theta,\lambda)$ offsets the leading bias in the score for $\theta$. Section (ref) makes this statement precise. Moreover, the penalty is non-data-dependent. It is determined entirely by the model-implied probabilities and can be evaluated directly.

We compare the PL estimator with two alternative estimators for the reciprocal model:

1) Relation to Estimator-level Correction (EC). For comparison, it is useful to note that the reciprocal model also admits an analytical estimator-level correction when the fixed effects MLE exists. Using the bias expansion derived in Section (ref), one can form $$ \widehat\theta_{\mathrm{EC}} = \widehat\theta - \tfrac{1}{n-1}\widehat{\mathcal{I}}_n^{-1}\widehat{ \mathscr{B}}_n, $$ where $\widehat{\mathcal{I}}_n$ and $\widehat{ \mathscr{B}}_n$ are consistent estimates of the information matrix and the leading bias term in the score, respectively. However, $\widehat\theta_{\mathrm{EC}}$ is unavailable whenever the underlying fixed effects MLE does not exist. We therefore treat $\widehat\theta_{\mathrm{EC}}$ as a benchmark estimator.

2) Relation to Conditional Logit. Conditional likelihood methods provide another benchmark. In the undirected case, graham2017econometric established the tetrad logit method, while in the directed case without reciprocity, jochmans2018semiparametric proposed a quadruple logit. In our reciprocal case, a conditional likelihood construction can in principle be built from six-node cycles, which we discuss in Appendix (ref) of Supplemental Material. For the reciprocity setting that motivates this paper, the PL approach is computationally simpler.

Average Partial Effects

Empirical network research typically interested in structural functions, such as partial effects or counterfactual link probabilities, not only in coefficients. Let $ p_{ij}(\theta,\lambda) = \Pr(g_{ij}=1\mid X,Z,\theta,\lambda) , $ denote the model-implied directed link probability, and let $W_{ij,k}$ denote a generic dyadic regressor, which may be a component of either $X_{ij}$ or $Z_{ij}$, and let $\Delta_{ij}(\theta,\lambda)$ denote the corresponding dyad-specific effect on $p_{ij}(\theta,\lambda)$.

For a continuous regressor $W_{ij,k}$, we define $ \Delta_{ij}(\theta,\lambda) = \frac{\partial p_{ij}(\theta,\lambda)}{\partial W_{ij,k}}. $ For a binary regressor $W_{ij,k}\in\{0,1\}$, we define $ \Delta_{ij}(\theta,\lambda) = p_{ij}(\theta,\lambda; W_{ij,k}=1) - p_{ij}(\theta,\lambda; W_{ij,k}=0), $ holding all other covariates fixed. When $W_{ij,k}$ is a component of $Z_{ij}$, the partial effect is understood as a dyad-level change applied symmetrically to $(i,j)$ and $(j,i)$, since $Z_{ij}=Z_{ji}$ in the reciprocal specification. In the trade application, these partial effects summarize how bilateral frictions affect unilateral and mutual trade links. We define the APE by $$ \Delta(\theta,\lambda) = \tfrac{1}{n(n-1)} \textstyle\sum_{i\neq j}^{n} \Delta_{ij}(\theta,\lambda). $$

The plug-in APE estimator, $\Delta(\widehat\theta_{\mathrm{PL}},\widehat\lambda_{\mathrm{PL}}) = \frac{1}{n(n-1)} \sum_{i\neq j} \Delta_{ij}(\widehat\theta_{\mathrm{PL}},\widehat\lambda_{\mathrm{PL}})$, suffers from asymptotic bias. As established in the two-way fixed effects nonlinear panel literature fernandez2016individual, yan2026robust, the first-order bias in APEs stems from three sources: (i) bias in the common parameters; (ii) bias in the fixed-effects point estimates; and (iii) variance of the fixed-effects estimates. Evaluating the APE at the PL estimates thus automatically eliminates the first two sources of bias. To remove the remaining first-order bias, we define the following bias-corrected APE estimator

This correction uses only the second-order derivative of the APE function and the fixed effects Hessian. No third-order derivatives are required.

This is one of the practical advantages of the PL approach. Estimator-level APE corrections can also be constructed, including for the reciprocal model, but they require more higher order derivative calculations. By contrast, (ref) keeps coefficient correction and APE correction within the same likelihood-based framework. The same logic applies to the nested directed and undirected benchmark models. Section (ref) establishes the asymptotic validity of (ref), and Section (ref) shows that it performs well in finite samples.

Asymptotic Theory for Penalized Likelihood Estimator

Notation and Assumptions

We first introduce notation used throughout the theoretical analysis.

For $a_1,a_2\in\{0,1\}$, let $p_{ij}^{(a_1a_2)}(\theta,\lambda):=\Pr(g_{ij}=a_1,g_{ji}=a_2\mid \theta,\lambda),$ and define the marginal link probability $p_{ij}(\theta,\lambda):=\Pr(g_{ij}=1\mid \theta,\lambda)=p_{ij}^{(11)}(\theta,\lambda)+p_{ij}^{(10)}(\theta,\lambda).$ When no confusion can arise, we suppress arguments evaluated at the true values $(\theta_0,\lambda_0)$, for example, we write $\ell=\ell(\theta_0,\lambda_0),$ and $p_{ij}=p_{ij}(\theta_0,\lambda_0).$

Let $H_{\theta\theta}=\partial_{\theta\theta'}\ell$, $H_{\theta\lambda}=\partial_{\theta\lambda'}\ell$ and $H_{\lambda\lambda}=\partial_{\lambda\lambda'}\ell$, where $\partial_x f$ denotes the derivative of $f$ with respect to $x$, and additional subscripts indicate higher-order derivatives. We write the Hessian of the log-likelihood as $$ H=

, \qquad H_{\theta\theta}=

, \qquad H_{\lambda\lambda}=

. $$

Let $\theta_0=(\beta_0',\rho_0')'$, $\lambda_0=(\alpha_{0,1},\gamma_{0,1},\alpha_{0,2},\gamma_{0,2},\ldots,\alpha_{0,n-1},\gamma_{0,n-1})'$ denote the true parameter vector under the identifying restriction $\alpha_n=\gamma_n=0$. An equivalent normalization is $\sum_{i=1}^n \alpha_i=\sum_{i=1}^n \gamma_i.$ We adopt $\alpha_n=\gamma_n=0$ because it allows both $X_{ij}$ and $Z_{ij}$ to retain a constant term. In the present model, the constant in $X_{ij}$ has a natural interpretation as baseline directed utility, while the constant in $Z_{ij}$ captures baseline reciprocal utility. Under this identification, the true parameter values are defined as the solution to the population maximization problem $$(\theta_0,\lambda_0)=\mathop{\rm arg\,max}\limits_{\theta\in\Theta,\lambda\in\Lambda_n}\mathbb E_0\bigl[\ell(\theta,\lambda)\bigr], $$ where $\mathbb E_0$ denotes expectation under the true distribution of the data conditional on the exogenous covariates.

Throughout, $\Theta$ is compact, while $\Lambda_n\subset\mathbb{R}^{2n-2}$ is closed but need not be bounded and hence not compact. The support of the covariates is bounded. We do not impose compactness on $\Lambda_n$, the asymptotic analysis later allows $\|\lambda_0\|_\infty$ to grow at a logarithmic rate. This accommodates sparse networks in which the number of observed links per agent does not grow proportionally with the network size.\footnote{A compact parameter space for the fixed effects would effectively force nodal degrees to grow proportionally with network size in the limit, thereby inducing a dense network asymptotically. This restriction can be strong for many empirical networks of interest graham2024sparse.} We collect the baseline conditions below:

Assumption (ref) is standard for nonlinear estimation problems newey1994large. The only point that differs from many fixed-effects formulations is that $\Lambda_n$ is assumed closed, but not compact. Condition (iv) is a higher level condition on the uniqueness of $\theta_0$ and $\lambda_0$, provided that they exist.\footnote{Since the model belongs to the exponential family, Assumption (ref)(iv) can equivalently be formulated as a rank condition on the sufficient statistic for large $n$. By Corollary 6.16 in lehmann1983theory, full rank of the sufficient statistic implies that the model is minimal, and guarantees uniqueness of the parameter vector.}

Existence

Two different existence questions arise in this paper. The first is finite-sample existence. In a realized network, the fixed effects MLE may fail to exist because the observed bi-degree sequence lies on the boundary of the convex support, as discussed in condition (ref) of Section (ref). The second is asymptotic existence. Even when an estimator is well defined in a given sample, its large sample behavior still depends on how quickly the true fixed effects are allowed to grow with the network size. The penalty resolves the first problem directly. The second remains a large-sample issue. We discuss finite-sample existence first.

Theorem (ref) is a finite sample result. It does not require the convex-hull condition (ref), and therefore does not require trimming isolated or fully connected nodes before estimation. This is the practical existence property that distinguishes the penalized estimator from estimator-level corrections that are defined only when the underlying MLE exists.

Finite-sample existence alone does not imply asymptotic existence. We still need to specify how the true fixed effects behave as the network grows. Following yan2016asymptotics,yan2019statistical, we allow the maximum entry of the true fixed effects $\|\lambda_0\|_\infty$ diverges at a rate proportional to $\tau \ln n$ for a sufficiently small constant $\tau>0$. While a compact parameter space for the fixed effects implies a dense limiting network, this logarithmic growth condition relaxes this restriction. It does not force the whole limiting network to be sparse, but it allows link probabilities to shrink as $n$ increases.\footnote{To see the intuition, consider the link probability $p_{ij} \approx e^{\alpha_i + \gamma_j} / [1 + e^{\alpha_i + \gamma_j}]$. Under the compactness assumption, $p_{ij}$ is strictly bounded away from zero to one. Expected nodal degrees such as $\sum_{j \ne i} p_{ij}$ grow at rate $\mathcal{O}(n)$ and generate a dense limiting network. Under $\|\lambda_0\|_\infty \le \tau \ln n$, the link probabilities can instead decay at a polynomial rate $p_{ij}\approx e^{-2\tau \ln n} = n^{-2\tau}$ See graham2024sparse for an in-depth discussions on the asymptotically sparsity.}

In the global trade application studied in Section (ref), expanding the network sample $n$ typically involves including smaller or more remote economies. Consequently, the probability that two countries form a trade link may become small. A standard compactness assumption would force two countries to trade when sample increases. Our asymptotic framework accommodates this empirical concern.

Approximation of the Inverse Hessian

A central technical challenge in this paper is deriving an explicit approximation for the inverse fixed-effects Hessian, $-H_{\lambda\lambda}^{-1}$. This task is substantially more demanding in the reciprocal directed model than in standard nested models.

In the undirected model, the directed model without reciprocity, and standard two-way fixed-effects panels, the fixed-effects Hessian exhibits a strictly diagonal or near-block-diagonal structure. Consequently, an approximation based on its diagonal is often adequate: $-H_{\lambda\lambda}^{-1}\approx -\mathrm{diag}(H_{\lambda\lambda})^{-1}$ fernandez2016individual, hughes2025estimating.\footnote{ Related graph-based fixed effects results show that precision depends on connectivity and on the small eigenvalues of the graph Laplacian chung1997spectral,jochmans2019fixed. } Once reciprocity is introduced, however, this simplification breaks down. Because an agent's outgoing and incoming links are mutually dependent, the off-diagonal entries within the node-specific blocks ($H_{\alpha_i\gamma_i}$) and the cross-node blocks ($H_{\alpha_i\gamma_j}$) are strictly non-zero and scale as $\mathcal{O}(n)$. As a result, the standard diagonal approximation no longer applies.

To bypass the computationally prohibitive high-dimensional inversion, we derive a closed-form matrix $S$ that analytically approximates the inverse Hessian. For each $i=1,\ldots,n-1$, we capture the node-specific information via the $2 \times 2$ matrix

and define the block diagonal matrix $D:=\operatorname{blkdiag}\bigl(D_1,\ldots,D_{n-1}\bigr)$. To account for global network movements, we define the aggregate basis vectors $u_+:=(1,\ldots,1)'$ and $u_-:=(1,-1,\ldots,1,-1)'$, stacked as $U:=[u_+,u_-]$. We then define the approximate inverse as

The approximation \(S\) is explicit and easy to compute. The first term, $D^{-1}$, requires only the inversion of $n-1$ independent $2\times 2$ blocks. In the trade network application, it accommodates the local dependence between outward and inward trade propensities. The second term is a low-rank aggregate correction requiring the inverse of only a single $2 \times 2$ matrix, $U'[-H_{\lambda\lambda}]U$. This decomposition provides a computationally tractable closed-form substitute for the full $(2n-2)\times(2n-2)$ matrix inversion.

To formalize the approximation precision of $S$ under sparse network asymptotics, we define a standardized Hessian $\widetilde V:= D^{-1/2} [-H_{\lambda\lambda}] D^{-1/2}$. In this standardized space, the aggregate subspace is spanned by $\widetilde U:= D^{1/2}U$. Let $P:= \widetilde U\bigl(\widetilde U'\widetilde U\bigr)^{-1}\widetilde U'$ denote the orthogonal projector onto the column space of $\widetilde U$, and let $P^\perp:=I-P$ be the residual projector.

Assumption (ref) formally describes the algebraic structure of the standardized Hessian. Condition (i) isolates the two aggregate directions $u_{+}$ and $u_{-}$. In a global trade network, $u_+$ represents a network-wide parallel shift in trade density, while $u_-$ raises outward trade propensities and lowers inward trade propensities by the same amount, or vice versa. Since these global movements are weakly informative for identifying relative country differences, they are absorbed by the low-rank projection $U$.

Condition (ii) imposes regularity on the remaining directions. After the two aggregate movements are removed, the inverse standardized Hessian is close to the identity on the orthogonal subspace. In the global trade example, countries differ fundamentally in their numbers of export and import partners. \(D_i\) measures the information available for country \(i\)'s exporter and importer effects. Countries with more trading partners contribute more variations, while countries close to isolation contribute less. The condition allows such degree heterogeneity. After this degree adjustment, condition (ii) requires that the likelihood still has enough local variation to distinguish countries' relative outward and inward trade propensities.

The rate in condition (ii) is compatible with Assumption (ref). Under the logarithmic growth condition used below, \(\|\lambda_0\|_\infty\le \tau\log n\), the right-hand side is of order $n^{-1+4\tau}$. Hence the condition does not require the fixed effects parameter space to be compact.

The next lemma gives the rate bounds for the inverse-Hessian approximation used below. We state the result for a generic vector \(\mathbf f\) because the same calculation is applied to the fixed effects score, Hessian-score products, and higher order remainders. The lemma uses two restrictions on \(\mathbf f\). The largest coordinate of \(\mathbf f\) and the aggregate sender-receiver difference are of the same order. This formulation avoids applying a matrix norm to \(-H_{\lambda\lambda}^{-1}-S\), which would give a bound too crude for the subsequent asymptotic expansions. In Lemma (ref), the first bound controls the approximation error, and the second controls the bound of \(S\mathbf f\).

Large-Sample Properties

We now turn to the large-sample properties of the estimator. For each $\theta$, let $$ \widehat\lambda(\theta)=\mathop{\rm arg\,max}_{\lambda\in\Lambda_n}\ell(\theta,\lambda), \qquad \bar\lambda(\theta)=\mathop{\rm arg\,max}_{\lambda\in\Lambda_n}\mathbb E_0[\ell(\theta,\lambda)]. $$ The common parameter estimator can then be viewed as: $ \widehat\theta = \mathop{\rm arg\,max}_{\theta\in\Theta} \ell\bigl(\theta,\widehat\lambda(\theta)\bigr). $ The incidental parameter problem arises because each component of $\widehat\lambda(\theta)$ is estimated from only $n-1$ observations, so the noise in $\widehat\lambda(\theta)$ contaminates the concentrated score for $\theta$. Hence the limit distribution of $\widehat{\theta}$ does not center at $\theta_0$, unless the quantity $\widehat{\lambda}(\theta)$ was replaced by its infeasible population quantity $\bar{\lambda}(\theta)$. This is the same mechanism that generates first-order bias in nonlinear panel models with two-way fixed effects.

The first step is to establish asymptotic existence and consistency. The proof combines a Newton iteration for $\widehat\lambda(\theta)$ with the approximation in Lemma (ref). The log growth condition on $\|\lambda_0\|_\infty$ guarantees that the iteration remains in a neighborhood where the Hessian approximation is valid.

The next result gives the asymptotic distribution.

Theorem (ref)(i) makes the incidental parameter problem explicit. The common parameter estimator is asymptotically normal, but centered at a nonzero bias term of order $1/n$. This is the asymptotic benchmark for the estimator-level correction (EC) discussed in Section (ref). Next result shows that the penalized likelihood correction ($\widehat{\theta}_{\mathrm{PL}}$) recenters the asymptotic distribution without inflating the variance.

Bias Correction for Average Partial Effects

We next study average partial effects. The same argument also applies to generic Average Structural Functions (ASF) that are smooth averages of dyad level link probabilities. We state a high-level condition that is sufficient for the PL-based bias correction mechanism.

Assumption (ref) imposes standard smoothness and boundedness on the generic dyadic structural functions. For the APEs studied in this paper, the derivative bounds follow from a simple count. The function $\Delta(\theta,\lambda)$ averages over $n(n-1)$ directed dyads. A common parameter enters every term in this average. A fixed effect for node $i$ enters only the $2(n-1)$ directed dyads involving that node. Thus, uniformly in a neighborhood of $(\theta_0,\lambda_0)$, we can verify that $ \|\partial_\theta\Delta\|=\mathcal{O}(1), $ $ \|\partial_\lambda\Delta\|=\mathcal{O}(n^{-1/2}), $ $ \|\partial_{\lambda\lambda'}\Delta\|=\mathcal{O}(n^{-1}), $ $ \|\partial_{\theta\theta'}\Delta\|=\mathcal{O}(1), $ and $ \|\partial_{\theta\lambda'}\Delta\|=\mathcal{O}(n^{-1/2}). $\footnote{ These rates follow from a simple count by the APE construction. For all $i$, we have $\max_i|\partial_{\lambda_i}\Delta | = \mathcal{O}(n^{-1})$. Hence, $\|\partial_\lambda\Delta\|\leq n^{1/2} \max_i|\partial_{\lambda_i}\Delta | =O(n^{-1/2})$. Moreover, because $\|\partial_{\lambda_i\lambda_i} \Delta \| = \mathcal{O}(n^{-1})$ and $\|\partial_{\lambda_i\lambda_j} \Delta \| = \mathcal{O}(n^{-2})$ for all $i\neq j$, it follows that $\|\partial_{\lambda\lambda'} \Delta\|_{\infty} = \|\partial_{\lambda\lambda'} \Delta\|_{1} =\mathcal{O}(n^{-1})$. This gives that $\|\partial_{\lambda\lambda'} \Delta\| \leq \sqrt{\|\partial_{\lambda\lambda'} \Delta\|_{\infty}\|\partial_{\lambda\lambda'} \Delta\|_{1}} = \mathcal{O}(n^{-1}).$ Finally, each entry of $\partial_{\theta\lambda'}\Delta$ is $O(n^{-1})$. Since $\dim(\theta)$ is fixed, $\|\partial_{\theta\lambda'}\Delta\| \leq n^{1/2}O(n^{-1})=O(n^{-1/2})$. } Taking more derivatives with respect to $\lambda$ does not enlarge the set of affected dyads. These structural limits guarantee that the Taylor remainders of the structural function shrink at appropriate rates.

Bias correction for coefficients does not automatically imply bias correction for APEs. As explained in Section (ref), the PL estimator removes the first-order bias coming from the common parameters and the fixed-effects point estimates. The only remaining first-order term is induced by fixed-effects uncertainty. Equation (ref) subtracts exactly this bias, and the next corollary establishes the resulting asymptotic distribution.

Corollary (ref) clarifies why the penalized likelihood approach is especially convenient for structural functions. After penalizing the likelihood, the only remaining bias term can be removed by subtracting the trace term above. It does not require third-order derivatives. In this sense, the PL approach keeps coefficient correction and APE correction within a single likelihood-based framework. From Corollary (ref), the bias-corrected APE estimator $\widetilde{\Delta}(\widehat\theta_{\mathrm{PL}},\widehat\lambda_{\mathrm{PL}})$ defined in equation (ref) is asymptotically centered at the true APE. Consequently, $ \sqrt{N}\bigl(\widetilde{\Delta}(\widehat\theta_{\mathrm{PL}},\widehat\lambda_{\mathrm{PL}})-\Delta\bigr)\overset{d}{\rightarrow}\mathcal N(0,V_{\Delta,\infty}). $

Monte Carlo Experiments

This section studies the finite sample performance of the PL estimator in designs tailored to the existence problem. We find that in sparse networks, where isolated nodes occur frequently, standard fixed effects methods often become unavailable, while the PL estimator remains computable and continues to deliver accurate bias correction. We benchmark the PL estimator against the uncorrected MLE, the estimator-level correction (EC), and conditional logit methods. We evaluate both the common parameters and the APEs.

Designs

Our data-generating processes adapt the simulation framework in graham2017econometric to accommodate asymmetric links and reciprocity. We consider three model classes: (i) the reciprocal directed model that is the leading case of the paper, (ii) the directed model without reciprocity, and (iii) the undirected model.

The directed covariate $X_{ij}$ is drawn from a discrete uniform distribution with possible values $\{0, 1\}$. We construct the symmetric mutual utility covariate as $Z_{ij} = \tilde{Z}_{i}\tilde{Z}_{j}$, where $\tilde{Z}_{i}$ is drawn from a discrete uniform distribution over $\{-1, 1\}$. The idiosyncratic error $\epsilon_{ij}$ follows a standard logistic distribution. To capture the correlation between observables and unobserved degree heterogeneity, we specify the sender and receiver fixed effects as $\alpha_i = \varrho_L \mathbbm{1}_{\{\tilde{Z}_i=-1\}} + \varrho_H \mathbbm{1}_{\{\tilde{Z}_i=1\}} + \vartheta_i^{\alpha}$ and $\gamma_i = \varrho_L \mathbbm{1}_{\{\tilde{Z}_i=-1\}} + \varrho_H \mathbbm{1}_{\{\tilde{Z}_i=1\}} + \vartheta_i^{\gamma}$ with $\varrho_L\leq \varrho_H$. Both $\vartheta_i^{\alpha}$ and $\vartheta_i^{\gamma}$ are drawn from the centered Beta distribution, $\text{Beta}(\varpi_0, \varpi_1) - \varpi_0/(\varpi_0 + \varpi_1)$. For the directed models, we simulate the network using a stochastic best-response dynamic, iterating 1,000 times for each dyad to ensure equilibrium convergence.\footnote{ The directed network graph with reciprocity is simulated via an iterative updating scheme based on $ g_{ij} = \mathbbm{1}{\{ X_{ij}\beta_0 + g_{ji}Z_{ij}\rho_0 +\alpha_i + \gamma_j \geq \epsilon_{ij} \}}. $ We begin by initializing $g_{ji}=0$. The directed network graph without reciprocity uses the same iterative scheme based on $ g_{ij} \mathbbm{1}{\{Z_{ij}\beta_0 +\alpha_i + \gamma_j \geq \epsilon_{ij} \}}. $ The undirected network graph is generated directly by $ g_{ij} = g_{ji}= \mathbbm{1}{\{ Z_{ij}\rho_0 +\alpha_i + \gamma_j \geq \epsilon_{ij} \}}. $ } We set the true common parameters to $\beta_0 = \rho_0 = 1$. For each design, we report results for $n \in \{100, 200\}$ based on 1,000 Monte Carlo replications.

Table (ref) summarizes six designs. In Designs A.3 and B.3, the fixed effects MLE succeeds in only 0.2%–0.3% of replications in the reciprocal model and never succeeds in the nested models. In these same designs, the PL estimator is available in every replication and maintains near-nominal coverage for both coefficients and APEs.

Monte Carlo Results

Tables (ref) and (ref) contain the Monte Carlo evidence for the reciprocal directed model. In the moderate dense network designs (A.1, A.2, B.1, B.2), the uncorrected MLE and estimator based correction (EC) are available throughout. The bias of the uncorrected MLE is visible in both $\widehat{\beta}$ and $\widehat{\rho}$, and EC removes most of that bias. The penalized-likelihood estimator (PL) performs excellent in these designs.

In the sparse designs (A.3, B.3), existence becomes the primary issue. The uncorrected MLE is undefined and therefore reported as “n.a.” throughout the tables. The EC method is unavailable for the same reason, since it is built on the underlying MLE. By contrast, the PL estimator is available in every replication. This is the central finite-sample implication of Theorem (ref).

The PL estimator remains accurate once sparsity becomes severe. In A.3 and B.3, the empirical coverage of PL estimates stays close to the nominal level for the common parameters. Panels C and D show the same pattern for the APEs. In sparse networks, PL is useful not only because the estimator exists. It also supports inference on coefficients and APEs when MLE-based corrections are infeasible. The same conclusion is even clearer when $n$ increased to 200 (Table (ref)).

Nested Benchmark Models. The supplemental appendix (ref) reports parallel results for the directed model without reciprocity and for the undirected model (Tables (ref) and (ref)). For the directed model without reciprocity, we also report the quadruple-logit (QL) estimates of jochmans2018semiparametric. For the undirected model, we report the tetrad-logit (TL) benchmark of graham2017econometric. The main pattern, however, is unchanged. As the network becomes sparse, the same existence issue arises, while the PL method remains implementable and continues to perform well for both common parameters and APEs bias corrections.

\begingroup \fontsize{8.5pt}{9.0pt}\selectfont

\endgroup

\begingroup \fontsize{8.5pt}{9.0pt}\selectfont

\endgroup

Empirical Application: The Global Textiles Trade Network

We apply our penalized likelihood estimator to model the formation of the 1986 global textiles trade network. Using data from the CEPII TradeProd database, we observe directed export links among 133 countries, alongside standard gravity covariates from the trade literature helpman2008estimating, including geographic, institutional, and cultural proximity measures.\footnote{The covariates include log distance, indicators for a shared land border, island status, landlocked status, common legal origin, common language, currency union, common free trade agreement, and religious proximity. Variable definitions and descriptive statistics are reported in Appendix (ref) of the supplemental material.} A directed link $g_{ij}=1$ indicates that country $i$ exports textiles to country $j$.

The application is useful for two reasons. First, international trade is naturally directed and highly reciprocal. To consistently estimate the effects of bilateral frictions, researchers may control for unobserved degree heterogeneity via exporter and importer fixed effects.

Second, the network contains boundary-degree countries even though it is not globally sparse. The standard MLE cannot accommodate such nodes and therefore requires iterative trimming. In this application, it also removes central economies once the initial zero in-degree countries are dropped. The PL estimator avoids this sample attrition and estimates the model on the full observed network.

Network Topology and Existence Failure

Figure (ref) summarizes the topology of the textiles trade network. The network contains 133 countries and has density 39.6 percent. It is therefore not globally sparse. Reciprocity is high. Mutual dyads account for 30.1 percent of all country pairs, almost twice the rate implied by an independent network formation benchmark. This high reciprocity rate validates the inclusion of the mutual utility component in our specification.

Panels A and B show substantial degree heterogeneity. 11 countries (e.g., Albania, Bulgaria, and Vietnam) import no textiles from any partner in the network ($b_i=0$). They export textiles to at least one destination but import textiles from no country in the observed network. Standard estimation protocols require dropping these 11 countries from the sample.

The trimming problem is iterative, and triggers a cascading selection problem. Removing the 11 importing countries alters the network sample for the remaining exporters. Conditional on this restricted 122 country sample, five major economies—France, Germany, Italy, Japan, and the UK—become full out-degree nodes, meaning they export to every remaining country in the data. Consequently, the standard MLE forces the deletion of these 5 global trade powerhouses.

The final MLE sample shrinks to 117 countries. Estimating a global trade model without France, Germany, Italy, Japan, and the UK severely distorts the economic reality of the network and removes the most informative nodes. In contrast, the PL estimator retains the full 133 countries network and avoid this structural sample selection bias.

Estimation Results

Table (ref) reports the coefficient estimates and average partial effects for the reciprocal directed model. In Panel A, the directed utility parameters follow established gravity patterns. Distance strictly reduces the probability of a unilateral trade link (PL APE: $-0.098$), whereas sharing a common language (PL APE: $0.032$) and participating in a free trade agreement (PL APE: $0.222$) increase it.

Panel B reports the mutual utility parameters, and demonstrates that bilateral variables affect unilateral and mutual trade through distinct channels. Sharing a land border strongly drives mutual trade (PL APE: $0.078$) but exhibits a negative direct effect (PL APE: $-0.076$). Conversely, island status reduces the probability of reciprocal trade (PL APE: $-0.067$) despite increasing the likelihood of unilateral links. These diverging partial effects confirm that explicitly modeling reciprocity is critical for recovering the correct mechanisms of network formation.

The differences between MLE and PL are economically meaningful. Relative to the MLE, the PL estimates are systematically less pronounced, while the standard errors remain stable. For example, the directed utility APE for a land border falls in magnitude from $-0.103$ to $-0.076$, and the mutual utility APE falls from $0.108$ to $0.078$. This pattern aligns precisely with our asymptotic prediction: the incidental parameter problem primarily shifts the center of the distribution rather than scaling the variance. Failure to correct this bias, combined with the mechanical exclusion of core trading counties, may lead researchers to overstate the effects of trade frictions.

Concluding Remarks

This paper studies fixed-effects estimation in dyadic network formation models with degree heterogeneity when the likelihood is affected by boundary-degree nodes and incidental parameter bias. In empirical networks, the fixed effects MLE frequently fails to exist due to boundary-degree nodes. Even when it exists, the incidental parameter problem biases both common parameters and average partial effects. We address these problems with a penalized likelihood estimator for a reciprocal directed network formation model that nests the existing undirected and nonreciprocal directed specifications. The penalty keeps the estimator well defined without trimming isolated or boundary-degree nodes and yields bias corrections for both coefficients and APEs.

The asymptotic analysis does not assume compactness of the fixed-effects parameter space. We allow the true sender and receiver effects to grow at a logarithmic rate with network size, so some link probabilities may shrink as the network expands. The reciprocal model also requires a different fixed-effects Hessian approximation because sender and receiver effects for the same node are locally coupled. We derive an explicit inverse-Hessian approximation that accounts for this coupling and is sufficient for the coefficient and APE expansions.

In the empirical application, the textiles network has moderate density, yet local zero- or full-degree countries still force the standard MLE onto a selected subnetwork. The PL estimator instead uses the full observed network and yields bias-corrected coefficients and APEs.

Several extensions are natural. One is likelihood-based inference built directly on the penalized objective. Another is to adapt the same logic to richer forms of dyadic dependence beyond reciprocity. These questions are left for future work.