(Debiased) Inference for Fixed Effects Estimators with Three-Dimensional Panel and Network Data

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

92,905 characters · 12 sections · 99 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\thispagestyle{empty}

abstractInference for fixed effects estimators is often unreliable due to Nickell- and incidental parameter biases. While these issues are well understood for classical two-dimensional panels, little is known about three-dimensional panel structures (e.g., sender $\times$ receiver $\times$ time). We develop inferential theory for a broad class of linear and nonlinear fixed effects M-estimators in this setting, covering bipartite, directed, and undirected network panel data, multiple specifications of additively separable unobserved effects, and both strictly exogenous and predetermined regressors. Our analysis reveals fundamentally different asymptotic properties compared to two-dimensional panels. In particular, we find a sharp dichotomy across specifications: (i) when unobserved effects vary along a single panel dimension, the estimator is asymptotically unbiased; (ii) when they vary along two panel dimensions, the estimator suffers from a severe inference problem characterized by a degenerate asymptotic distribution. We resolve the latter by deriving explicit bias formulas and proposing analytically debiased estimators with nondegenerate, correctly centered asymptotic distributions. An empirical application studies dynamic network formation in a directed panel of bilateral trade relationships.\\[1em] JEL Classification: C13, C23\\ Keywords: panel data, network data, dynamic model, multi-way fixed effects, incidental parameter problem, asymptotic bias correction.

\onehalfspacing

\setcounter{page}{1}

Introduction

The inconsistency of fixed effects estimators for linear and nonlinear panel models, first noted by ns1948 and n1981, remains a highly studied topic in panel data econometrics. Although empirical researchers increasingly use more granular, multi-dimensional panel data with layers of variation beyond the classical cross-sectional and time dimensions, available solutions to these inference problems typically focus on two-dimensional panels.

This paper contributes to closing this gap by developing an inferential theory for fixed effects M-estimators in a generic class of (non)linear models with three-dimensional panels, i.e., panels with two cross-sectional dimensions and one time dimension, covering bipartite, directed, and undirected network panels. An example is data on bilateral network activities observed over time, widely used in international trade research to study flows between each $i = 1, \dots, N_{1}$ exporting country and each $j = 1, \dots, N_{1}$ importing country for each $t = 1, \dots, T$ year. This multidimensionality allows researchers to control for multiple, richer sources of unobserved heterogeneity varying along one or two panel dimensions. We call heterogeneity arising from the interaction of two panel indices interacted and the one-dimensional counterpart non-interacted. For example, empirical analyses in international trade typically control for interacted unobserved heterogeneity of the form exporter-by-year, importer-by-year, and importer-by-exporter, using the three-way fixed effects specification $\alpha_{it} + \gamma_{jt} + \rho_{ij}$.\footnote{Controlling for this type of unobserved heterogeneity is, among others, recommended in hm2014.}

