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.
66,102 characters · 15 sections · 106 citation commands
Debiased Fixed Effects Estimation of Binary Logit Models with Three-Dimensional Panel Data
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1}
Even after more than 75 years since its discovery by ns1948, the incidental parameter problem, originally a specific inconsistency problem that occurs in many fixed effects estimators for nonlinear models, remains a highly studied topic in panel data econometrics. Most work on the incidental parameter problem focuses on classical panel data sets, i.e. two-dimensional panels, where cross-sectional units $N$ are observed over several time periods $T$. Here, the incidental parameter problem arises under asymptotics where $N$ tends to infinity, while $T$ is fixed. One approach to tackle the incidental parameter problem is to rely on a different asymptotic framework, usually $N, T \rightarrow \infty$. This fixes the inconsistency problem but introduces an asymptotic bias problem. If this asymptotic bias problem is not properly addressed, e.g.\ by applying suitable bias corrections, the fixed effects estimator becomes unreliable for drawing inferences. In this paper, we focus on bias correction approaches developed under large-$T$ asymptotics. However, it is worth mentioning another important strand of literature that focuses on developing fixed-$T$ consistent estimators. These estimators typically rely on eliminating the unobserved effects from the model by differencing or conditioning on sufficient statistics.\footnote{Examples for binary logit models are, among others, r1960, a1970, c1980, hk2000, or hw2022.}
Multi-dimensional panel data is becoming increasingly common in empirical research, as data becomes more granular. One example of a three-dimensional panel is data on bilateral network activities observed over time, which is often used in international trade research. More specifically, researchers may study trade flows between $I$ exporting countries and $J$ importing countries over $T$ years. This multi-dimensionality allows to control for richer sources of unobserved heterogeneity, leading to model specifications with multi-way fixed effects. For example, researchers in international trade often control for unobserved heterogeneity at the exporter-time, importer-time, and importer-exporter levels in their empirical analyses.\footnote{Controlling for unobserved heterogeneity at the exporter-time, importer-time, and importer-exporter levels is, among others, recommended in hm2014.} However, the asymptotic properties of the corresponding fixed effects estimators are largely unknown, with a few exceptions.
In this paper, we derive the asymptotic properties of fixed effects estimators for static binary logit models for three-dimensional panels, where the three sets of unobserved effects enter additively into the linear index as $\alpha_{it} + \gamma_{jt} + \rho_{ij}$. Under asymptotics, where all panel dimensions grow large and $I \sim J \sim T$, we use expansions to characterize the leading bias term and suggest an appropriate bias correction. We show that the order of the bias is $1 / I + 1 / J + 1 / T$, confirming a conjecture of fw2018 that has not been proven yet.\footnote{The conjecture is based on a heuristic formula developed by fw2018 as part of a review of recent advances in fixed effects estimation.} Moreover, we confirm the correctness of the expressions conjectured and proposed by hsw2020 for bias-corrected estimators.
Our main finding, which distinguishes our case from all other cases studied in the bias correction literature, is that the inference problem is more severe. In most cases, the order of the bias and standard deviation of fixed effects estimators are the same, allowing authors to derive non-degenerate asymptotic distributions for the uncorrected estimators, which are centered around distorted expected values. Bias corrections then center the asymptotic distributions correctly around zero. Our inference problem is particularly severe because the leading bias of our fixed effects estimator, $1 / I + 1 / J + 1 / T$, is of a higher order than its standard deviation, $1 / \sqrt{IJT}$, which leads to a degenerating asymptotic distribution. Therefore, developing a debiased estimator is particularly important. In simulation experiments we confirm the severity of the inference problem without bias correction, as confidence intervals constructed around the uncorrected estimator almost never cover the true model parameters, even in large samples. Our proposed bias correction is effective in improving the inferential accuracy of the fixed effects estimator. An empirical example from international trade shows that debiased estimates can differ substantially from uncorrected estimates in real-world applications. Thus, our findings have important implications for empirical researchers, highlighting the need for bias correction to obtain reliable inference. We expect our results to generalize to other link functions and dynamic models. To simplify the analysis, we study the asymptotic properties of fixed effects estimators using the binary logit model as an example, following arguments of cfw2020.\footnote{The simplification mainly comes from the fact that high-order derivatives of the log-likelihood function no longer depend on the outcome variable.} We plan to generalize our asymptotic analysis to (dynamic) nonlinear models with concave objective functions in the future.
A large part of the previous large-$T$ literature has focused on providing solutions for classical panel data models with individual effects, see among others, hk2002, l2002, w2002, s2003, hn2004, c2007, ab2009, bh2009, f2009, hk2011, dj2015, ks2016, p2019, sst2021, hj2022, or s2023. These proposed solutions differ in various ways, including the assumptions they make, the methods used to derive them, and the types of corrections they propose. We refer the reader to ah2007 and fw2018 for comprehensive reviews of this strand of literature. fw2016 advance the literature by developing solutions for nonlinear panel data models that can account for both individual and time effects.\footnote{The authors' analysis is not limited to panels with time as the second panel dimension. It can also be applied to other two-dimensional panels, such as panels where the second panel dimension is another cross-section, such as countries or industries. For example, their analysis could be used to study a cross-section of bilateral trade flows between countries or a cross-section of patent citations between industries.} This is a major contribution, as accounting for both types of effects is challenging. We will discuss their contribution in more detail in Section (ref), as it is essential for the derivation of our results. jo2019 also studied nonlinear panel data models with individual and time effects. However, instead of correcting the bias in the asymptotic distribution, as fw2016, they proposed a likelihood correction approach.\footnote{Recently, lms2023 presented a related approach to jo2019 and additionally proved the asymptotic properties of the corrected likelihood and the test statistics of the trinity tests of maximum likelihood estimation (Wald, Lagrange-multiplier, and Likelihood-ratio test).} Contemporaneously with the development of fw2016's bias correction, c2017 extended the conditional logit estimator of r1960 and c1980 to handle two-way fixed effects.\footnote{She applied her estimator to a cross-sectional model of bilateral export probability. Her estimator requires both panel dimensions to grow large, as proven later by j2018, who derived its asymptotic properties. Because her approach is already computationally demanding in bilateral cross-sections, applying it to bilateral panels or extending it to three-way fixed effects may not be feasible. Moreover, it is not possible to generalize her approach to other nonlinear models and weakly exogenous regressors.} There is only little research on fixed effects estimators for multi-dimensional panel models with multiple unobserved effects, although multi-dimensional panels are increasingly common in empirical studies. wz2021 is the only other paper apart from ours that has theoretically analyzed the properties of a fixed effects estimator for a three-dimensional panel model with three additive and overlapping unobserved effects. In particular, they studied the properties of the fixed effects Pseudo-Poisson estimator for the gravity model, the workhorse model of international trade.\footnote{yh2023 theoretically analyzed an estimator for an alternative gravity model with three multiplicative, instead of additive, and overlapping unobserved effects by extending the generalized method of moments (GMM) estimation strategy proposed by j2017.} Contrary to us, wz2021 can exploit a unique property of the Poisson model to eliminate $\rho_{ij}$ from the linear index. This essentially reduces the problem to analyzing a two-way fixed effects model, which allows them to rely on the asymptotic analysis of fw2016. As a consequence, they only require $I$ and $J$ to grow to infinity, while we require all three panel dimensions to grow.
The rest of the paper is organized as follows. Section (ref) introduces the model and the fixed effects estimator. Section (ref) presents the asymptotic theory and Section (ref) discusses the key differences from previous studies. Sections (ref) and (ref) report results of simulation experiments and an empirical example. Section (ref) concludes.
We observe three-dimensional panel data $\{(y_{ijt}, x_{ijt}) \colon i \in \mathcal{I}, \, j \in \mathcal{J}, \, t \in \mathcal{T}\}$, where $y_{ijt}$ is a binary outcome variable, $x_{ijt}$ is a vector of strictly exogenous explanatory variables, $\mathcal{I} = \{1, \ldots, I\}$, $\mathcal{J} = \{1, \ldots, J\}$, and $\mathcal{T} = \{1, \ldots, T\}$. We consider the following semi-parametric binary logit model with additive unobserved effects:
where $\operatorname{\mathbbold{1}}\{\cdot\}$ is an indicator function, $x = (x_{111}, \ldots, x_{IJT})$, $\alpha = (\alpha_{11}, \ldots, \alpha_{IT} )$, $\gamma = (\gamma_{11}, \ldots, \gamma_{JT} )$, $\rho = (\rho_{11}, \ldots, \rho_{IJ})$, $\epsilon_{ijt}$ is an idiosyncratic error term, and $F_{\epsilon}$ is the logistic cumulative distribution function. Further, $\beta$ is a $K$-dimensional vector of model parameters, and $\alpha$, $\gamma$, and $\rho$ are $IT$-, $JT$-, and $IJ$-dimensional vectors of unobserved effects, respectively. We interpret the model as semi-parametric because we do not make assumptions about the relationship between the unobserved effects and the explanatory variables, nor do we make assumptions about the distributions of the unobserved effects.
Below we present two examples of three-dimensional panel data sets in which our model could be applied.
We collect the incidental parameters in the vector $\phi = (\alpha, \gamma, \rho)$ and estimate them along with the model parameters $\beta$ by minimizing the following constrained negative log-likelihood function:
where $0 < c_{1} < \infty$, $\mu_{ijt}(\beta, \phi) = \mu(x_{ijt}^{\prime} \beta + w_{ijt}^{\prime} \phi)$, $\mu(z) = (1 + \exp(- z))^{- 1}$ is the logistic cumulative distribution function, and the matrix $w$ is a collection of $IT + JT + IJ$ indicator variables arising from “dummy encoding” the following interactions of the three panel indices: $i \times t$, $j \times t$, and $i \times j$. The matrix $v$ imposes constraints on the incidental parameters $\phi$ to ensure uniqueness of the solution of the optimization problem and therefore the second term in $L(\beta, \phi)$ acts as “penalty” term. Essentially, the penalty term prevents $w$ from being rank-deficient and therefore plays an important role in ensuring the invertibility of the incidental parameter Hessian,
Finally, note that specific choices for $c_{1}$ and the scaling factor are important for our asymptotic analysis.
To understand the rank deficiency problem problem and the derivation of the constraints, it is instructive to have a closer look at the linear index, $x_{ijt}^{\prime} \beta + \alpha_{it} + \gamma_{jt} + \rho_{ij}$. The incidental parameters enter additively into the linear index which makes the log-likelihood invariant to certain parameter transformations. For example, the linear index is invariant to adding a constant $c_{t}$ to all $\alpha_{it}$ while subtracting it from all $\gamma_{jt}$. Therefore, we introduce $T$ constraints $\sum_{i = 1}^{I} \alpha_{it} = \sum_{j = 1}^{J} \gamma_{jt}$ for $t = \{1, \ldots, T\}$, or in matrix notation $(1_{I} \otimes \operatorname{\mathbb{I}}_{T})^{\prime} \alpha = (1_{J} \otimes \operatorname{\mathbb{I}}_{T})^{\prime} \gamma$. Similarly, 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 unaffected, leading to $I$ constraints $\sum_{t = 1}^{T} \alpha_{it} = \sum_{j = 1}^{J} \rho_{ij}$ for $i = \{1, \ldots, I\}$ and $J$ constraints $\sum_{t = 1}^{T} \gamma_{it} = \sum_{i = 1}^{I} \rho_{ij}$ for $j = \{1, \ldots, J\}$. In matrix notation, these constraints translate to $(\operatorname{\mathbb{I}}_{I} \otimes 1_{T})^{\prime} \alpha = (\operatorname{\mathbb{I}}_{I} \otimes 1_{J})^{\prime} \rho$ and $(\operatorname{\mathbb{I}}_{J} \otimes 1_{T})^{\prime} \gamma = (1_{I} \otimes \operatorname{\mathbb{I}}_{J})^{\prime} \rho$. To impose all constraints simultaneously, we define
such that $v^{\prime} \phi = 0$ characterizes the system of linear equality constraints. Because one of the constraints in $v$ is implied by all other constraints, the rank of $v$ reduces to $T + I + J - 1$. However, for our asymptotic analysis, it is more convenient to work with the $(IT + JT + IJ) \times (T + I + J)$ matrix $v$. In practice there can be several choices for $v$ that work. Perhaps the most familiar way is to set specific incidental parameters to zero, like excluding one time effect in models with individual and time effects without common intercept for classical panels. As in fw2016, for our asymptotic analysis it is however important to choose a specific normalization which is easier to work with.\footnote{More precisely, our normalization ensures that the inverse of the incidental parameter Hessian, defined in (ref), becomes block diagonal which helps us to bound its spectral norm in Lemma (ref).}
Since our primary interest is the estimation of the model parameters $\beta$, i.e.\ we treat the incidental parameters $\phi$ as high-dimensional nuisance parameters, we define the (profile) maximum likelihood estimator as
In this section, we derive the asymptotic properties of the maximum likelihood estimator $\hat{\beta}$, defined in (ref), using an asymptotic framework where all three panel dimensions, $I$, $J$, and $T$, simultaneously grow to infinity. To simplify the notation and make our asymptotic analysis more concise, we follow wz2021 and set $N = I = J$.
We make the following assumptions.
Before presenting the asymptotic distribution, we first need to introduce some additional notation. Let $\mu^{\langle 1 \rangle}(z) = \partial_{z} \mu(z)$, $\mu^{\langle 2 \rangle}(z) = \partial_{z^2} \mu(z)$, and $\mu^{\langle 3 \rangle}(z) = \partial_{z^3} \mu(z)$ denote the first-, second-, and third-order derivatives of the logistic cumulative distribution function $\mu(\cdot)$. Further, we define $\mu_{ijt} = \mu(x_{ijt}^{\prime} \beta^{0} + w_{ijt}^{\prime} \phi^{0})$. The definitions of $\mu_{ijt}^{\langle 1 \rangle}$, $\mu_{ijt}^{\langle 2 \rangle}$, and $\mu_{ijt}^{\langle 3 \rangle}$ follow accordingly. For every regressor $x_{k}$, we define $\tilde{x}_{k} = x_{k} - w^{\prime} \phi_{k}^{\ast}$, where
are the coefficients of a weighted least-squares problem. The residuals $\tilde{x}_{k}$ stem from Legendre transforms that we use to project out the incidental parameters from the asymptotic expansions (details about the transformation are provided in Appendix (ref)). Furthermore, we define the leading asymptotic bias
where
is the normalized profile Hessian, $\overline{W} = \mathbb{E}\left[ W \right]$, and
are bias components that arise from estimating the incidental parameters $\alpha$, $\gamma$, and $\rho$, respectively.
We establish in the following Theorem that $\hat{\beta}$ has a degenerating asymptotic distribution.
We proof Theorem (ref) in Appendix (ref).
The following theorem states that after correcting the leading asymptotic bias, we obtain a correctly centered non-degenerate asymptotic distribution.
The proof of Theorem (ref) is provided in Appendix (ref).
Theorem (ref) shows that the uncorrected estimator has a degenerating asymptotic distribution. Consequently, standard maximum likelihood inference, i.e.\ confidence regions constructed around the uncorrected estimator are in general invalid. However, the inference problem can be resolved, as shown in Theorem (ref), by subtracting the leading asymptotic bias $b$ defined in (ref).
In the following, we use plug-in estimates of $B_{\alpha}$, $B_{\gamma}$, $B_{\rho}$, and $\overline{W}$ to construct a bias-corrected estimator $\tilde{\beta}$. Let $\hat{\mu}_{ijt} = \mu(x_{ijt}^{\prime} \hat{\beta} + w_{ijt}^{\prime} \hat{\phi})$, $\hat{\mu}_{ijt}^{\langle 1 \rangle} = \mu^{\langle 1 \rangle}(x_{ijt}^{\prime} \hat{\beta} + w_{ijt}^{\prime} \hat{\phi})$, and $\hat{\mu}_{ijt}^{\langle 2 \rangle} = \mu^{\langle 2 \rangle}(x_{ijt}^{\prime} \hat{\beta} + w_{ijt}^{\prime} \hat{\phi})$. Further, for every regressor $x_{k}$, we define $\hat{\tilde{x}}_{k} = x_{k} - w \hat{\phi}_{k}^{\ast}$, where
are the coefficients of a weighted least-squares problem. Then, a bias-corrected estimator is constructed as
where
and $\hat{\beta}$ is the uncorrected maximum likelihood estimator defined in (ref).
The following Lemma shows that the estimators for the various bias components and the expected normalized profile Hessian, proposed in (ref), are consistent.
We proof Lemma (ref) in Appendix (ref).
To better align the results presented in the previous section with the results from the previous literature, we compare our results with the ones from the two most related papers, fw2016 and wz2021, and discuss where the differences come from.
Compared to classical panel models with only individual fixed effects, e.g.\ hn2004, f2009, and hk2011, models with additional fixed effects add further complications to the asymptotic analysis. First, it is not possible to express the log-likelihood function as a sum of individual log-likelihood contributions, where each log-likelihood contribution depends only on a fixed-dimensional set of parameters. This strategy was proposed by hn2004 to deal with the infinite-dimensional parameter space and is a common strategy in the panel data literature. Second, the incidental parameter Hessian $\partial_{\phi \phi^{\prime}} L(\beta, \phi)$ is no longer diagonal, which complicates, for example, bounding some quantities in the asymptotic expansion. fw2016 solve both issues for classical panel data models with individual and time effects. For the first issue, they propose a projection method based on Legendre transforms of the log-likelihood function to obtain asymptotic expansions which do not depend on the incidental parameters. For the second issue, they establish an approximation argument for the inverse of the incidental parameter Hessian, in which they show that asymptotically the inverse is a (weakly) diagonally dominant matrix (see Lemma D.1 in fw2016). Thus, it can be uniformly approximated by a diagonal matrix and has off-diagonal elements that are sufficiently small. This approximation argument is particularly important to show that the asymptotic bias can be “decoupled”, i.e.\ expressed as the sum of two bias components, one for each set of fixed effects in the model specification. Finally, the most related paper is wz2021, that analyzes the properties of a fixed effects (Pseudo-)Poisson estimator for three-dimensional panels with the same linear index specification as in this paper under similar assumptions. However, contrary to us, the authors can exploit a unique property of the (Pseudo-)Poisson model that allows them to profile-out $\rho$ from the (pseudo-)log-likelihood function. This essentially turns their three-way model into a stacked two-way model, with only $\alpha$ and $\gamma$ as incidental parameters, and allows them to rely on the results of fw2016 for their asymptotic analysis. As a consequence, wz2021 only need $N = I = J$ to grow to infinity and $T$ can be fixed.
Both, fw2016 and wz2021, derive non-degenerate asymptotic distributions of $\hat{\beta}$,
where $r_{n}$ is the convergence rate of $\hat{\beta}$, $b$ is the constant leading bias, and $V$ is an asymptotic covariance matrix. Here, the order of the normalized asymptotic bias $b / r_{n} = \mathcal{O}_{P}(r_{n}^{- 1})$ and the convergence rate of $\hat{\beta}$ are exactly balanced yielding a non-degenerate but distorted asymptotic distribution.
Figure (ref) illustrates the non-degenerate asymptotic distributions derived in both papers. It shows the empirical densities of the normalized differences of the uncorrected estimators and the true parameter values, $r_{n} (\hat{\beta} - \beta^{0})$, for a logit model with individual and time effects (left panel), as studied in fw2016, and for a Pseudo-Poisson model with exporter-time, importer-time, and exporter-importer effects (right panel), as studied in wz2021. The figure is based on simulated data for different sample sizes. We use the data generating processes from the corresponding papers. The figure shows that as the sample size increases, the distribution of the uncorrected estimator converges to a normal distribution centered around the bias. Both papers propose bias corrections to re-center the asymptotic distribution properly to ensure reliable inference.
In contrast, Theorem 2 reveals that in our case, the order of the normalized asymptotic bias $b / \sqrt{NT} = \mathcal{O}_{P}((NT)^{- 1 / 2})$ and the convergence rate of $\hat{\beta}$, $r_{n} = N\sqrt{T}$, are not balanced. More precisely, the normalized bias shrinks slower than the standard deviation of the uncorrected estimator, resulting in a degenerating asymptotic distribution.
Figure (ref) illustrates the balancing problem using two different normalizing constants for the differences of the uncorrected estimators and the true parameter values: (left panel) $N \sqrt{T} (\hat{\beta} - \beta^{0})$ and (right panel) $\sqrt{NT} (\hat{\beta} - \beta^{0})$. The figure shows empirical densities of the normalized differences based on simulated data for different sample sizes. The data generation process is introduced in equation (ref) in Section (ref). The left panel shows that the mean of the normalized differences increases with the sample size, while the variance converges to a constant. The right panel shows the opposite: the mean converges to a constant, but the variance decreases to zero. This illustrates that, unlike Figure (ref), there is no appropriate normalization that yields a non-degenerate asymptotic distribution. Therefore, unlike fw2016 and wz2021, the order of bias and variance cannot be exactly balanced. However, as shown in Theorem (ref), we can construct an estimator with a correctly centered non-degenerate asymptotic distribution.
Properly handling the incidental parameter Hessian, which enters asymptotic expansions through its inverse, is key to the strategy of fw2016 and therefore also to wz2021. It is important to bound certain quantities in asymptotic expansions and to ensure that the asymptotic bias can be decoupled into separate bias components. Figure (ref) shows the structure of the incidental parameter Hessians without constraints, $\partial_{\phi \phi^{\prime}} L_{u}(\beta, \phi)$, from fw2016 (left panel) and from this paper (right panel), for a data set with $N = T = 5$. For ease of exposition, we look at the unconstrained incidental parameter Hessians, as the difference is already apparent here.
We do not additionally show the Hessian of wz2021, as their proof strategy is fundamentally based on fw2016. The left Hessian is of dimension $(N + T) \times (N + T)$ and the right Hessian is of dimension $(2NT + N^2) \times (2NT + N^2)$. Higher order values (red dots) are located on the diagonals of the matrices, while all non-zero off-diagonal values are of lower order (green dots). Although the right Hessian has a much higher dimension than the left Hessian, both matrices have asymptotically the same number of non-zero off-diagonal elements. The main difference between our Hessian and the one derived by fw2016 is that our Hessian has sparse off-diagonal blocks, while theirs has dense off-diagonal blocks. This sparsity pattern is due to the overlapping fixed effects in our model specification. For example, $\partial_{\alpha_{it} \alpha_{i^{\prime} t^{\prime}}} L_{u}(\beta, \phi) = \sum_{j = 1}^{N} \mu^{\langle 1 \rangle}_{ijt} / \sqrt{NT}$ if $i = i^{\prime}$ and $t = t^{\prime}$ and zero otherwise, because $\alpha_{it}$ and $\alpha_{i^{\prime} t^{\prime}}$ only enter in the same linear index if $i = i^{\prime}$ and $t = t^{\prime}$. This leads to the $NT \times NT$ diagonal block, $\partial_{\alpha \alpha^{\prime}} L_{u}(\beta, \phi)$, with elements of order $N / \sqrt{NT}$. The two other diagonal blocks, $\partial_{\gamma \gamma^{\prime}} L_{u}(\beta, \phi)$ and $\partial_{\rho \rho^{\prime}} L_{u}(\beta, \phi)$, follow analogously. Additionally, $\partial_{\alpha_{it} \gamma_{j t^{\prime}}} L_{u}(\beta, \phi) = \mu^{\langle 1 \rangle}_{ijt} / \sqrt{NT}$ if $t = t^{\prime}$ and zero otherwise, because $\alpha_{it}$ and $\gamma_{j t^{\prime}}$ only enter in the same linear index if $t = t^{\prime}$. This results in the sparse $NT \times NT$ off-diagonal block $\partial_{\alpha \gamma^{\prime}} L_{u}(\beta, \phi)$ with $NT$ elements of order $1 / \sqrt{NT}$. The other off-diagonal blocks, $\partial_{\alpha \rho^{\prime}} L_{u}(\beta, \phi)$, $\partial_{\gamma \alpha^{\prime}} L_{u}(\beta, \phi)$, $\partial_{\gamma \rho^{\prime}} L_{u}(\beta, \phi)$, $\partial_{\rho \alpha^{\prime}} L_{u}(\beta, \phi)$, and $\partial_{\rho \gamma^{\prime}} L_{u}(\beta, \phi)$, follow analogously. Intuitively, although a three-dimensional panel is much larger than a classical panel, the number of observations that can be used to estimate the incidental parameters is asymptotically the same as in a classical panel with individual and time effects. Thus, the increased sample size does not improve the convergence rates of the corresponding estimators, which are still $\sqrt{N}$ or $\sqrt{T}$. Importantly, the sparsity is also reflected in the inverse of the incidental parameter Hessian. Thus, properly handling this sparsity is one of the main challenges in deriving our results (see Appendixes (ref) and (ref) for further details).
In this section, we conduct simulation experiments to study the finite sample behaviour of the uncorrected and debiased maximum likelihood estimators of the model parameters defined in (ref) and (ref), respectively. We analyze biases and the reliability of the derived asymptotic distributions for inference. In particular, we consider the following statistics for our analysis: relative bias in percent, bias relative to standard deviation, and coverage rates of confidence intervals with 95% nominal level. We adapt the static data generating process of fw2016 to logit models for bilateral panels with three sets of overlapping unobserved effects,
where $i, j = \{1,\ldots, N\}$, $t = \{1, \ldots, T\}$, $\alpha_{it}, \gamma_{jt}, \rho_{ij} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 24)$, $u_{ijt} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{U}}(0, 1)$, $v_{ijt} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1 / 2)$, and $x_{ij0} \sim \operatorname{\text{iid.}\;} \operatorname{\mathcal{N}}(0, 1)$. We set $\beta = 1$ and generate data sets with $N \in \{50, 75, \dots, 225, 250 \}$ senders and receivers observed for $T(N) = N / 5 $ time periods. Our study design ensures that $N$ and $T$ grow at a constant rate and is therefore in line with our asymptotic analysis. All results presented are based on $5{,}000$ simulated samples for each $N$.
The left panel of Table (ref) shows the simulation results for the uncorrected estimator. For the smallest sample size, (50, 10), the relative bias is substantial at 18.465%, but decreases steadily with increasing sample size. This is as expected, since the theory predicts that the bias is of order $1 / N + 1 / T$ and should therefore decrease as the panel dimensions increase. For the largest sample size, (250, 50), the relative bias reduces to 3.113%. Although the bias may seem small, it is still large relative to the dispersion of the estimator. This can be seen from the second column, which shows the ratio of bias to standard deviation. More precisely, the ratio actually increases with the sample size, i.e.\ the bias problem gets worse in relative terms as the sample size increases. This bias problem is accordingly reflected in the zero coverage rates shown in the third column. Thus, as predicted by our asymptotic theory, the uncorrected estimator for our model, (ref), exhibits a more severe form of asymptotic bias problem than, for example, the uncorrected estimators in fw2016 and wz2021. The right panel of Table (ref) shows the simulation results for the bias-corrected estimator. If we compare the bias-corrected and the uncorrected estimator, we find that the former outperforms the latter in every metric in every sample. For example, even for the smallest sample size, (50, 10), the bias of 18.465% is reduced to less than 1% and coverage rates are improved from zero to the desired nominal level of 95%. The same applies to all other analyzed sample sizes. Overall, the simulation experiments add numerical evidence that our asymptotic results provide a reasonable approximation for samples with sufficiently large $N$ and $T$.
In the following, we apply the uncorrected and our debiased estimator to real data, using an example from international trade.
To construct a panel data set on bilateral relationships, as described in Example (ref), we combine two data sources. The first data source is the CEPII Gravity Database, provided by ccm2022.\footnote{\url{http://www.cepii.fr/CEPII/en/bdd_modele/bdd_modele_item.asp?id=8}} This database provides information on bilateral trade flows between countries over time, from different sources such as UNSD's Comtrade, IMF DOTS, or CEPII's BACI, as well as other trade cost variables that are frequently used for gravity estimation. The second data source is the Regional Trade Agreements Database, provided by el2008.\footnote{\url{https://www.ewf.uni-bayreuth.de/de/forschung/RTA-daten/index.html}} This database contains additional information about regional trade agreements (RTA), allowing us to distinguish between different but not mutually exclusive types, such as customs unions (CU), free trade agreements (FTA), partial scope agreements (PSA), or economic integration agreements (EIA). Because we use trade flows from CEPII's BACI, which are only available from 1996, and restrict ourselves to the most recent year before the COVID-19 pandemic, our final sample consists of $N = 237$ countries observed between 1996 and 2019 (i.e.\ $T = 24$ years). After removing self-trade and incomplete observations, we are left with an unbalanced panel of $n = 1{,}306{,}232$ observations.
We estimate the following binary logit model,
where $\text{trade}_{ijt}$ is the trade flow from exporting country $i$ to importing country $j$ at time $t$,
is a set of RTA-type indicator variables, $\beta = (\beta_{1}, \ldots, \beta_{4})$ are the corresponding model parameters, $\alpha_{it}$, $\gamma_{jt}$, and $\rho_{ij}$ are three sets of fixed effects accounting for different sources of unobserved heterogeneity (e.g.\ market sizes, multilateral resistance, or other time-invariant trade costs), and $\epsilon_{ijt}$ is an idiosyncratic error term. We lag the RTA-type indicator variables by one period to account for the time it takes for firms to adjust to changes in trade agreements.
Table (ref) presents uncorrected and debiased estimation results for model (ref).
In addition to the model parameter estimates in panel A, we also report the odds ratios (or relative risks) in panel B. Odds ratios are calculated as $\exp(\hat{\beta}_{k})$ for each $k \in \{1, \ldots, 4\}$ and are a useful metric for interpreting the results of logit models. Unlike partial effects, which are another useful metric, odds ratios only depend on the model parameter estimates and therefore do not require further theoretical investigation. We are primarily interested in analyzing the differences between inferences drawn from the uncorrected and debiased estimators. Therefore, we also investigate the corresponding test statistics for typical two-sided hypothesis tests: $\mathbb{H}_{0} \colon \beta_{k} = 0$ for panel A and $\mathbb{H}_{0} \colon \exp(\beta_{k}) = 1$ for panel B, for each $k \in \{1, \ldots, 4\}$. Analyzing panel A, we find that debiasing the estimates substantially reduces the magnitude of the model parameter estimates. Relative to the corresponding standard errors, the reductions range between 0.3 and 1 times the standard error. The debiasing of the estimates also results in lower test statistics. For example, the estimate for EIA becomes insignificant at the 5% level after correcting for the bias. As the odds ratios are just a function of the estimated model parameters, the findings from panel A also carry over to panel B. For example, the uncorrected estimate suggests that forming a free trade agreement increases the probability to trade by 43.4%. However, after debiasing the estimate, we find that the increase is reduced to 37%, which is a 6.4 percentage point reduction. Similarly, the uncorrected estimate suggests that forming a partial scope agreement reduces the probability to trade by 45.7%. After debiasing, the decrease is reduced to 36.2%, which is a 9.5 percentage point reduction.
The empirical example illustrates that, although the panel data set is quite large, with $N = 237$ countries observed for $T = 24$ years, debiasing the estimates significantly impacts the results and the inferences drawn.
We studied the asymptotic behavior of fixed effects estimators for logit models with three additive and overlapping unobserved effects in three-dimensional panels, under asymptotic sequences where all three panel dimensions grow large. To address the asymptotic bias problem of the uncorrected estimator, we proposed a debiasing procedure. The inference problem we identify is more severe than in previous studies, highlighting the need for further research on the properties of fixed effects estimators for nonlinear models with multiple unobserved effects in multi-dimensional panels. Therefore, empirical researchers should be aware of the potential pitfalls of these estimators before using them in practice.
Several interesting topics remain for future research. For instance, our results could be extended to average partial effects and other (potentially dynamic) nonlinear models, as well as to panels with more than three dimensions. Additionally, it could be useful to derive fixed-$T$ consistent fixed effects estimators, as not every panel spans a sufficiently long time period. We plan to explore some of these topics in future work.
{\LARGEAppendix}