We make the following main contributions. First, we develop the first generic asymptotic framework for inference on fixed effects M-estimators in three-dimensional panel models in which all dimensions grow jointly. This framework covers different data structures (bipartite, directed, and undirected network panel data), several types of fixed effects, and accommodates both strictly exogenous and predetermined regressors. In our main analysis, we derive the asymptotic properties of fixed effects estimators for models with three sets of unobserved heterogeneity in interacted and non-interacted form.\footnote{In a bipartite panel, the fixed effects enter the linear index as $\alpha_{it} + \gamma_{jt} + \rho_{ij}$ in the interacted case and as $\alpha_{i} + \gamma_{j} + \rho_{t}$ in the non-interacted case.} These two specifications are the most general cases and nest a wide class of one- and two-way fixed effects models, allowing us to characterize the asymptotic properties of the corresponding estimators across all such settings. For cases in which inference would otherwise be invalid, we develop appropriate debiased estimators.\footnote{The estimators are computationally tractable and are implemented in our R-package \href{https://cran.r-project.org/web/packages/alpaca/index.html}{alpaca}.} Second, we develop new analytical tools that address the specific technical challenges of the multi-dimensional panel setting but which also apply to two-dimensional settings. Treating unobserved effects as incidental parameters creates a non-standard inference problem because the dimension of the parameter vector grows with all three panel dimensions. Moreover, in multi-way interacted specifications, shared indices of two unobserved effects induce a high-dimensional incidental parameter Hessian with a nontrivial sparsity pattern. These features render standard approaches from two-dimensional panels inadequate and motivate new analytical tools, which may also prove useful for future work on multi-dimensional panels. Specifically, we extend the asymptotic expansions of fw2016 from third to fourth order, derive a new bound on the $p$-th moment of the spectral norm of a random matrix whose rows are independent mean-zero weakly dependent processes, and develop a generic strategy for approximating the inverse of the high-dimensional expected incidental parameter Hessian in sparse settings. These challenges may explain the substantial gap between the literature on fixed effects estimators for two-dimensional versus multi-dimensional panels. Third, our findings clarify that standard results for two-dimensional panel data models do not simply carry over to three-dimensional settings. Finally, our insights raise serious concerns about the increasingly common empirical practice of using even more granular panel data, for example with four dimensions.

Our first main finding concerns models with three-way interacted unobserved heterogeneity. We identify an unusually severe inference problem caused by a stark imbalance between the order of the leading bias and the order of the standard deviation. Since the former exceeds the latter, the asymptotic distribution degenerates and the inference problem grows worse asymptotically.\footnote{For bipartite panel data, the order of the leading bias is $1 / N_{2} + 1 / N_{1} + 1 / T$ and the order of the standard deviation is $1 / \sqrt{N_{1} N_{2} T}$, where $N_{1}$ and $N_{2}$ denote the lengths of the two cross-sections and $T$ the number of time periods.} This contrasts with the typical result in the two-dimensional literature, where the orders of bias and standard deviation can be balanced by choosing appropriate rates, yielding a nondegenerate asymptotic distribution with a distorted first moment. Although the inference problem is more severe here, the bias can be estimated at a sufficient rate, and we develop an analytically debiased estimator with a nondegenerate asymptotic distribution centered at zero.

Our second set of results concerns models with three-way non-interacted unobserved heterogeneity. This is the opposite extreme: the fixed effects estimator is asymptotically free of any incidental parameter bias, even though the dimension of the fixed effects grows with the sample size. This is because the leading bias of the fixed effects estimator always shrinks faster than the standard deviation.\footnote{For bipartite panel data, the order of the bias is $1 / (N_{1} T) + 1 / (N_{2} T) + 1 / (N_{1} N_{2})$ and the order of the standard deviation is $1 / \sqrt{N_{1} N_{2} T}$.}

Our simulation experiments confirm the severity of the inference problem for the interacted specification: confidence intervals around the uncorrected estimator almost never cover the true parameter. Our proposed bias correction substantially improves inferential accuracy. For the non-interacted specification, we confirm empirically that the fixed effects estimator is asymptotically unbiased.

In our empirical application, we apply our debiased estimator to study the dynamic formation of export and import relationships between countries, using a probit model with three-way interacted unobserved heterogeneity. We adapt the dynamic network formation model of g2016 to a directed network setting. We find strong evidence for state dependence: any past trade relationship, regardless of direction, increases the probability of an export link today. We also find evidence for two triadic transitivity channels, indirect trade linkages through intermediary countries and shared upstream suppliers, each of which raises the probability of a direct bilateral trade link forming.

\noindentRelated literature. Our work contributes to the large-$T$ literature on bias corrections for the inconsistency of fixed effects estimators identified by n1981 and ns1948.\footnote{Another important strand focuses on fixed-$T$ asymptotics, where estimators typically eliminate unobserved effects by differencing or conditioning on sufficient statistics. Examples for binary logit models include r1960, a1970, c1980, hk2000, and hw2025.} The largest part of this literature has focused on classical panel data models with individual effects. For linear models, see among others hk2002 and dj2015; for nonlinear models, see among others l2002, w2002, s2003, hn2004, c2007, ab2009, bh2009, f2009, hk2011, dj2015, ks2016, p2019, sst2021, hj2022, and schumann2023. These solutions differ in their assumptions, techniques, and types of corrections. We refer the reader to ah2007 and fw2018 for comprehensive reviews. For linear models, hm2006 show that additional time effects do not introduce a further inference problem, so standard Nickell-bias corrections developed for individual effects carry over. In nonlinear models, time effects introduce an additional incidental parameter bias. fw2016 and yjfl2018 advance this literature by developing asymptotic theory and bias corrections for nonlinear models with two types of fixed effects. Whereas yjfl2018 focus on logit models with sender and receiver fixed effects for directed network data, fw2016 cover nonlinear models with individual and time effects, accommodating time dependence and predetermined regressors.\footnote{The analysis of fw2016 also applies to directed networks where the second panel dimension is not time but another cross-section, such as countries or industries. Their theory therefore nests the structural parameter results of yjfl2018.} Using a likelihood correction approach, jo2019 propose an alternative to the ex-post distribution correction of fw2016. Building on fw2016, d2019 extends the framework to directed network models with sender and receiver fixed effects, developing specific inference procedures and specification tests, and h2026 introduces an improved jackknife bias correction for network models. For undirected network data, g2017 develops a bias-corrected estimator for network formation with degree heterogeneity.\footnote{Another strand develops conditional logit estimators in the spirit of r1960 and c1980; see for example g2017 for undirected network data, or c2017 and j2018 for directed network data.} Despite the growing use of multi-dimensional panels in empirical work, wz2021 is the only other paper in this literature that studies fixed effects estimators for three-dimensional panel data. They derive properties of the three-way fixed effects Pseudo-Poisson estimator for the gravity model. Unlike us, wz2021 assume strict exogeneity of regressors and exploit a property specific to the Poisson estimator that reduces the problem essentially to a two-way fixed effects analysis, allowing them to rely substantially on fw2016. As a consequence, they only require the two cross-sectional dimensions to grow, while we require all three panel dimensions to grow. h2026_jackknife develops a jackknife $t$-statistic for a broad class of fixed effects models, including multi-dimensional panels; applying it in our setting requires knowledge of the asymptotic bias structure, which we provide in this paper.

A growing interest in methods for three-dimensional panel models has also emerged in other strands of the econometrics literature. For example, g2016 develops a conditional fixed effects logit estimator for dynamic network formation in a directed network panel model with pair-specific unobserved heterogeneity. mp2025 recently introduce a conditional fixed effects logit estimator for a triadic network model with three-way interacted unobserved heterogeneity. yh2023 provide an estimator for a gravity model with three multiplicative unobserved effects, extending the GMM strategy of j2017. f2022 and jls2025 propose interactive fixed effects estimators for linear panel models with predetermined regressors.

\noindentOutline. Section (ref) introduces the data structures, model, and fixed effects estimators. Section (ref) presents the asymptotic theory. Section (ref) reports simulation results. Section (ref) presents the empirical application. Section (ref) concludes.

\noindentNotation. Throughout the paper, $\mathbb{P}\left(\cdot\right)$ and $\mathbb{E}\left[ \cdot \right]$ denote probability and expectation. A superscript “0” on a parameter denotes its true population value. We write $\text{a.\,s.}$ and $\text{wpa1}$ for almost surely and with probability approaching one, respectively. We use $o_{P}(1)$ to denote a sequence of random variables that converges in probability to zero, and $\mathcal{O}_{P}(1)$ to denote a sequence that is bounded in probability. We use $\xrightarrow{d}$ and $\xrightarrow{p}$ to denote convergence in distribution and in probability, respectively. Unless otherwise stated, all stochastic statements are understood as almost sure statements conditional on $\Phi$, the sigma-algebra generated by the unobserved effects and initial conditions.

Data Structures, Models, and Estimators

Data Structures

We observe three-dimensional panel data $\{(y_{ijt}, x_{ijt}) \colon (i, j) \in \mathcal{D}_{s}, \, t \in \{1, \ldots, T\}\}$, where $y_{ijt}$ is an outcome variable, $x_{ijt}$ is a $K$-dimensional vector of explanatory variables, and $\mathcal{D}_{s}$ is the set of observed pairs $(i, j)$ corresponding to one of the three panel structures indexed by $s \in \{1, 2, 3\}$. Let $\mathcal{I} \coloneqq \{1, \ldots, N_{1}\}$ and $\mathcal{J} \coloneqq \{1, \ldots, N_{2}\}$ denote the index sets for the two cross-sectional dimensions.

In bipartite panel data ($s = 1$), the indices $i$ and $j$ refer to two distinct samples of cross-sectional units, and the set of observed pairs is $\mathcal{D}_{1} = \{(i, j) \colon i \in \mathcal{I}, j \in \mathcal{J}\}$. In network panel data ($s \in \{2, 3\}$), both indices refer to the same set of $N_{1}$ agents, and self-ties are excluded ($i \neq j$). Network panels can be directed or undirected. In directed network panels ($s = 2$), the relationship between agents is asymmetric: $(y_{ijt}, x_{ijt})$ need not equal $(y_{jit}, x_{jit})$. The set of observed pairs is $\mathcal{D}_{2} = \{(i, j) \colon i \in \mathcal{I}, j \in \mathcal{I}, i \neq j\}$. In undirected network panels ($s = 3$), the relationship is symmetric: $(y_{ijt}, x_{ijt}) = (y_{jit}, x_{jit})$. To avoid duplicates, the observed set is restricted to $\mathcal{D}_{3} = \{(i, j) \colon i \in \mathcal{I}, j \in \mathcal{I}, i < j\}$. The total number of observations is $n_{s} \coloneqq \lvert \mathcal{D}_{s} \rvert \, T$, giving sample sizes $n_{1} = N_{1} N_{2} T$, $n_{2} = N_{1} (N_{1} - 1) T$, and $n_{3} = N_{1} (N_{1} - 1) T / 2$.

example[Bipartite Panel Data] Firm-level panel data often contain information beyond the standard firm and time identifiers that can be used to form a bipartite panel. Examples include data on products or industries, export destinations, and the locations of subsidiaries or investments.

example[Network Panel Data] Data on countries observed over time can often be used to form a network panel. Country pairs may represent origin and destination in bilateral trade or migration data, initiator and signatory in data on bilateral agreements, or aggressor and target in data on (military) conflicts.

remark[Tripartite Data and Unbalancedness] The time dimension can also be replaced by a third cross-sectional dimension. The exclusion of self-ties is one source of unbalancedness in network panels. All three data structures can be further unbalanced along the cross-sectional and time dimensions (see fw2018). For notational simplicity, we treat self-tie exclusion as the only source of unbalancedness and otherwise work with balanced data.

Model

We consider the following generic semiparametric unobserved effects model for $(i, j) \in \mathcal{D}_{s}$, $s \in \{1, 2, 3\}$, and $t \in \{1, \dots, T\}$,

equation[equation omitted — 157 chars of source]

where $\mathcal{X}_{ij}^{t^{\prime}} \coloneqq \sigma(\{x_{ijt^{\prime\prime}} \colon t^{\prime\prime} \in \{1, \ldots, t^{\prime}\}\})$, $g(\cdot)$ is a known link function, and $\pi_{ijt}(\beta, \mu) \coloneqq x_{ijt}^{\prime} \beta + \mu$ is the linear index. Here, $\beta$ is a $K$-dimensional parameter vector and $\mu$ is a scalar. Let $\phi$ be an $L$-dimensional (block) vector stacking $\mathcal{M} \in \{1, 2, 3\}$ vectors of unobserved effects. The function $\mu_{ijt}(\phi)$ sums a specific subset of up to $\mathcal{M}$ component vectors of $\phi$, where the subset may vary across observations. We impose no restrictions on the relationship between the unobserved effects and the explanatory variables, and we do not assume a parametric distribution for the unobserved effects.

Table (ref) presents 17 identifiable linear combinations through which unobserved heterogeneity can enter the linear index $\pi_{ijt}(\beta, \mu_{ijt}(\phi))$.

table[table omitted — 1,446 chars of source]

The specifications in Table (ref) differ along two dimensions. First, they differ in $\mathcal{M}$, the number of unobserved-effect vectors stacked in $\phi$. Second, each component vector may capture heterogeneity varying along one panel dimension (non-interacted) or two (interacted). For undirected network panels, the symmetry condition $(y_{ijt}, x_{ijt}) = (y_{jit}, x_{jit})$ imposes additional restrictions on the unobserved effects. For instance, in a directed network panel one might specify $\mu_{ijt}(\phi) = \alpha_{it} + \gamma_{jt}$, with $\alpha_{it}$ and $\gamma_{jt}$ as distinct effects. In an undirected panel, symmetry requires $\gamma_{jt} = \alpha_{jt}$, so that $\mu_{ijt}(\phi) = \alpha_{it} + \alpha_{jt}$.

The following three examples illustrate model (ref) across different specifications and data structures.

example[Linear Model for Bipartite Panel Data] Let $y$ be a dependent variable with a conditional mean: \begin{equation*} \mathbb{E}\left[ y_{ijt} \mid \mathcal{X}_{ij}^{t}, \beta, \phi \right] = x_{ijt}^{\prime} \beta + \mu_{ijt}(\phi) \, , \quad \left( (i,j) \in \mathcal{D}_1, t \in \{1, \dots, T\}\right) \, . \end{equation*} For example, ss2015technology use a country-industry-time panel to analyze how technological characteristics interact with business cycle contractions to affect industry growth. They employ the interacted three-way specification (3.b) in Table (ref).

example[Exponential Model for Directed Network Panel Data] Let $y$ be a positive outcome variable with a conditional mean: \begin{equation*} \mathbb{E}\left[ y_{ijt} \mid \mathcal{X}_{ij}^{t}, \beta, \phi \right] = \exp\big(x_{ijt}^{\prime} \beta + \mu_{ijt}(\phi)\big) \, , \quad \left( (i,j) \in \mathcal{D}_2, t \in \{1, \dots, T\}\right) \, . \end{equation*} In international trade research, this specification is commonly used to estimate the impact of trade policy variables, such as joint WTO membership, on bilateral trade flows. Researchers using directed network panel data typically model the relationship between exporter $i$ and importer $j$ at time $t$ with an exponential conditional mean, as recommended by st2006 and lsy2025, among others. Standard practice is to account for exporter-by-time, importer-by-time, and exporter-by-importer unobserved heterogeneity, as in specification (3.b) of Table (ref).

example[Binary Response Model for Undirected Network Panel Data] Let $y$ be a binary dependent variable with a conditional mean: \begin{equation*} \mathbb{E}\left[ y_{ijt} \mid \mathcal{X}_{ij}^{t}, \beta, \phi \right] = F_{\epsilon}(x_{ijt}^{\prime} \beta + \mu_{ijt}(\phi)) \, , \quad \left( (i,j) \in \mathcal{D}_3, t \in \{1, \dots, T\}\right) \, , \end{equation*} where $F_{\epsilon}(\cdot)$ is a suitable cumulative distribution function, such as the standard normal CDF in the probit case. This specification can be used to study the formation of regional trade agreements between country pairs over time, in the spirit of bb2007. Since both the outcome and typical regressors, such as country-pair similarity measures, are symmetric in $i$ and $j$, the data naturally form an undirected network panel. The richest unobserved heterogeneity structure available for this model is the symmetric counterpart of specification (3.b), namely $\alpha_{it} + \alpha_{jt} + \rho_{ij}$.

Fixed Effects Estimators

We follow a fixed effects approach and treat the unobserved effects $\phi$ as nuisance parameters to be estimated jointly with $\beta$.

For each $s \in \{1, 2, 3\}$, we assume the true parameter values are identified by the population problem:

equation[equation omitted — 272 chars of source]

where the objective function,

equation[equation omitted — 285 chars of source]

consists of a criterion function $\psi_{ijt}(\pi) \coloneqq \psi(y_{ijt}, \pi)$ and a penalty term that ensures unique identification of $\phi$ by imposing suitable linear constraints. For instance, $\psi_{ijt}(\pi) = (y_{ijt} - \pi)^2$ corresponds to OLS estimation (as in Example (ref)). The objective function is normalized by a model-specific scaling factor $w_{n, s}$ that grows with the sample size $n_{s}$ at a suitable rate. We restrict attention to objective functions that are locally convex near the true parameter values and sufficiently smooth in all parameters. This covers a wide range of estimators, including many popular (pseudo-)maximum likelihood and (non)linear least squares estimators, but excludes those with non-smooth criterion functions, such as those used in quantile regression.

We estimate $\beta_{s}^{0}$ and $\phi_{s}^{0}$ by minimizing the sample analogue of (ref). For each $s \in \{1, 2, 3\}$, the M-estimator is:

equation[equation omitted — 258 chars of source]

remark[Penalty Term] To understand the identification problem and the role of the penalty term in (ref), consider a bipartite panel with three-way interacted unobserved heterogeneity (specification (3.b) in Table (ref)), so that $\pi_{ijt}(\beta, \mu_{ijt}(\phi)) = x_{ijt}^{\prime} \beta + \alpha_{it} + \gamma_{jt} + \rho_{ij}$. Because the incidental parameters enter the linear index additively, the objective function is invariant to certain parameter transformations. For example, adding a constant $c_{t}$ to all $\alpha_{it}$ while subtracting it from all $\gamma_{jt}$ leaves the linear index unchanged. Likewise, subtracting a constant $c_{i}$ from all $\alpha_{it}$ while adding it to all $\rho_{ij}$, or adding a constant $c_{j}$ to all $\gamma_{jt}$ while subtracting it from all $\rho_{ij}$, leaves the linear index unchanged. To resolve this identification problem, we impose a system of linear constraints on the incidental parameters. In particular, we require: $\sum_{i = 1}^{N_{1}} \alpha_{it} = \sum_{j = 1}^{N_{2}} \gamma_{jt}$ for all $t \in \{1, \ldots, T\}$; $\sum_{t = 1}^{T} \alpha_{it} = \sum_{j = 1}^{N_{2}} \rho_{ij}$ for all $i \in \{1, \ldots, N_{1}\}$; and $\sum_{t = 1}^{T} \gamma_{jt} = \sum_{i = 1}^{N_{1}} \rho_{ij}$ for all $j \in \{1, \ldots, N_{2}\}$. Together, these $T + N_{1} + N_{2}$ restrictions are encoded as $V^{\prime} \phi = \mathbf{0}_{T + N_{1} + N_{2}}$, where $V$ is an $L \times (T + N_{1} + N_{2})$ matrix defined as \begin{equation*} V = \begin{pmatrix} \iota_{N_{1}} \otimes I_{T} & I_{N_{1}} \otimes \iota_{T} & 0 \\ - \iota_{N_{2}} \otimes I_{T} & 0 & I_{N_{2}} \otimes \iota_{T} \\ 0 & - I_{N_{1}} \otimes \iota_{N_{2}} & - \iota_{N_{1}} \otimes I_{N_{2}} \end{pmatrix} \, . \end{equation*}

remark[Undirected Network Panels] An undirected network panel model, for example with linear index $\pi_{ijt}(\beta, \mu_{ijt}(\phi)) = x_{ijt}^{\prime} \beta + \alpha_{it} + \alpha_{jt}$, is typically estimated on a sample of unique pairs $(i, j)$. Alternatively, one can use a doubled sample that includes both $(i, j)$ and $(j, i)$ and estimate the more general specification $\pi_{ijt}(\beta, \mu_{ijt}(\phi)) = x_{ijt}^{\prime} \beta + \alpha_{it} + \gamma_{jt}$. Both approaches yield identical estimates of $\beta^{0}$.\footnote{A similar equivalence was recently stated by h2026 and lsw2025 in the context of network data (without time dimension).} We exploit this equivalence to derive our asymptotic results for this data structure.

remark[Computation] In practice, joint estimation of $\beta$ and $\phi$ may involve a very large number of parameters, making (ref) a high-dimensional optimization problem. Standard software routines that rely on dummy variables are therefore often computationally infeasible, even for moderately sized panels. We recommend algorithms such as those of gp2010, g2013, b2018, and s2018, which are designed for high-dimensional problems of this type and can handle unbalanced data.\footnote{When (ref) reduces to OLS, i.e., $g(\cdot)$ is the identity function and $\psi(\cdot)$ is the quadratic loss function, closed-form within-group transformations exist in two specific cases: (i) one-way fixed effects structures; and (ii) two- or three-way fixed effects structures for balanced panels.} We provide computationally efficient implementations of the (debiased) estimators studied in this paper in the R package \href{https://cran.r-project.org/web/packages/alpaca/index.html}{alpaca}.

Asymptotic Theory

Assumptions

Before stating our assumptions, we introduce additional notation. For each $s \in \{1, 2, 3\}$, let $X$ denote the $n_{s} \times K$ matrix of explanatory variables, where $x_{ijt}^{\prime}$ is its $ijt$-th row. Let $(d^{r} \psi(\beta, \phi))_{ijt}$ denote the $r$-th derivative of the criterion function with respect to the linear index, $\partial^{r} \psi_{ijt}(\pi_{ijt}(\beta, \mu_{ijt}(\phi))) / \left(\partial \pi_{ijt}\right)^{r}$. When evaluated at the true parameter values, we suppress the arguments and write $(d^{r} \psi)_{ijt}$. Let $(d_{\mathcal{X}}^{r} \psi(\beta, \phi))_{ijt} \coloneqq \mathbb{E}\big[(d^{r} \psi(\beta, \phi))_{ijt} \mid \mathcal{X}_{ij}^{t}\big]$ denote the corresponding conditional expectation.

We further define population weighted least-squares projections. For each $s \in \{1, 2, 3\}$ and $k \in \{1, \ldots, K\}$, let

align[align omitted — 533 chars of source]

where $x_{ijt, k}$ denotes the $k$-th element of $x_{ijt}$. The resulting $n_{s} \times K$ matrix of fitted values is denoted by $\mathfrak{X}$, with its $ijt$-th row defined as

equation[equation omitted — 155 chars of source]

Thus, the fitted values in (ref) are the weighted least-squares projection of $\mathbb{E}\left[ x_{ijt} (d_{\mathcal{X}}^{2} \psi)_{ijt} \right] / \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$ onto the subspace spanned by the incidental parameters. Under OLS, where $g(\pi) = \pi$, $\psi_{ijt}(\pi) = (y_{ijt} - \pi)^2 / 2$, and $(d_{\mathcal{X}}^{2} \psi)_{ijt} = 1$, (ref) reduces to the familiar population between-transformation. For example, for a bipartite panel ($s = 1$) with three-way interacted unobserved heterogeneity (specification (3.b) in Table (ref)), the corresponding between-transformation follows from bmw2018:

align*[align* omitted — 859 chars of source]

We define the population residuals as $\ddot{X} \coloneqq X - \mathfrak{X}$.

We make the following assumptions.

assumption[Sampling and regularity conditions] Let $z_{ijt} = (y_{ijt}, x_{ijt})$, $\nu > \check{\nu} > 0$, and $\varphi > 10 (20 + \nu) / \nu > 10$. Furthermore, for each $s \in \{1, 2, 3\}$, let $\varepsilon > 0$ and let $\Theta_{s}^{0}(\varepsilon)$ be a subset of $\operatorname{\mathbb{R}}^{K + 1}$ that contains an $\varepsilon$-neighborhood of $(\beta_{s}^{0}, \mu_{ijt}(\phi_{s}^{0}))$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. \begin{enumerate}[i)] • Asymptotics: We consider joint limits in which all panel dimensions diverge proportionally. For $s = 1$: $N_{1}, N_{2}, T \rightarrow \infty$ with $N_{1} / T \rightarrow \tau_{1} \in (0, \infty)$ and $N_{2} / T \rightarrow \tau_{2} \in (0, \infty)$. For $s \in \{2, 3\}$: $N_{1}, T \rightarrow \infty$ with $N_{1} / T \rightarrow \tau_{1} \in (0, \infty)$. • Sampling: For each $s \in \{1, 2, 3\}$, conditional on $\Phi$, $\{ \{z_{ijt}\}_{t = 1}^{T} \colon (i, j) \in \mathcal{D}_{s}\}$ is independent across $(i, j)$, and, for each $(i, j)$, $\{z_{ijt}\}_{t = 1}^{T}$ is $\alpha$-mixing with mixing coefficients satisfying $\sup_{ij} a_{ij}(q) = \mathcal{O}(q^{- \varphi})$ a.\,s. as $q \rightarrow \infty$, where $\mathcal{A}_{ij}^{t}$ is the sigma-algebra generated by $(z_{ijt}, z_{ij(t - 1)}, \ldots)$, $\mathcal{B}_{ij}^{t}$ is the sigma-algebra generated by $(z_{ijt}, z_{ij(t + 1)}, \ldots)$, and \begin{equation*} a_{ij}(q) \coloneqq \sup_{t} \sup_{A \in \mathcal{A}_{ij}^{t}, B \in \mathcal{B}_{ij}^{t + q}} \left\lvert \mathbb{P}\left(A \cap B\right) - \mathbb{P}\left(A\right) \mathbb{P}\left(B\right) \right\rvert \quad a.\,s. \end{equation*} • Model: For each $s \in \{1, 2, 3\}$, we assume that for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$, \begin{equation*} \mathbb{E}\left[ y_{ijt} \mid \mathcal{X}_{ij}^{t} \right] = g(\pi_{ijt}(\beta_{s}^{0}, \mu_{ijt}(\phi_{s}^{0}))) \, . \end{equation*} The unobserved effects $\phi_{s}^{0}$ are normalized to $V V^{\prime} \phi_{s}^{0} = \mathbf{0}_{L}$. • Smoothness and moments: For each $s \in \{1, 2, 3\}$, we assume that $(\beta, \mu) \rightarrow \psi_{ijt}(\pi_{ijt}(\beta, \mu))$ is five times continuously differentiable over $\Theta_{s}^{0}(\varepsilon)$ a.\,s. Each element of $x_{ijt}$, and the partial derivatives of $\psi_{ijt}(\pi_{ijt}(\beta, \mu))$ with respect to the elements of $(\beta, \mu)$ up to fifth order, are bounded in absolute value uniformly over $(\beta, \mu) \in \Theta_{s}^{0}(\varepsilon)$ by a function $\Psi(z_{ijt}) > 0$ a.\,s. In addition, $\sup_{ijt} \mathbb{E}\left[ \Psi(z_{ijt})^{20 + \nu} \right]$ is a.\,s.\ uniformly bounded over $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. • Convexity: For each $s \in \{1, 2, 3\}$, there exists a constant $c_{H, s}$ such that $\mathbb{E}\left[ (d^{2} \psi)_{ijt} \right] \geq c_{H, s} > 0$ a.\,s.\ uniformly over $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. Furthermore, there exists a constant $c_{W, s} > 0$ such that for all $s \in \{1, 2, 3\}$, \begin{equation*} \underset{\{v \in \mathbb{R}^{K} \colon \lVert v \rVert_{2} = 1\}}{\min} \; \frac{1}{n_{s}} \sum_{(i, j) \in \mathcal{D}_{s}} \sum_{t = 1}^{T} \mathbb{E}\left[ (d^{2} \psi)_{ijt} \big\{(\ddot{X} v)_{ijt} \big\}^{2} \right] \geq c_{W, s} \quad \text{a.\,s.} \end{equation*} \end{enumerate}

remark[Assumption (ref)] We comment on each part in turn, noting differences from Assumption 4.1 of fw2016 where relevant. \begin{enumerate}[i)] • The asymptotics condition specifies a joint limit in which all panel dimensions diverge at proportional rates. This extends the approach of fw2016 to the multi-dimensional panel setting. For bipartite panels ($s = 1$), two proportionality conditions are required: $N_{1} / T \rightarrow \tau_{1}$ and $N_{2} / T \rightarrow \tau_{2}$. For network panels ($s \in \{2, 3\}$), both cross-sectional dimensions index the same set of $N_{1}$ agents, so the single condition $N_{1} / T \rightarrow \tau_{1}$ suffices. • The sampling condition restricts dependence in two ways: it imposes independence across dyads conditional on $\Phi$, and $\alpha$-mixing within each dyad with coefficients decaying at a polynomial rate. We use $\alpha$-mixing because it is the weakest among the standard mixing conditions in econometrics and is preserved under measurable transformations. The $\alpha$-mixing condition could be replaced by any other form of weak dependence that ensures the implied regularity conditions underlying our results remain satisfied. Compared with fw2016, who impose independence across individuals conditional on $\Phi$, we impose independence across dyadic pairs $(i, j)$. This formulation accommodates both bipartite and network panel data. In the directed network case ($s = 2$), the framework can be extended to allow for reciprocity, that is, pairwise dependence between the dyads $(i, j)$ and $(j, i)$. • The model condition imposes a conditional mean restriction. Conditioning on $\mathcal{X}_{ij}^{t}$ allows the regressors to be predetermined, accommodating dynamic models. Beyond the first moment, the conditional distribution of $y_{ijt}$ is left unrestricted. The normalization $V V^{\prime} \phi_{s}^{0} = \mathbf{0}_{L}$ is required for identification, since constants can be shifted between the components of $\phi$ without affecting $\pi_{ijt}$. If one additionally specifies a fully parametric model, for example, \begin{equation*} y_{ijt} \mid \mathcal{X}_{ij}^{t} \sim f_{Y}\big(y_{ijt} \mid \pi_{ijt}(\beta_{s}^{0}, \mu_{ijt}(\phi_{s}^{0})) \big) \end{equation*} for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$, with $f_{Y}(\cdot)$ a known conditional density or probability mass function, then one can also exploit Bartlett identities (see b1953). These identities can simplify bias and variance expressions and improve the finite-sample performance of debiased estimators, as noted by f2009. We return to this point in Remark (ref) ii). In contrast, fw2016 assume a fully parametric model for their main results (see Assumption 4.1 (iii)) and present results for conditional mean models only in Remark 3. We therefore adopt the more general formulation as our baseline. • The smoothness and moment condition supports the higher-order asymptotic expansion underlying our results. Our expansion is of higher order than that of fw2016, so we require $\psi_{ijt}(\pi_{ijt}(\beta, \mu))$ to be five times continuously differentiable, rather than four. The dominating function $\Psi(z_{ijt})$ uniformly bounds both the regressors and all derivatives of the criterion function up to fifth order. Requiring $\sup_{ijt} \mathbb{E}\left[ \Psi(z_{ijt})^{20 + \nu} \right]$ to be a.\,s.\ uniformly bounded ensures these envelopes have sufficiently many finite moments: $20 + \nu$ here, compared with $8 + \nu$ in fw2016, reflecting our higher-order expansion. • The convexity condition has two components. The first, $\mathbb{E}\left[ (d^{2} \psi)_{ijt} \right] \geq c_{H, s} > 0$, requires the expected second derivative of the criterion function to be uniformly bounded away from zero. This ensures local strict convexity in $\pi_{ijt}$, which is needed for identification and consistency. The second component is a generalized non-collinearity condition: after projecting out the incidental parameters, the within-transformed regressors $\ddot{X}$ must display sufficient variation in all directions $v$. Together, the two components imply strict convexity of the expected objective function over the relevant part of the parameter space, guaranteeing that (ref) has a unique solution. Rather than directly assuming global strict convexity as in fw2016, we impose this more primitive non-collinearity condition. \end{enumerate}

Three-way Fixed Effects Specifications

Before presenting our results for the three-way fixed effects specifications (3.a) and (3.b) in Table (ref), we introduce $\delta_{a}$ to denote a degenerate distribution concentrated at $a$. We also write $\overline{\mathbb{E}}\left[ \cdot \right]$ for the probability limit of the enclosed expression, where the limit is taken as $N_{1}, N_{2}, T \rightarrow \infty$ for $s = 1$ and as $N_{1}, T \rightarrow \infty$ for $s \in \{2, 3\}$. By construction, $\mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$ denotes the weighted residuals of the weighted least-squares problem (ref). As a consequence, for specification (3.b) and all $(i, j) \in \mathcal{D}_{s}$, for example, $\sum_{t = 1}^{T} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{2} \psi)_{ijt} \right] = \mathbf{0}_{K}$ for each $s \in \{1, 2, 3\}$.

\noindentInteracted specification. Theorem (ref) establishes that, for specification (3.b), estimators $\hat{\beta}_{s}$ have degenerate asymptotic distributions for any $s \in \{1, 2, 3\}$.\footnote{More precisely, $\hat{\beta}_{s}$ does not have a nondegenerate asymptotic distribution for any $s \in \{1, 2, 3\}$. We use degenerate asymptotic distributions to avoid double negation.}

theorem[Asymptotic distributions of $\hat{\beta}_{s}$ - Interacted specification] Let Assumptions (ref) hold. Then, for any $s \in \{1, 2, 3\}$, \begin{equation*} \frac{\sqrt{n_{s}}}{\sqrt{T}} \, (\hat{\beta}_{s} - \beta_{s}^{0}) \xrightarrow{d} \delta_{\overline{b}_{s, \infty}} \, , \end{equation*} where \begin{align*} \overline{b}_{1, \infty} \coloneqq& \, - \overline{W}_{1, \infty}^{- 1} \big( \tau_{1}^{\frac{1}{2}} \, \tau_{2}^{- \frac{1}{2}} \, \overline{B}_{1, \alpha, \infty} + \tau_{1}^{- \frac{1}{2}} \, \tau_{2}^{\frac{1}{2}} \, \overline{B}_{1, \gamma, \infty} + \tau_{1}^{\frac{1}{2}} \, \tau_{2}^{\frac{1}{2}} \, \overline{B}_{1, \rho, \infty} \big) \, , \\ \overline{b}_{2, \infty} \coloneqq& \, - \overline{W}_{2, \infty}^{- 1} \big( \overline{B}_{2, \alpha, \infty} + \overline{B}_{2, \gamma, \infty} + \tau_{1} \, \overline{B}_{2, \rho, \infty} \big) \, , \\ \overline{b}_{3, \infty} \coloneqq& \, - \overline{W}_{3, \infty}^{- 1} \big( \overline{B}_{3, \alpha, \infty} + \tau_{1} \, \overline{B}_{3, \rho, \infty} \big) \, , \end{align*} with \begin{align*} & \overline{B}_{1, \alpha, \infty} \coloneqq \overline{\mathbb{E}}\left[ - \frac{1}{N_{1} T} \sum_{i = 1}^{N_{1}} \sum_{t = 1}^{T} \frac{\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ \ddot{x}_{ijt} (d^{2} \psi)_{ijt} (d^{1} \psi)_{ijt} \right]}{\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]} \, + \right. \\ & \quad \left. \frac{1}{2 \, N_{1} T} \sum_{i = 1}^{N_{1}} \sum_{t = 1}^{T} \frac{\big\{\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right]\big\} \sum_{j = 1}^{N_{2}} \mathbb{E}\left[ \big\{(d^{1} \psi)_{ijt}\big\}^{2} \right]}{\big\{\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]\big\}^{2}} \right] \, , \\ & \overline{B}_{1, \gamma, \infty} \coloneqq \overline{\mathbb{E}}\left[ - \frac{1}{N_{2} T} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \frac{\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ \ddot{x}_{ijt} (d^{2} \psi)_{ijt} (d^{1} \psi)_{ijt} \right]}{\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]} \right. \, + \\ & \quad \left. \frac{1}{2 \, N_{2} T} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \frac{\big\{\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right]\big\} \sum_{i = 1}^{N_{1}} \mathbb{E}\left[ \big\{(d^{1} \psi)_{ijt}\big\}^{2} \right]}{\big\{\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]\big\}^{2}} \right] \, , \\ & \overline{B}_{1, \rho, \infty} \coloneqq \overline{\mathbb{E}}\left[ - \frac{1}{N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\sum_{t = 1}^{T} \sum_{t^{\prime} = t}^{T} \mathbb{E}\left[ \ddot{x}_{ijt^{\prime}} (d^{2} \psi)_{ijt^{\prime}} (d^{1} \psi)_{ijt} \right]}{\sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]} \right. \, + \\ & \quad \left. \frac{1}{2 \, N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right]\big\} \Big\{\sum_{t = 1}^{T} \mathbb{E}\left[ \big\{(d^{1} \psi)_{ijt}\big\}^{2} \right]\Big\}}{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]\big\}^{2}} \right] \, , \\ & \overline{W}_{1, \infty} \coloneqq \overline{\mathbb{E}}\left[ \frac{1}{N_{1} N_{2} T} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \, \ddot{x}_{ijt} \, \ddot{x}_{ijt}^{\prime} \right] \right] \, , \\ & \overline{\Sigma}_{1, \infty} \coloneqq \overline{\mathbb{E}}\left[ \frac{1}{N_{1} N_{2} T} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \mathbb{E}\left[ \big\{(d^{1} \psi)_{ijt}\big\}^{2} \ddot{x}_{ijt} \, \ddot{x}_{ijt}^{\prime} \right] \right] \, . \end{align*} The expressions for $s = 2$ and $s = 3$ are stated in Appendix (ref).

remark[Theorem (ref)] To keep the following discussion concise, we focus on bipartite panel data ($s = 1$). \begin{enumerate}[i)] • Each component of the leading incidental parameter bias terms, $\overline{B}_{\alpha, \infty}, \overline{B}_{\gamma, \infty}, \overline{B}_{\rho, \infty}$, originates from estimating the incidental parameters $\alpha$, $\gamma$, and $\rho$, respectively. The corresponding estimators $\hat{\alpha}, \hat{\gamma}$, and $\hat{\rho}$ converge at rates much slower than $\sqrt{N_{1} N_{2} T}$, which causes the inference problem. The biases $\overline{B}_{\alpha, \infty}$ and $\overline{B}_{\gamma, \infty}$ occur even under strict exogeneity. By contrast, $\overline{B}_{\rho, \infty}$ contains additional components that arise when regressors are predetermined, besides those that also arise under strict exogeneity. Given the similarity to Nickell (1981), we refer to this additional bias component as Nickell-type bias. We show \begin{equation*} \sqrt{N_{1} N_{2} T}\big(\hat{\beta} - \beta^{0}\big) \approx \overline{W}_{\infty}^{- 1} \partial_{\beta} \mathcal{L}_{n}\big(\beta^{0}, \phi^{0}\big) - \overline{W}_{\infty}^{- 1} \overline{B}_{\infty} \, , \end{equation*} where \begin{equation*} \overline{B}_{\infty} \coloneqq \overline{B}_{\alpha, \infty} \sqrt{(N_{1} T) / N_{2}}+\overline{B}_{\gamma, \infty} \sqrt{(N_{2} T) / N_{1}} + \overline{B}_{\rho, \infty} \sqrt{(N_{1} N_{2}) / T} \, . \end{equation*} The first term obeys a central limit theorem, i.e., \begin{equation*} \partial_{\beta} \mathcal{L}_{n}\big(\beta^{0}, \phi^{0}\big) \xrightarrow{d} \mathcal{N}\big(0, \overline{\Sigma}_{\infty}\big) \, . \end{equation*} Under Assumption (ref) i), $N_{1} / T \rightarrow \tau_{1} \in (0, \infty)$ and $N_{2} / T \rightarrow \tau_{2} \in (0, \infty)$ as $N_{1}, N_{2}, T \rightarrow \infty$, the second term is of order $\sqrt{N_{1}} + \sqrt{N_{2}} + \sqrt{T}$. Hence, one cannot balance the orders of the bias and variance to obtain a nondegenerate asymptotic distribution. This is a key difference from the existing large-$T$ bias-correction literature, in which uncorrected estimators have nondegenerate asymptotic distributions. Intuitively, the incidental parameters in a three-dimensional panel are estimated from the same number of effective observations as in a standard two-dimensional panel with individual and time fixed effects, even though the overall sample is much larger. Hence, the increased sample size does not improve the convergence rates of the estimators of the incidental parameters. • Using the fully parametric model assumption (e.g., as outlined in Remark (ref)) would allow us to simplify the bias and variance components using Bartlett identities. In particular, by the second Bartlett identity, $\mathbb{E}\left[ \{(d^{1} \psi)_{ijt}\}^{2} \right] = \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$ and $\mathbb{E}\big[\{(d^{1} \psi)_{ijt}\}^{2} \ddot{x}_{ijt} \, \ddot{x}_{ijt}^{\prime}\big] = \mathbb{E}\big[(d_{\mathcal{X}}^{2} \psi)_{ijt} \ddot{x}_{ijt} \, \ddot{x}_{ijt}^{\prime}\big]$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. Hence, $\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ \{(d^{1} \psi)_{ijt}\}^{2} \right] = \sum_{j = 1}^{N_{2}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$, $\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ \{(d^{1} \psi)_{ijt}\}^{2} \right] = \sum_{i = 1}^{N_{1}} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$, and $\sum_{t = 1}^{T} \mathbb{E}\left[ \{(d^{1} \psi)_{ijt}\}^{2} \right] = \sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. In addition, $\overline{\Sigma}_{\infty} = \overline{W}_{\infty}$. \end{enumerate}

remark[(Pseudo-) ML Estimators for Binary Response Models] When $y_{ijt} \in \{0, 1\}$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$, the simplifications from the second Bartlett identity in Remark (ref) ii) follow automatically from the conditional mean assumption. Let $F_{\epsilon}(\cdot)$ be a cumulative distribution function, and define $g_{ijt}(\beta, \mu) \coloneqq F_{\epsilon}(\pi_{ijt}(\beta, \mu))$, $g_{ijt}^{\prime}(\beta, \mu) \coloneqq \partial_{\mu} g_{ijt}(\beta, \mu)$, $h_{ijt}(\beta, \mu) \coloneqq g_{ijt}^{\prime}(\beta, \mu) / (g_{ijt}(\beta, \mu) (1 - g_{ijt}(\beta, \mu)))$, and $h_{ijt}^{\prime}(\beta, \mu) \coloneqq \partial_{\mu} h_{ijt}(\beta, \mu)$. Then, \begin{equation*} \psi_{ijt}(\pi_{ijt}(\beta, \mu_{ijt}(\phi))) = - \big(y_{ijt} \log(g_{ijt}(\beta, \mu_{ijt}(\phi))) + (1 - y_{ijt}) \log(1 - g_{ijt}(\beta, \mu_{ijt}(\phi)))\big) \, , \end{equation*} so that $(d^{1} \psi)_{ijt} = - h_{ijt} (y_{ijt} - g_{ijt})$ and $(d^{2} \psi)_{ijt} = h_{ijt} g_{ijt}^{\prime} - h_{ijt}^{\prime} (y_{ijt} - g_{ijt})$. By the tower property of conditional expectations, the conditional mean assumption $\mathbb{E}\big[y_{ijt} \mid \mathcal{X}_{ij}^{t}\big] = g_{ijt}$, and the fact that $y_{ijt}^{2} = y_{ijt}$, it follows that \begin{align*} \mathbb{E}\left[ \{(d^{1} \psi)_{ijt}\}^{2} \right] =& \, \mathbb{E}\left[ h_{ijt}^{2} (y_{ijt}^{2} - g_{ijt}^{2}) \right] = \mathbb{E}\left[ h_{ijt}^{2} (y_{ijt} - g_{ijt}^{2}) \right] = \mathbb{E}\left[ h_{ijt}^{2} g_{ijt} (1 - g_{ijt}) \right] \\ =& \, \mathbb{E}\left[ h_{ijt} g_{ijt}^{\prime} \right] = \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right] \, . \end{align*}

remark[No predetermined regressors] When none of the regressors is predetermined (i.e., all regressors are strictly exogenous), the conditional mean assumption in Assumption (ref) iii) can be strengthened to \begin{equation*} \mathbb{E}\left[ y_{ijt} \mid \mathcal{X}_{ij}^{T} \right] = g(\pi_{ijt}(\beta_{s}^{0}, \mu_{ijt}(\phi_{s}^{0}))) \end{equation*} for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. Accordingly, the bias term $\overline{B}_{1, \rho, \infty}$ simplifies to \begin{align*} & \overline{B}_{1, \rho, \infty} = \overline{\mathbb{E}}\left[ - \frac{1}{N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\sum_{t = 1}^{T} \mathbb{E}\left[ \ddot{x}_{ijt} (d^{2} \psi)_{ijt} (d^{1} \psi)_{ijt} \right]}{\sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]} \right. \, + \\ & \quad \left. \frac{1}{2 \, N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right]\big\} \Big\{\sum_{t = 1}^{T} \mathbb{E}\left[ \big\{(d^{1} \psi)_{ijt}\big\}^{2} \right]\Big\}}{\big\{\sum_{t = 1}^{T} \mathbb{E}\left[ (d_{\mathcal{X}}^{2} \psi)_{ijt} \right]\big\}^{2}} \right] \, . \end{align*} Strict exogeneity is a natural assumption when the time dimension is replaced by a third cross-sectional dimension, as in a tripartite data structure.

remark[OLS and (Pseudo-)Poisson ML Estimators] For OLS and (pseudo-)Poisson maximum likelihood (PPML) estimators, the bias expressions in Theorem (ref) simplify. \begin{enumerate}[i)] • Several terms vanish because $\sum_{j = 1}^{N_{2}} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right] = \mathbf{0}_{K}$, $\sum_{i = 1}^{N_{1}} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right] = \mathbf{0}_{K}$, and $\sum_{t = 1}^{T} \mathbb{E}\left[ \ddot{x}_{ijt} (d_{\mathcal{X}}^{3} \psi)_{ijt} \right] = \mathbf{0}_{K}$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. For OLS, this follows because $(d^{2} \psi)_{ijt} = 1$ and hence $(d^{3} \psi)_{ijt} = 0$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. For PPML, it follows from the weighted least-squares problem (ref) together with $(d_{\mathcal{X}}^{3} \psi)_{ijt} = (d_{\mathcal{X}}^{2} \psi)_{ijt}$ for all $i, j, t, \lvert \mathcal{D}_{s} \rvert, T$. These results also hold for $s = 2$ and $s = 3$. • If none of the regressors is predetermined (see Remark (ref)), we obtain correctly centered, nondegenerate asymptotic distributions. The first terms of each of the three bias expressions, $\overline{B}_{1, \alpha, \infty}$, $\overline{B}_{1, \gamma, \infty}$, and $\overline{B}_{1, \rho, \infty}$, vanish by the law of iterated expectations, since $\mathbb{E}\big[\ddot{x}_{ijt} (d_{\mathcal{X}}^{2} \psi)_{ijt} \mid \mathcal{X}_{ij}^{T}\big] = \ddot{x}_{ijt} (d_{\mathcal{X}}^{2} \psi)_{ijt}$ and $\mathbb{E}\big[(d^{1} \psi)_{ijt} \mid \mathcal{X}_{ij}^{T}\big] = 0$. These results also hold for $s = 2$ and $s = 3$. Hence, for OLS and PPML estimators under strict exogeneity and any $s \in \{1, 2, 3\}$, \begin{equation*} \sqrt{n_{s}} (\hat{\beta}_{s} - \beta_{s}^{0}) \xrightarrow{d} \operatorname{\mathcal{N}}\big(0, \overline{W}_{s, \infty}^{- 1} \overline{\Sigma}_{s, \infty} \overline{W}_{s, \infty}^{- 1}\big) \, . \end{equation*} The expressions for $s = 2$ and $s = 3$ are stated in Appendix (ref). This also verifies Remark 2 in wz2021. \end{enumerate}

Theorem (ref) shows that uncorrected estimators under the interacted specification have degenerate asymptotic distributions under Assumption (ref), except for the special cases in Remark (ref). Hence, any inference based on these uncorrected estimators is generally invalid. To address this, we propose debiased estimators that enable asymptotically valid hypothesis tests and confidence intervals.

For each $s \in \{1, 2, 3\}$, we adopt an analytical bias correction: the uncorrected estimator $\hat{\beta}_{s}$ is adjusted by an estimate of the bias terms. This estimate is constructed using sample analogues of the expressions in Theorem (ref), with the true parameter values replaced by the corresponding fixed effects estimates. In particular, for $r \in \{1, 2\}$, let $(\widehat{d^{r} \psi})_{ijt}$ denote the sample analogue of $(d^{r} \psi)_{ijt}$. For $r \in \{2, 3\}$, let $(\widehat{d_{\mathcal{X}}^{r} \psi})_{ijt}$ denote the sample analogue of $(d_{\mathcal{X}}^{r} \psi)_{ijt}$. Moreover, for each $k \in \{1, \ldots, K\}$, let

equation[equation omitted — 371 chars of source]

be the sample analogue of (ref). The resulting $n_{s} \times K$ matrix of fitted values is denoted by $\widehat{\mathfrak{X}}$, with its $ijt$-th row defined as $\hat{\mathfrak{x}}_{ijt} \coloneqq (\mu_{ijt}(\hat{\xi}_{1}), \ldots, \mu_{ijt}(\hat{\xi}_{K}))$. We set $\hat{\ddot{X}} \coloneqq X - \widehat{\mathfrak{X}}$, so that $\hat{\ddot{x}}_{ijt} = x_{ijt} - \hat{\mathfrak{x}}_{ijt}$. We denote the debiased estimator by $\tilde{\beta}_{s}$.

Theorem (ref) establishes that the debiased estimators have nondegenerate asymptotic distributions and are asymptotically unbiased. It also establishes consistency of the variance estimators. To estimate the spectral expectations in $\overline{B}_{s, \rho, \infty}$ (from Theorem (ref)), we adapt the truncated spectral density estimator of \textcites{hk2007}{hk2011} and include the finite-sample adjustment of fw2016. The corresponding bandwidth parameter is denoted by $h$.

theorem[Asymptotic distributions of $\tilde{\beta}_{s}$ - Interacted specification] Let Assumptions (ref) hold. Then, \begin{equation*} \widehat{\Sigma}_{s} \xrightarrow{p} \overline{\Sigma}_{s, \infty} \quad and \quad \widehat{W}_{s}^{- 1} \xrightarrow{p} \overline{W}_{s, \infty}^{- 1} \, . \end{equation*} If, in addition, $h / T^{1 / 10} \rightarrow \tau_{h}$ with $0 < \tau_{h} < \infty$ as $h \rightarrow \infty$, then, for any $s \in \{1, 2, 3\}$, \begin{equation*} \sqrt{n_{s}} \, (\tilde{\beta}_{s} - \beta_{s}^{0}) \xrightarrow{d} \operatorname{\mathcal{N}}\big(0, \overline{W}_{s, \infty}^{- 1} \overline{\Sigma}_{s, \infty} \overline{W}_{s, \infty}^{- 1}\big) \, , \end{equation*} where \begin{align*} \tilde{\beta}_{1} =& \, \hat{\beta}_{1} + \widehat{W}_{1}^{- 1} \big( N_{2}^{- 1} \, \widehat{B}_{1, \alpha} + N_{1}^{- 1} \, \widehat{B}_{1, \gamma} + T^{- 1} \, \widehat{B}_{1, \rho} \big) \, , \\ \tilde{\beta}_{2} =& \, \hat{\beta}_{2} + \widehat{W}_{2}^{- 1} \big( N_{1}^{- 1} \, \widehat{B}_{2, \alpha} + N_{1}^{- 1} \, \widehat{B}_{2, \gamma} + T^{- 1} \, \widehat{B}_{2, \rho} \big) \, , \\ \tilde{\beta}_{3} =& \, \hat{\beta}_{3} + \widehat{W}_{3}^{- 1} \big( N_{1}^{- 1} \, \widehat{B}_{3, \alpha} + T^{- 1} \, \widehat{B}_{3, \rho} \big) \, , \end{align*} with \begin{align*} & \widehat{B}_{1, \alpha} \coloneqq - \frac{1}{N_{1} T} \sum_{i = 1}^{N_{1}} \sum_{t = 1}^{T} \frac{\sum_{j = 1}^{N_{2}} \hat{\ddot{x}}_{ijt} (\widehat{d^{2} \psi})_{ijt} (\widehat{d^{1} \psi})_{ijt}}{\sum_{j = 1}^{N_{2}} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}} \, + \\ & \quad \frac{1}{2 \, N_{1} T} \sum_{i = 1}^{N_{1}} \sum_{t = 1}^{T} \frac{\big\{\sum_{j = 1}^{N_{2}} \hat{\ddot{x}}_{ijt} (\widehat{d_{\mathcal{X}}^{3} \psi})_{ijt} \big\} \sum_{j = 1}^{N_{2}} \big\{(\widehat{d^{1} \psi})_{ijt}\big\}^{2}}{\big\{\sum_{j = 1}^{N_{2}} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}\big\}^{2}} \, , \\ & \widehat{B}_{1, \gamma} \coloneqq - \frac{1}{N_{2} T} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \frac{\sum_{i = 1}^{N_{1}} \hat{\ddot{x}}_{ijt} (\widehat{d^{2} \psi})_{ijt} (\widehat{d^{1} \psi})_{ijt}}{\sum_{i = 1}^{N_{1}} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}} \, + \\ & \frac{1}{2 \, N_{2} T} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \frac{\big\{\sum_{i = 1}^{N_{1}} \hat{\ddot{x}}_{ijt} (\widehat{d_{\mathcal{X}}^{3} \psi})_{ijt} \big\} \sum_{i = 1}^{N_{1}} \big\{(\widehat{d^{1} \psi})_{ijt}\big\}^{2}}{\big\{\sum_{i = 1}^{N_{1}} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}\big\}^{2}} \, , \\ & \widehat{B}_{1, \rho} \coloneqq - \frac{1}{N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\sum_{q = 0}^{h} T / (T - q) \sum_{t = q + 1}^{T} \hat{\ddot{x}}_{ijt} (\widehat{d^{2} \psi})_{ijt} (\widehat{d^{1} \psi})_{ij(t - q)}}{\sum_{t = 1}^{T} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}} \, + \\ & \quad \frac{1}{2 \, N_{1} N_{2}} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \frac{\big\{\sum_{t = 1}^{T} \hat{\ddot{x}}_{ijt} (\widehat{d_{\mathcal{X}}^{3} \psi})_{ijt} \big\} \Big\{\sum_{t = 1}^{T} \big\{(\widehat{d^{1} \psi})_{ijt}\big\}^{2}\Big\}}{\big\{\sum_{t = 1}^{T} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt}\big\}^{2}} \, , \\ & \widehat{W}_{1} \coloneqq \frac{1}{N_{1} N_{2} T} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} (\widehat{d_{\mathcal{X}}^{2} \psi})_{ijt} \, \hat{\ddot{x}}_{ijt} \, \hat{\ddot{x}}_{ijt}^{\prime} \, , \\ & \widehat{\Sigma}_{1} \coloneqq \frac{1}{N_{1} N_{2} T} \sum_{i = 1}^{N_{1}} \sum_{j = 1}^{N_{2}} \sum_{t = 1}^{T} \big\{(\widehat{d^{1} \psi})_{ijt}\big\}^{2} \, \hat{\ddot{x}}_{ijt} \, \hat{\ddot{x}}_{ijt}^{\prime} \, . \end{align*} The expressions for $s = 2$ and $s = 3$ are stated in Appendix (ref).

remark[Theorem (ref)] The condition $h / T^{1 / 10} \rightarrow \tau_{h}$ can be relaxed at the cost of stronger moment restrictions. Following fw2016, we recommend reporting results for several bandwidths. Importantly, constructing the debiased estimator does not require knowing whether regressors are predetermined, i.e., the proposed estimators for the bias and variance components are consistent in either case. Finally, the bias formulas in Theorem (ref) formally justify the conjecture of hsw2020, who proposed a bias correction for three-way fixed effects binary choice models based on the heuristic of fw2018 without formal derivation.

remark[Jackknife Inference] As an alternative to the analytical bias correction, one may conduct inference using the jackknife $t$-statistic of h2026_jackknife, at the cost of an additional unconditional homogeneity assumption. Under this assumption, the asymptotic expansion underlying Theorems (ref) and (ref) satisfies both Assumption AD$^{\dagger}$ and Assumption JK$^{\dagger}$ of h2026_jackknife, and the result $J_{q} \xrightarrow{d} t_q$ for $q \in \{1, 2, 3\}$ follows from his Example 3, where $J_{q}$ denotes the jackknife $t_{q}$-statistic. We refer to Appendix (ref) for further details. The jackknife approach is tuning-parameter-free, i.e., no bandwidth selection is necessary, and does not require the explicit bias formulas derived in this paper. However, it still requires knowledge of the asymptotic bias structure to form suitable subsamples.

\noindentNon-interacted specification. Theorem (ref) establishes that, for specification (3.a), estimators $\hat{\beta}_{s}$ have nondegenerate asymptotic distributions and are asymptotically unbiased for any $s \in \{1, 2, 3\}$.

theorem[Asymptotic distributions of $\hat{\beta}_{s}$ - Non-interacted specification] Let Assumptions (ref) hold. Then, for any $s \in \{1, 2, 3\}$, \begin{equation*} \sqrt{n_{s}} \, (\hat{\beta}_{s} - \beta_{s}^{0}) \xrightarrow{d} \operatorname{\mathcal{N}}\big(0, \overline{W}_{s, \infty}^{- 1} \overline{\Sigma}_{s, \infty} \overline{W}_{s, \infty}^{- 1}\big) \end{equation*} and \begin{equation*} \widehat{\Sigma}_{s} \xrightarrow{p} \overline{\Sigma}_{s, \infty} \quad and \quad \widehat{W}_{s}^{- 1} \xrightarrow{p} \overline{W}_{s, \infty}^{- 1} \, . \end{equation*}

Theorem (ref) shows that debiasing is unnecessary for non-interacted specifications. Notably, this result holds even for nonlinear models such as binary choice models, where incidental parameter bias is typically an issue in two-dimensional panels, and for models with predetermined regressors that would otherwise require a Nickell-type bias correction. To explain this finding heuristically, we consider the case of a bipartite panel ($s = 1$). Although $\hat{\alpha}^{\star}$, $\hat{\gamma}^{\star}$, and $\hat{\rho}^{\star}$ converge at rates slower than $\sqrt{N_{1} N_{2} T}$, the bias of $\hat{\beta}$ is of order $1 / (N_{2} T) + 1 / (N_{1} T) + 1 / (N_{1} N_{2})$, which shrinks faster than $1 / \sqrt{N_{1} N_{2} T}$. The bias is therefore asymptotically negligible. This is consistent with the conjecture of fw2018 (Example 11).

Other Fixed Effects Specifications

The three-way fixed effects specifications are the most general cases; the asymptotic results for the remaining 15 one- and two-way structures in Table (ref) follow as special cases. We discuss each in turn.

\noindentInteracted specifications (one- and two-way). For specifications with only interacted unobserved heterogeneity (Table (ref) (1.1.b)--(1.3.b), (2.1.b)--(2.3.b)), the asymptotic distribution remains degenerate, and the bias components from Theorem (ref) induced by the respective incidental parameters carry over. For illustration, consider specification (2.1.b), $\alpha_{it} + \gamma_{jt}$. Let Assumptions (ref) hold. For any $s \in \{1, 2, 3\}$,

equation*[equation* omitted — 142 chars of source]

where

align*[align* omitted — 577 chars of source]

These expressions are the same as in Theorem (ref), except that the $\rho$-component is absent and quantities such as $\xi_{k}^{0}$ are adjusted for the new specification. Bias correction and inference proceed as in Theorem (ref), with all components adjusted accordingly.

\noindentNon-interacted specifications (one- and two-way). For specifications with only non-interacted unobserved heterogeneity (Table (ref) (1.1.a)--(1.3.a), (2.1.a)--(2.3.a)), the asymptotic distributions are nondegenerate and centered at zero, as in Theorem (ref). For any $s \in \{1, 2, 3\}$,

equation*[equation* omitted — 222 chars of source]

\noindentMixed specifications. For specifications with both non-interacted and interacted unobserved heterogeneity (Table (ref) (2.1.c)--(2.3.c)), the slowly decaying bias from the interacted component dominates the faster-decaying bias from the non-interacted component. Hence, the estimators of $\beta_{s}^{0}$ have a degenerate asymptotic distribution, as in Theorem (ref). As an example, consider specification (2.3.c), $\rho_{t}^{\star} + \rho_{ij}$, for which the bias component for $\rho$ from Theorem (ref) remains. Let Assumptions (ref) hold. For any $s \in \{1, 2, 3\}$,

equation*[equation* omitted — 142 chars of source]

where

equation*[equation* omitted — 466 chars of source]

We again follow Theorem (ref) to form the debiased estimator, with bias and variance components adjusted for the respective specification.

Simulation Experiments

We conduct simulation experiments to study the finite-sample behavior of the uncorrected and debiased three-way fixed effects maximum likelihood estimators. We examine relative bias (in percent, relative to the truth), bias relative to standard deviation, and coverage rates of 95% nominal confidence intervals.\footnote{We calculate the bias as $R^{-1} \sum_{r=1}^{R} (\hat{\beta}_{r} - \beta^{0})$, where $R$ is the number of Monte Carlo replications and $\hat{\beta}_{r}$ is the coefficient in the $r$-th replication.} We adapt the dynamic data generating process of fw2016 to a probit model for bipartite panels with three sets of unobserved effects. For $i \in \{1, \ldots, N_{1}\}$, $j \in \{1, \ldots, N_{2}\}$, and $t \in \{0, \ldots, T\}$,

equation*[equation* omitted — 319 chars of source]

where $u_{ijt}, v_{ijt} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1)$ and $x_{ijt} = \mu_{ijt}(\phi) + v_{ijt}$. In DGP I, we consider interacted unobserved effects, $\mu_{ijt}(\phi) = \alpha_{it} + \gamma_{jt} + \rho_{ij}$, with $\alpha_{it}, \gamma_{jt}, \rho_{ij} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 24)$. In DGP II, we consider non-interacted fixed effects, $\mu_{ijt}(\phi) = \alpha_{i}^{\star} + \gamma_{j}^{\star} + \rho_{t}^{\star}$, with $\alpha_{i}^{\star}, \gamma_{j}^{\star}, \rho_{t}^{\star} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 24)$. We set $\beta_{y} = 0.5$ and $\beta_{x} = 1$, and generate samples with $N_{1} \in \{60, 120, 240\}$, $N_{2} = N_{1} / 2$, and $T(N_{1}) = N_{1} / 5$. This design ensures that $N_{1}$, $N_{2}$, and $T$ grow at constant relative rates, consistent with our asymptotic analysis. All results are based on $10{,}000$ simulated samples for each $N_{1}$.

The results for DGP I are presented in Table (ref). For the debiased estimator, we report results for bandwidth values $h \in \{0, 1, 2, 3, 4\}$. A bandwidth of $h = 0$ corrects only for the incidental parameter bias; $h > 0$ additionally corrects for the Nickell-type bias.

table[table omitted — 2,766 chars of source]

We first examine the results for the strictly exogenous regressor $\beta_x$, shown in the left panel of Table (ref). For the smallest sample size ($N_1 = 60, N_2 = 30, T = 12$), the uncorrected estimator has a relative bias of $21.868\%$. This bias decreases with sample size, as predicted by theory: the bias is of order $1 / N_{1} + 1 / N_{2} + 1 / T$ and falls to $4.229\%$ at the largest sample size ($N_1 = 240, N_2 = 120, T = 48$). Despite this decrease, the bias remains large relative to the estimator's dispersion. The bias-to-standard-deviation ratio in fact increases with sample size, so the problem worsens in relative terms. Accordingly, coverage rates are zero throughout. As predicted by our asymptotic theory, the uncorrected estimator therefore exhibits a more severe bias problem than those studied in, for example, fw2016 and wz2021.

The left panel also shows results for the debiased estimator. Since $x_{ijt}$ is strictly exogenous, correcting only for the incidental parameter bias ($h = 0$) suffices in principle. Nevertheless, we also report results for $h > 0$ to reflect the more realistic scenario in which the researcher does not know whether regressors are strictly exogenous. The debiased estimator with $h = 0$ outperforms the uncorrected estimator on every metric in every sample. At the smallest sample size, the bias falls from $21.868\%$ to $-1.584\%$ and coverage improves from zero to $82.3\%$. The debiased estimator with $h > 0$ performs reliably for strictly exogenous regressors as well. Notably, even for the strictly exogenous regressor, $h = 1$ outperforms $h = 0$ at every sample size. For example, at ($N_{1} = 60, N_{2} = 30, T = 12$), coverage improves from $82.3\%$ to $92.8\%$. This might be because, in finite samples, the incidental parameter bias and the Nickell-type bias are not as cleanly separated as asymptotic theory suggests, due to correlation between the strictly and predetermined regressor. This provides further motivation for the recommendation to use $h > 0$ in practice.

The right panel of Table (ref) presents results for the predetermined regressor $y_{ij(t-1)}$. The qualitative pattern is the same as for $\beta_x$, but the biases of the uncorrected estimator are considerably larger due to the additional Nickell-type bias. Debiasing with $h > 0$ substantially reduces these biases and improves coverage. Debiasing with $h = 0$, which corrects for the incidental parameter bias but not the Nickell-type bias, yields unreliable inference, similar to no correction at all.

We turn next to DGP II, with results in Table (ref). As established by our theory, the uncorrected estimator with non-interacted fixed effects is consistent and has a correctly centered asymptotic distribution, even for predetermined regressors, so no bias correction is needed. The left and right panels confirm this: in both cases, biases relative to standard deviation are negligible and coverage rates are close to the nominal $95\%$ level.

table[table omitted — 1,005 chars of source]

Overall, the simulation experiments support the validity of our asymptotic results in samples of sufficient size. Based on the DGP I results, we recommend that practitioners use a debiasing method with $h > 0$ whenever it is unclear whether regressors are strictly exogenous.

Empirical Application

We illustrate our inferential procedure using bilateral trade data on plastic articles (HS6 code 392690), a broad category of standardized plastic components. We adapt the dynamic network formation model of g2016 to a directed network setting. This product category is well-suited for two reasons. First, it has the highest country participation in the data, as nearly every country in the world produces or consumes products in this category. Second, and more importantly, plastic articles of this type can in principle be produced by any country, implying that the probability of a trade link forming is bounded away from zero for all country pairs. This satisfies Assumption (ref) v), which requires $0 < \mathbb{P}\big(y_{ijt} \mid \mathcal{X}_{ij}^{t}\big) < 1$ for all $i, j, t, \lvert \mathcal{D}_{2} \rvert, T$, i.e., the network is dense.

\noindentDynamic model of network formation. g2016 proposes a dynamic model of network formation for undirected network panel data in which the probability of a pair $(i,j)$ forming a link at time $t$ increases if (i) the pair was already directly connected in the previous period (state dependence), (ii) the pair shared many links in common in the previous period (a taste for transitivity), and (iii) the pair shares unobserved time-invariant characteristics (homophily). Channels (i) and (ii) are observed and enter the utility function directly through the lagged link indicator $y_{ij(t - 1)}$ and the lagged count of two-paths $r_{ij(t - 1)}$, respectively.

Trade networks are inherently directed: exporting to a country is structurally distinct from importing from it, so the undirected framework of g2016 is not directly applicable. Directed network data has two consequences for the model. First, the single lagged two-path count $r_{ij(t - 1)}$ of g2016 decomposes into four structurally distinct triadic variables:

align*[align* omitted — 332 chars of source]

In trade networks, these triadic variables can be interpreted as supply chain transitivity, cyclic closure, shared upstream suppliers, and shared export markets, respectively. Second, the directed setting introduces a natural analog of undirected link persistence: the lagged reverse-direction link $y_{ji(t-1)}$, which we refer to as dynamic reciprocity. If country $j$ exported to country $i$ in the previous period, existing bilateral infrastructure, such as shipping routes, lowers the cost of establishing a link in the reverse direction.

Based on these considerations, we model whether exporter $i$ and importer $j$ have a trading relationship at time $t$ as:

equation[equation omitted — 275 chars of source]

where $\alpha_{it}$ and $\gamma_{jt}$ capture time-varying out-degree and in-degree unobserved heterogeneity, $\rho_{ij}$ is pair-specific unobserved heterogeneity (homophily), and $u_{ijt} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0,1)$ is a pair- and time-specific shock. We include $r_{1, ij(t-1)}$ and $r_{3,ij(t-1)}$ as the two triadic variables with the clearest economic motivation in a trade setting: supply chain transitivity, where an indirect export chain $i \to k \to j$ raises the probability of a direct link $i \to j$, and shared upstream suppliers, where two countries sourcing from a common supplier ($k \to i$ and $k \to j$ ) are more likely to establish direct bilateral trade. We exclude the cyclic closure term $r_{2, ij(t-1)}$, which has no compelling trade interpretation for a homogeneous good, and the shared-export-market term $r_{4, ij(t-1)}$, which is theoretically ambiguous due to competing effects.\footnote{When we include $r_{2,ij(t-1)}$ and $r_{4,ij(t-1)}$ in our empirical specification, we find insignificant coefficients close to zero.}

\noindentData. We use the BACI international trade database, which provides harmonized bilateral trade flows at the HS6 level.\footnote{BACI reconciles import and export declarations reported to UN Comtrade. The reconciliation procedure effectively discards flows that are very small or inconsistently reported; see baci_database for details.} Our sample covers the period 1995--2019 ($T = 25$), includes $N_{1} = 233$ countries, and defines $y_{ijt} = 1$ if country $i$ reports positive exports to country $j$ in year $t$, and zero otherwise. The panel is unbalanced.

The network becomes markedly denser over the sample period: in 1995, 5,796 of 45,156 directed country pairs traded (12.8%), rising to 11,461 of 49,506 pairs (23.2%) in 2006 and 14,823 of 50,850 pairs (29.2%) in 2019. This near-doubling of network density over 25 years reflects the broad globalization of trade in this product category. It also motivates the dynamic specification in (ref): the gradual densification suggests that new trade relationships do not form in isolation, but build on existing ones through state dependence and triadic closure, as new exporters enter markets where indirect trade chains or shared suppliers are already in place.

\noindentResults. We estimate (ref) by probit maximum likelihood. Table (ref) reports the estimates, comparing the uncorrected estimator against the debiased estimator for bandwidths $h \in \{1, \dots , 4\}$.

table[table omitted — 1,533 chars of source]

The debiased estimates are stable across all values of $h$, which serves as a sensitivity check on the bandwidth choice. All estimates are positive and statistically significant under both estimators. The estimate of $\beta_{1}$ indicates strong persistence in bilateral trade relationships, consistent with sunk costs at the firm level generating state dependence hsw2020. The uncorrected estimate of $\beta_{1}$ is substantially smaller than the debiased estimates across all bandwidths, consistent with a downward Nickell-type bias. The estimate of $\beta_{2}$ confirms that an existing import relationship from $j$ to $i$ raises the propensity of a subsequent export link from $i$ to $j$, consistent with bilateral trade infrastructure lowering costs in both directions. The estimates of $\beta_{3}$ and $\beta_{4}$ provide evidence for both triadic channels: supply chain transitivity and shared upstream suppliers each raise the probability of a direct trade link forming, in line with global value chain mechanisms. The comparable magnitudes of the estimates for $\beta_{3}$ and $\beta_{4}$ suggest that both channels contribute similarly to network densification. For $\beta_{2}$, $\beta_{3}$, and $\beta_{4}$, the difference between the uncorrected and debiased estimates is less pronounced but still relatively large compared to the very small standard errors.

Concluding Remarks

We provide the first extensive theoretical analysis of fixed effects M-estimators for three-dimensional panel data across different data structures and specifications of multi-way unobserved heterogeneity. Our analysis reveals a sharp dichotomy: non-interacted specifications are free of asymptotic bias, whereas interacted specifications suffer from a severe inference problem manifested in degenerate asymptotic distributions. The choice of heterogeneity specification therefore has first-order consequences for inference. We resolve the inference problem with analytical debiased estimators that have nondegenerate, correctly centered asymptotic distributions. Simulation evidence confirms good finite-sample performance, and our empirical application shows that the correction is economically meaningful.

In ongoing work, we explore the following extensions. For some models, in particular binary choice models, we extend our analysis to average partial effects. In addition, we conjecture that for some specifications, it is possible to relax the asymptotics, e.g., to derive results under fixed $T$ asymptotics. Moreover, drawing on the connection between fw2016 and wz2021, we conjecture that our results will be helpful in studying the properties of (pseudo-)poisson maximum likelihood estimators for four-dimensional panel data models with triple interacted unobserved heterogeneity under a strict exogeneity assumption.

\printbibliography