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.
95,328 characters · 16 sections · 79 citation commands
Nonlinear Factor Models for Network and Panel Data
\abstract{ Factor structures or interactive effects are convenient devices to incorporate latent variables in panel data models. We consider fixed effect estimation of nonlinear panel single-index models with factor structures in the unobservables, which include logit, probit, ordered probit and Poisson specifications. We establish that fixed effect estimators of model parameters and average partial effects have normal distributions when the two dimensions of the panel grow large, but might suffer from incidental parameter bias. We show how models with factor structures can also be applied to capture important features of network data such as reciprocity, degree heterogeneity, homophily in latent variables, and clustering. We illustrate this applicability with an empirical example to the estimation of a gravity equation of international trade between countries using a Poisson model with multiple factors.
{ Keywords:} Panel data, network data, interactive fixed effects, factor models, bias correction, incidental parameter problem, gravity equation
{ JEL:} C13, C23.
Factor structures or interactive effects are convenient devices to incorporate latent variables in panel data models. They are commonly used to capture aggregate shocks that might have heterogeneous impacts on the agents in macroeconomic models, and multidimensional individual heterogeneity that might have time varying effects in microeconomic models. More generally, the inclusion of these structures serves to account for dependences along the cross-section and time series dimensions in a parsimonious fashion. While methods for linear factor models are well-established, there are very few studies that develop methods for nonlinear factor models. (We provide a literature review at the end of this section.) Nonlinear models are commonly used when the outcome variable is discrete or has a limited support. In this paper we introduce factor structures in single-index nonlinear specifications such as the logit, probit, ordered probit and Poisson models.
The model that we consider is semiparametric. It includes an outcome, strictly exogenous covariates, and a fixed number of factors and factor loadings. The parametric part is the distribution of the outcome conditional on the covariates, factors and loadings, which is specified up to a finite dimensional parameter. The nonparametric part is the distribution of the factors and loadings conditional on the covariates. In other words, our model is of the “fixed effects” type because we do not impose any restriction on the relationship between the observed covariates and the unobserved factors and loadings. This flexibility allows us to capture features of economic behavior more realistically, but poses important challenges to estimation and inference. The objects of interest are the model parameter and average partial effects (APEs), which are averages of functions of the data, parameter, factors and loadings. The APEs measure the effect of covariates on moments of the outcome conditional on the covariates, factors and loadings. We consider a fixed effects estimation approach that treats the factors and loadings as parameters to be estimated. As it is well-known in the panel data literature, the resulting estimators generally suffer from the incidental parameter problem coming from the high-dimensionality of the estimated parameter NeymanScott1948.
We derive asymptotic theory for our estimators of the model parameter and APEs under sequences where the two dimensions of the panel pass to infinity with the sample size. Even establishing consistency is complicated in our setting because the dimension of the estimated parameters increases with the sample size. We develop a new proof of consistency that relies on concavity of the log-likelihood function on a single-index that captures the dependence on covariates, parameter, factors and loadings. However, unlike FW16, we need to deal with the complication that our log-likelihood function is not concave in all the estimated parameters because the factors and loadings enter multiplicatively in the index. We also establish that our estimators are normally distributed in large samples, but might have biases of the same order as their standard deviations. For example, we find that the estimator of the model parameter is asymptotically unbiased in the Poisson model, but is biased in logit and probit models. Following the recent panel data literature, we develop analytical and split-sample corrections for the case where the estimator has asymptotic bias. One specific feature of our estimator is that the bias depends on the number of factors. In particular, we show that the bias grows proportionally with the number of factors in examples.
We discuss implementation details of our methods including the computation of the estimator and selection of the number of factors. Thus, we propose an EM-type algorithm based on C14 and a concrete proposal to estimate the number of factors based on the eigenvalue ratio test of AhnHorenstein2013. The estimator of the number of factors requires to specify an upper bound for the number of factors, but does not rely on any arbitrary choice of penalty function or other tuning parameter. We do not provide asymptotic theory for this estimator, but show that it performs well in numerical simulations. Formally deriving the theory is rather challenging, because it requires to study the asymptotic properties of the initial fixed effects estimators of the parameters and factor structure obtained from a specification with too many factors, which is a difficult problem even in linear panel factor model (MoonWeidner2015). We leave this analysis to future research.
We also introduce factor structures as practical tools to model network data. We show how the inclusion of latent factors is useful to incorporate important features of the network such as reciprocity, degree heterogeneity, homophily on latent variables, and clustering S11, G15. We focus on directed networks with unweighted and weighted outcomes. These cover binary response models for network formation where the outcome is an indicator for the existence of a link between sender and receiver, and count data models for network flows where the outcome is a measure of the volume of flow between sender and receiver. As we shall discuss, our factor model provides a parsimonious reduced-form specification that captures the important network features mentioned above. The statistical treatment of the network factor model is identical to the panel factor model after noticing that a network is isomorphic to a panel after labeling the senders as individuals and the receivers as time periods.
We illustrate the use of the factor structure in network data with an application to gravity equations of trade between countries. We estimate a Poisson model where the outcome is the volume of trade and the covariates include typical gravity variables such as the distance between the countries or whether the country pair belongs to a currency union or a free trade area. The unobserved factors and loadings serve to account for scale and multilateral resistance effects, unobserved partnerships, presence of multinational firms, and differences in natural resources or industrial composition. We find that accounting for these multiple unobserved factors changes the effects of the gravity variables, making all of them to have the expected signs while keeping most of them to be statistically significant.
\paragraph{Literature review: } This paper contributes to the econometric panel data and network data literatures. Regarding the panel data literature, our statistical analysis relies on the recent developments in fixed effects methods. We refer to ARE for a recent review on fixed effects estimation of nonlinear panel models with additive individual and time effects, and to BP16 for a recent review on fixed effects estimation of linear factor or interactive effects panel models. Since the first draft of this paper appeared in CFW14, BL17 and AB16 have considered special cases of nonlinear factor models. BL17 analyzed a probit model using the common correlated random effects approach of Pesaran2006, and AB16 a logit model using a Bayesian approach with data augmentation. Our analysis is different in the modeling assumptions and estimation method.\footnote{We refer to BL17 and AB16 for more detailed comparisons with our analysis.} The most closely related work is FaWang2018. This paper derives the asymptotic distribution of the estimators of the factors and loadings in non-linear single index models without covariates. By contrast, we focus on covariate coefficients and average partial effects and treat the factors and loadings as nuisance parameters. Accordingly, we view our results as complementary to the results in FaWang2018.
In terms of the network literature, our paper is related to the recent work on the application of panel fixed effects methods to network data including FW16, Yan2016statistical, Stata2017, Dzemski2017, Graham17, and yan2018. These papers account for degree heterogeneity by including additive unobserved sender and receiver effects. Additive effects, however, do not capture other network features such as homophily in latent factors and clustering. Graham2016 considered a binary response model of network formation with all these features plus state dependence, for the case where the network is observed at multiple time periods. Compared to Graham2016, our method can capture all these features, except for state dependence, applies to ordered and count outcomes in addition to binary outcomes, and only requires observing the network at one time period. A stream of the statistic literature has considered nonlinear factor network models using a random effects approach including HRH02, hoff05, KHRH09, and HRT07. Unlike the fixed effects approach that we adopt, the random effects approach assumes independence between covariates and factors and between covariates and loadings. This assumption is regarded as implausible for most economic applications where the loadings reflect unobserved individual heterogeneity and some of the covariates are individual choice variables. There is also a recent econometric literature on structural models of strategic network formation where the main focus is on identification. We refer to depaula2017 for an excellent up-to-date review on this topic. The focus of our paper is on estimation and inference.
Finally, there is an extensive literature in international economics on the estimation of the gravity equation including Harrigan1994, EatonKortum2001, AndersonWincoop2003, SantosSilvaTenreyro2006, Helpman01052008, Charbonneau2011 and jochmans2017two. We refer to HEAD2014 for a recent review on this literature. These papers estimate models with additive unobserved sender and receiver country effects to account for scale or multilateral resistence effects. Our innovation to this literature is the inclusion of multiple unobserved factors to account for not only scale effects, but also unobserved partnerships, and homophily induced by differences in natural resources, industrial composition or other country characteristics.
To sum-up, our paper makes the following contributions. First, we derive asymptotic theory for fixed effects estimators of model parameters and APEs in a class of nonlinear single-index factor models that include logit, probit, ordered probit and Poisson models. Second, we provide bias corrections for fixed effects estimators of model parameters and APEs. Third, we propose an estimator of the number of factors in nonlinear single-index models with factor structure. Fourth, we bring in the factor structure to model important features of network data such as reciprocity, degree heterogeneity, homophily in latent factors and clustering in a reduced form fashion. Fifth, we apply our methods to the estimation of a gravity equation of trade between countries and confirm the importance of the gravity variables even after conditioning on multiple unobserved latent factors.
\paragraph{Outline: } In Section (ref), we introduce the model and estimators. Section (ref) discusses the statistical issues in the estimation and inference of factor models with a simple example. Section (ref) derives asymptotic theory for our estimators. Section (ref) provides implementation details for the estimators of the parameters and number of factors. Section (ref) describes the results of the empirical application to the gravity equation and a calibrated simulation. The proofs of the main results and other technical details are given in the Appendix.
We observe the data $\{(Y_{ij}, X_{ij}) : (i,j) \in \mathcal{D} \}$, where $Y_{ij}$ is a scalar outcome variable and $X_{ij}$ is a $d_x$-dimensional vector of covariates. The subscripts $i$ and $j$ index individuals and time periods in traditional panels, but they might index different dimensions in other data structures such as network data. In our empirical application, for example, we use country trade network data where $Y_{ij}$ is the volume of trade between country $i$ and country $j$, and $X_{ij}$ includes gravity variables such as the distance between country $i$ and country $j$. Both $i$ and $j$ index countries as exporters and importers respectively. The set $\mathcal{D}$ contains the indexes of the units that are observed. It is a subset of the set of all possible pairs $\mathcal{D}_0 := \{(i,j) : i = 1,\dots,I; j = 1, \dots, J \}$, where $I$ and $J$ are the dimensions of the data set. We introduce $\mathcal{D}$ to allow for missing data that are common in panel and network applications. For example, in the trade application $I=J$ and $\mathcal{D} = \mathcal{D}_0 \setminus \{(i,i) : i = 1,\dots,I \}$ because we do not observe trade of a country with itself. We denote the total number of observations by $n$, i.e. $n = |\mathcal{D}|$.
We assume that the outcome is generated by
where $f$ is a known density function with respect to some dominating measure, $\beta$ is $d_{x}$-dimensional parameter vector, and $\alpha_i$ and $\gamma_j$ are $R$-vectors of unobserved effects. We collect these effects in the $I \times R$ matrix $\alpha = (\alpha_1,\ldots, \alpha_I)'$, and the $J \times R$ matrix $\gamma = (\gamma_1,\ldots, \gamma_J)'$, which are further stacked in the $R(I+J)$-vector $\phi_n = ({\rm vec}(\alpha)', {\rm vec}(\gamma)')'$. We make explicit in $\phi_n$ that the number of unobserved effects changes with the sample size because it will have important effects on the asymptotic theory. We assume that the dimension of the unobserved effects $R$ is known, and provide a practical method to estimate $R$ in Section (ref). The effects $\alpha_i$ and $\gamma_j$ are unobserved factors and factor loadings. In panel data they represent individual and time effects that in economic applications capture individual heterogeneity and aggregate shocks, respectively. In network data $\alpha_i$ and $\gamma_j$ represent unobserved characteristics of senders and receivers that affect the network flow. The model is semiparametric because we do not specify the distribution of the unobserved effects nor their relationship with the covariates. This flexibility is important for economic applications where some of the covariates are choice variables with values determined in part by the unobserved effects. The conditional distribution $f$ represents the parametric part of the model.
The model has a single-index specification because the covariates and unobserved effects enter $f$ through the index $ z_{ij} = X_{ij}'\beta + \alpha_i' \gamma_j$. The parameter $\beta$ is a quantity of interest because it measures the effect of the covariates on the distribution of the outcome controlling for the unobserved effects. For example, in network data $\beta$ can measure homophily in an observable characteristic $W$ if $X_{ij}$ includes $(W_i - W_j)^2$ as one of its components. The unobserved effects have a factor or interactive structure because they enter the index $z_{ij}$ multiplicatively through $\pi_{ij} = \alpha_i' \, \gamma_j$. The standard additive structure $\alpha_{1i} + \gamma_{1j}$ can be seen as a special case of the factor structure with $R=2$, $\alpha_i = (\alpha_{1i},1)'$, and $\gamma_j = (1, \gamma_{1j})'$. More generally, in panel data applications the factor structure allows one to incorporate multiple aggregate shocks $\gamma_t$ with heterogeneous effects across agents $\alpha_i$, or multidimensional individual heterogeneity $\alpha_i$ with time-varying returns $\gamma_t$. For example, we can have productivity and monetary shocks with heterogeneous effects across industries, or multiple dimensions of individual ability and skills with time-varying returns in the labor market.
One of the contributions of the paper is to introduce factor structures to network data. In this case the factor structure serves to capture important network features in an unspecified or reduced-form fashion. For example, degree heterogeneity can be captured with the additive structure $\alpha_{1i} + \gamma_{1j}$ mentioned above, and reciprocity by allowing $Y_{ij}$ to be arbitrarily related to $Y_{ji}$ even after conditioning on the covariates and unobserved effects. Another important feature is homophily in latent factors, in addition to the homophily on observed factors captured by $X_{ij}$. Assume that there is a latent factor $\xi_i$ such that the flow between $i$ and $j$ increases or decreases with the distance between $\xi_i$ and $\xi_j$ as measured by $(\xi_i - \xi_j)^2$. This type of homophily can also be captured by a factor structure with $R=3$, $\alpha_i=(\xi_i^2,1,-2\xi_i)'$ and $\gamma_j = (1,\xi_j^2,\xi_j)$. The factor structure can also account for clustering or transitivity of links due to latent factors. Assume that there is a cluster of individuals such that there are more flows within the cluster. This would be captured by a factor structure with $R=1$, $\alpha_i = \xi_i I_i$ and $\gamma_j = \chi_j I_j$, where $\xi_i$ and $\chi_j$ are positive cluster effects on the sender and receiver, and $I_i$ is an indicator for cluster membership. The factor structure can also account for combinations of these network features. Indeed, one of its advantages is that the researcher has the flexibility of specifying some features and leaving other features unspecified. For example, in the trade application we use a specification that includes additive effects to account explicitly for degree heterogeneity and multiple interactive effects to account for the possibility of having homophily in latent factors and clustering without explicitly modelling any of them.
We consider three running examples throughout the analysis:
In addition to the model parameter $\beta$, we might be interested in average partial effects (APEs). These effects are averages of the data, parameters and unobserved effects. They measure the effect of the covariates on moments of the distribution of the outcome conditional on the covariates and unobserved effects. The leading case is the conditional expectation, $$ \mathbb{E} [Y_{ij} \mid X_{ij}, \alpha_i, \gamma_j, \beta] = \int y f(y \mid X_{ij}'\beta + \pi_{ij} ) dy, $$ where the partial effects are differences or derivatives of this expression with respect to the components of $X_{ij}$. We denote generically the partial effects by $\Delta(Y_{ij},X_{ij}, \beta, \alpha'_i \gamma_j) = \Delta_{ij}(\beta, \alpha_i' \gamma_j)$, where the restriction that they depend on $\alpha_i$ and $\gamma_j$ through $\pi_{ij}$ is natural given the model for the conditional density of $Y_{ij}$. We allow the partial effect to depend on $Y_{ij}$ to cover scale and other parameters not included in the single-index. The APE is
Example (ref) (Linear model). The variance $\sigma^2$ in the linear model can be expressed as an APE with
}
Example (ref) (Binary response model). If $X_{ij,k}$, the $k$th element of $X_{ij}$, is binary, its partial effect on the conditional probability of $Y_{ij}$ is
where $\beta_k$ is the $k$th element of $\beta$, and $X_{ij,-k}$ and $\beta_{-k}$ include all elements of $X_{ij}$ and $\beta$ except for the $k$th element. If $X_{ij,k}$ is continuous and $F$ is differentiable, the partial effect of $X_{ij,k}$ on the conditional probability of $Y_{ij}$ is
}
Example (ref) (Count response model). If $X_{ij,k}$, the $k$th element of $X_{ij}$, is binary, its partial effect on the conditional probability of $Y_{ij}$ in the Poisson model is
where $\beta_k$ is the $k$th element of $\beta$, and $X_{ij,-k}$ and $\beta_{-k}$ include all elements of $X_{ij}$ and $\beta$ except for the $k$th element. If $X_{ij,k}$ is continuous, the partial effect of $X_{ij,k}$ on the conditional expectation of $Y_{ij}$ is
}
We adopt a fixed effects approach and treat the unobserved effects $\phi_n$ as a vector of nuisance parameters to be estimated. Let
be the conditional log-likelihood function of the data constructed from the parametric part of the model. The fixed effects estimator is
This problem has a unique solution with probability one for $\beta$ under the assumption that $z \mapsto \log f(\cdot \mid z)$ is concave. This assumption holds for all the cases that we consider including logit, probit, ordered probit and Poisson models. The solution for $\phi_n$ is only unique up to normalization -- see Remark (ref) below. Obtaining the solution to (ref) can be computationally challenging because the objective function is not concave in the parameter $\phi_n$ and the high-dimensionality of the parameter space. In Section (ref) we provide an iterative method based on C14 to obtain the estimates. This method performs well in simulations.
Let $\widehat \phi_n = ({\rm vec}(\widehat \alpha)', {\rm vec}(\widehat \gamma)')' $, where $\widehat \alpha$ and $\widehat \gamma$ correspond to the components $\alpha$ and $\gamma$ such that $\widehat \alpha = (\widehat \alpha_1, \ldots, \widehat \alpha_I)'$ and $\widehat \gamma = (\widehat \gamma_1, \ldots, \widehat \gamma_J)'$. Plugging the estimator of $(\beta,\phi_n)$ in (ref) yields the estimator of the APE,
In Section (ref), we show that $\widehat \beta$ and $\widehat \delta$ are consistent and normally distributed in large samples, but might have incidental parameter bias because the dimension of the nuisance parameter $\phi_n$ grows with the sample size NeymanScott1948.
We illustrate the statistical issues that arise in the estimation of factor models with a simple example. This example is analytically tractable and might be of practical interest as it provides an estimator of the variance of a random variable in network and panel data allowing for flexible patterns of dependence. The analysis in this section is mainly heuristic leaving technical details such as the derivation of the orders of some remainder terms in the asymptotic expansions for Section (ref).
Consider a version of Example (ref) without covariates where $Y_{ij} \mid \phi_n \sim \mathcal{N}(\alpha'_i \gamma_j, \sigma^2)$. Assume that the observations $Y_{ij}$ are independent over $i$ and $j$, and that there is no missing data, i.e. ${\cal D}= \mathcal{D}_0$. The quantity of interest is the scale parameter $\sigma^2$, which can be treated as an APE. This is a linear factor model where $\widehat \phi_n$ can be obtained using the principal component algorithm of Bai:2009p3321. Then, the plug-in estimator of $\sigma^2$ is
To analyze the properties of $\widehat \sigma^2$, it is useful to consider an asymptotic expansion of $\widehat \alpha_i' \widehat \gamma_j$ around $\alpha_i' \gamma_j$ as $I,J \to \infty$. This yields
where $\approx$ means equal up to terms of lower order. Plugging this expansion in (ref) shows that $\widehat \sigma^2$ behaves asymptotically as a sample variance with $R(I+J)$ estimated fixed effects corresponding to the $\widehat \alpha_i$'s and $\widehat \gamma_t$'s. Then, standard degrees of freedom calculations give
which shows that $\widehat \sigma^2$ has an incidental parameter bias that grows proportionally to the number of factors $R$. The order of the bias corresponds to the number of estimated parameters, $R(I+J)$, divided by the number of observations, $IJ$, as predicted by the general formula in ARE for fixed effects estimators. We show in numerical examples that this expression produces a very accurate approximation to the bias even for small sample sizes.
We carry out 50,000 simulations with $\sigma^2 = 1$, and $\alpha_i$ and $\gamma_j$ drawn independently from multivariate normal distributions with mean zero and covariance function $\mathbb{I}_R$, the identity matrix of order $R$. Table (ref) compares the bias of $\widehat \sigma^2$ with the asymptotic approximation (ref) in datasets with $I,J \in \{10, 25, 50\},$ and $R \in \{1,2,3\}$. We only report the results for $J \leq I$ since all the expressions are symmetric in $I$ and $J$. Comparing the two rows in each panel of the table, we find that the asymptotic bias provides a very accurate approximation to the finite-sample bias of the estimator for all the sample sizes and numbers of factors.
The bias of $\widehat \sigma^2$ can be removed using analytical and split-sample methods. Thus, an analytical bias corrected estimator can be formed as $$ \widetilde \sigma_{\rm ABC}^2 = \frac{IJ}{(I-R)(J-R)} \widehat \sigma^2. $$ A split-sample bias corrected estimator can be formed as $$ \widetilde \sigma_{\rm SBC}^2 = 3 \widehat{\sigma}^2 - \bar{\sigma}^2_{I,J/2} - \bar{\sigma}^2_{I/2,J}, $$ where $\bar{\sigma}^2_{I,J/2}$ is the average of the estimators in the half-panels $ \{(i,j) : i = 1,\dots,I; j = 1, \dots, \lceil J/2 \rceil \}$ and $ \{(i,j) : i = 1,\dots,I; j = \lfloor J/2 + 1\rfloor, \dots, J \}$, and $\bar{\sigma}^2_{I/2,J}$ is the average of the estimators in the half-panels $ \{(i,j) : i = 1,\dots,\lceil I/2 \rceil; j = 1, \dots, J \}$ and $ \{(i,j) : i = \lfloor I/2 + 1\rfloor,\dots,I; j = 1, \dots, J \}$, where $\lceil \cdot \rceil$ and $\lfloor \cdot \rfloor$ are the ceil and floor functions. As in nonlinear panel data, we expect these corrections to remove most of the bias of the estimator without increasing dispersion. Moreover, constructing confidence intervals around the corrected estimators should help bring coverage probabilities close to their nominal levels. We confirm these predictions in a numerical simulation.
Table (ref) reports the bias, standard deviation and RMSE of the uncorrected and bias corrected estimators, together with coverage probabilities of 95% confidence interval constructed around them. The results are based on 50,000 simulations of datasets generated as in Table (ref) with $I, J \in \{10, 25, 50\}$, and $R = 3$. The confidence intervals around the estimator $\widetilde \sigma^2 \in \{ \widehat \sigma^2, \widetilde \sigma_{\rm ABC}^2, \widetilde \sigma_{\rm SBC}^2\} $ are constructed as $\widetilde \sigma^2(1 \pm 1.96 \sqrt{2/(IJ)})$, where we use that the asymptotic variance of all the estimators is $2\sigma^4/(IJ)$. We find that the corrections offer huge improvements in terms of bias reduction and coverage of the confidence intervals. The corrections increase the dispersion for small sample sizes, but always reduce the RMSE. In this case the analytical correction slightly outperforms the split-sample correction.
We derive the asymptotic distribution of the estimators of the model parameter and APEs under sequences where $I$ and $J$ grow with the sample size at the same rate. We focus on these sequences because they are the only ones that deliver a non-degenerate limit distribution. Moreover, they are very natural choices for network data where $I = J$. Throughout this section, all the stochastic statements are conditional on the realization of the unobserved effects $\phi_n$ and should therefore be qualified with almost surely. We shall omit this qualifier to lighten the notation.
We consider single-index models with strictly exogenous covariates and unobserved effects that enter the density of the outcome through $z_{ij} = X_{ij}'\beta + \pi_{ij}$, where $\pi_{ij} = \alpha_i'\gamma_j$. These models cover the linear, probit and Poisson specifications of Examples (ref)--(ref). We focus on strictly exogenous covariates because for some data structures of interest such as network data there is no natural ordering of the observations. The results can be extended to predetermined covariates when one of the dimensions is time, see the earlier version of the paper CFW14. Let
be the conditional log-likelihood coming from the parametric part of the model. We denote the derivatives of $z \mapsto \ell_{ij}(z)$ by $\partial_{z^q} \ell_{ij}(z) := \partial^q \ell_{ij}(z)/\partial z^q$, $q = 1,2, \ldots$. Let $\beta^0$, $\alpha_i^0$, $\gamma_j^0$, and $\pi_{ij}^0 = \alpha_i^{0 \prime} \gamma_j^0$ denote the values of $\beta$, $\alpha_i$, $\gamma_j$, and $\pi_{ij}$ that generated the data. We drop the argument $z_{ij}$ when the derivatives are evaluated at the true value of the index $z^0_{ij} := X_{ij}'\beta^0 + \pi_{ij}^0$, i.e., $\partial_{z^q} \ell_{ij} := \partial_{z^q}\ell_{ij}(z^0_{ij})$. Let $\boldsymbol{X}=\{X_{ij} : (i,j) \in {\cal D}\}$, $\alpha^0 = (\alpha_1^0,\ldots, \alpha_I^0)'$, and $\gamma^0 = (\gamma_1^0,\ldots, \gamma_J^0)'$ .
We make the following assumptions:
The two cases considered in Assumption (ref)$(i)$ are designed for different data structures. Case (b) is more suitable for network data because it allows for reciprocity between the observations $(i,j)$ and $(j,i)$, whereas case (a) is more suitable for panel data where there is no special relationship between these observations. Assumption (ref)$(i)$ also imposes that the number of factors is known. We provide a practical method to choose the number of factors in Section (ref). We also recommend checking the sensitivity to this number by reporting the maximum value of the average log-likelihood and the parameter estimates for multiple values of $R$. We provide an example in the empirical application of Section (ref). Assumption (ref)$(i)-(iii)$ are similar to FW16, so we do not discuss them further here. The concavity condition in Assumption (ref)$(iv)$ holds for the logit, probit, ordered probit and Poisson models. The strong factor and generalized noncollinearity conditions in Assumption (ref)$(v)-(vi)$ were previously imposed in Bai:2009p3321 and MoonWeidner2015,MoonWeidner2017 for linear models with interactive effects. Generalized noncollinearity rules out covariates that do not display variation in the two dimensions of the dataset. BL17 and AB16 impose very similar conditions to Assumption (ref), so we refer to these papers for further discussion.
We introduce more notation that is convenient to simplify the expressions in the asymptotic distribution. Let $\Xi_{ij}$ be a $d_x$-dimensional vector defined by the following population weighted least squares projection for each component of $\mathbb{E}( \partial_{z^2} \ell_{ij} X_{ij})$,
Also define the residual of the projection
Finally, let $\overline{\mathbb{E}} := \operatorname*{plim}_{I,J \to \infty}$, ${\cal D}_i := \{ j \, : \, (i,j) \in {\cal D} \}$ and ${\cal D}_j := \{ i \, : \, (i,j) \in {\cal D} \}$.
The following theorem establishes the asymptotic distribution of $\widehat \beta$ defined in (ref).
Theorem (ref) shows that $\widehat \beta$ is consistent and normally distributed, but can have bias of the same order as its standard deviation. The scaling factors in the expressions for $\overline B_{\infty}$ and $\overline D_{\infty}$ are such that those expressions are of order one, for example, we can express $\overline B_{\infty}$ equivalently as $$ - \overline{\mathbb{E}} \left\{ \frac {1} {I} \sum_{i=1}^I \frac 1 {|{\cal D}_i|} \sum_{j \in {\cal D}_i} \gamma^{0 \, \prime}_j \left[ \frac 1 {|{\cal D}_i|} \sum_{h \in {\cal D}_i} \gamma^{0}_h \gamma^{0 \prime}_h \, \mathbb{E}\left( \partial_{z^2} \ell_{ih} \right) \right]^{-1} \gamma^{0}_j \; \mathbb{E}\left( \partial_{z} \ell_{ij} \partial_{z^2} \ell_{ij} \tilde X_{ij} + \frac 1 2 \partial_{z^3} \ell_{ij} \tilde X_{ij} \right) \right\} , $$ where all sums explicitly appear as part of a sample average. We verify the presence of bias in our running examples.
Example (ref) (Linear model). In this case $$ \ell_{ij}(z) = - \frac{1}{2} \log (2 \pi \sigma^2) - \frac{(Y_{ij} -z_{ij})^2}{2\sigma^2}, $$ so that $\partial_{z} \ell_{ij} = (Y_{ij} -z_{ij}^0)/\sigma^2$, $ \partial_{z^2} \ell_{ij} = -1/\sigma^2$, and $ \partial_{z^3} \ell_{ij} = 0$. Substituting these values in the expressions of the bias of Theorem (ref) yields $ \overline B_{\infty} = \overline D_{\infty} = 0, $ which agrees with the result in Bai:2009p3321 of no asymptotic bias for $\beta$ in homoskedastic linear models with interactive effects and strictly exogenous covariates.
Example (ref) (Binary response model). In this case $$\ell_{ij}(z) = Y_{ij} \log F(z) + (1 - Y_{ij}) \log [1 - F(z)],$$ so that $\partial_{z} \ell_{ij} = H_{ij} (Y_{ij} - F_{ij}),$ $ \partial_{z^2} \ell_{ij} = - H_{ij} \partial F_{ij} + \partial H_{ij} (Y_{ij} - F_{ij})$, and $\partial_{z^3} \ell_{ij} = - H_{ij} \partial^2 F_{ij} - 2 \partial H_{ij} \partial F_{ij} + \partial^2 H_{ij} (Y_{ij} - F_{ij})$, where $H_{ij} = \partial F_{ij} / [F_{ij}(1-F_{ij})], and $ $\partial^{j} G_{ij} := \partial^{j} G(Z)|_{Z =z_{ij}^0}$ for any function $G$ and $j = 0,1,2$. Substituting these values in the expressions of the bias of Theorem (ref) for the probit model yields
The asymptotic bias is therefore a positive definite matrix weighted average of the true parameter value as in the case of the probit model with additive individual and time effects in FW16. The bias grows linearly with the number of factors because
and $\mathbb{E}\left( \partial_{z^2} \ell_{i j} \right)$ and $\mathbb{E}\left(\partial_{z^2} \ell_{ij} \tilde{X}_{ij} \tilde{X}_{ij}' \right) $ are bounded uniformly in $i,j$. }
Example (ref) (Count response model). In this case $$\ell_{ij}(z) = z Y_{ij} - \exp(z) - \log Y_{ij}!,$$ where the symbol $!$ denotes the factorial function, so that $\partial_{z} \ell_{ij} = Y_{ij} - \lambda_{ij}$ and $\partial_{z^2} \ell_{ij} = \partial_{z^3} \ell_{ij}= - \lambda_{ij} $, where $\lambda_{ij} = \exp(z_{ij}^0)$. Substituting these values in the expressions of the bias of Theorem (ref) yields $$\overline B_{\infty} = \overline D_{\infty} = 0,$$ which generalizes the result in FW16 of no asymptotic bias in the Poisson model with strictly exogenous covariates and additive individual and time effects to the Poisson model with strictly exogenous covariates and factor structure.
We use additional assumptions to derive the asymptotic distribution of the estimator of the APEs. They involve smoothness conditions on the partial effect function $(\beta,\pi) \mapsto \Delta_{ij}(\beta,\pi)$ needed to obtain the limit distribution of $\widehat \delta$ from the limit distribution of $(\widehat \beta, \widehat \phi_n)$ via delta method. For a vector of nonnegative integer numbers $v = (v_1, \ldots, v_{d_x})$, let $\partial_{\beta^v}:= \partial^{|v|} / \partial \beta_1^{v_1} \cdots \partial \beta_{d_x}^{v_{d_x}}$ and $|v| = v_1 + \ldots + v_{d_x}$.
It is convenient again to introduce notation to simplify the expressions in the asymptotic distribution. Let $\Psi_{ij}$ be the weighted least squares population projection
We denote the partial derivatives of $(\beta,\pi) \mapsto \Delta_{ij}(\beta,\pi)$ by $\partial_\beta \Delta_{ij}(\beta, \pi) := \partial \Delta_{ij}(\beta,\pi)/\partial \beta$, $\partial_{\beta \beta'} \Delta_{ij}(\beta, \pi) := \partial^2 \Delta_{ij}(\beta,\pi)/(\partial \beta \partial \beta')$, $\partial_{\pi^q} \Delta_{ij}(\beta, \pi) := \partial^q \Delta_{ij}(\beta,\pi)/\partial \pi^q$, $q = 1,2,3,\ldots$. We drop the arguments $\beta$ and $\pi$ when the derivatives are evaluated at the true values $\beta^0$ and $\pi^0_{ij}$, e.g. $\partial_{\pi^q} \Delta_{ij} := \partial_{\pi^q}\Delta_{ij}(\beta^0,\pi^0_{ij})$. We also define $D_{\pi} \Delta_{ij} := \partial_{\pi} \Delta_{ij} - \partial_{z^2} \ell_{ij} \Psi_{ij}$ and $D_{\pi^2} \Delta_{ij} := \partial_{\pi^2} \Delta_{ij} -\partial_{z^3} \ell_{ij} \Psi_{ij}$.
We are now ready to present the asymptotic distribution of $\widehat \delta$ defined in (ref).
Theorem (ref) shows that $\widehat \delta$ is consistent and normally distributed, but can have bias of the same order as its standard deviation. The first two terms of the bias come from the bias of $\widehat \beta$. They drop out when either $\widehat \beta$ does not have bias or the APE is estimated from a bias corrected estimator of $\beta$. We verify the presence of bias in two of the running examples.
Example (ref) (Linear model). In this case $\overline B_{\infty} = \overline D_{\infty} = 0$ and $$ \Delta_{ij}(\beta,\pi) = (Y_{ij} - X_{ij}'\beta - \pi)^2, $$ so that $\partial_z \Delta_{ij} = -2(Y_{ij} - X_{ij}'\beta^0 - \pi_{ij}^0)$ and $\partial_{z^2} \Delta_{ij} = 2$. Substituting these values in the expressions of the bias of Theorem (ref) yields $$ \overline B^\delta_{\infty} = \overline D^\delta_{\infty} = - R \sigma^2, $$ where we use (ref). This result formalizes the analysis in Section (ref)
Example (ref) (Binary response model). Let $\Delta_{ij}(\beta,\pi)$ be as defined in either ((ref)) or ((ref)). Using the notation previously introduced for this example, the expressions of $\overline B_{\infty}^{\delta}$ and $\overline D_{\infty}^{\delta}$ in Theorem (ref) yield
As for the model parameter, these bias terms grow linearly with the number of factors $R$.}
Example (ref) (Count response model). Let $\Delta_{ij}(\beta,\pi)$ be as defined in either (ref) or (ref). In this case $\overline B_{\infty} = \overline D_{\infty} = 0$, and $\partial_z \Delta_{ij} = \partial_{z^2} \Delta_{ij} = \Delta_{ij}$. Substituting these values in the expressions of the bias of Theorem (ref) yields $$ \overline B^\delta_{\infty} = \overline D^\delta_{\infty} = 0, $$ which generalizes the result in FW16 of no asymptotic bias for the estimators of the APEs in the Poisson model with strictly exogenous covariates and additive individual and time effects to the Poisson model with strictly exogenous covariates and factor structure.
Theorems (ref) and (ref) establish that the estimators of the model parameter and APEs have a bias of the same order as their standard deviations in some models. In this section, we describe how to apply recent developments in nonlinear panel data to correct the bias from the estimators. To simplify the notation we assume that there is no missing data.\footnote{ We refer to ARE for a discussion on how to modify the corrections to deal with missing data.} We consider a generic estimator $\widehat \theta$ of the parameter $\theta$, which may correspond to the model parameter or an APE. In this notation, Theorems (ref) and (ref) show that $\widehat \theta$ can have a bias $\mathcal{B}_{\infty} = \overline{\mathbb{E}}[ \mathcal{B}(\beta^0,\phi_n^0)]$ with structure $$ \mathcal{B}(\beta,\phi_n) = \frac{B(\beta,\phi_n)}{J} + \frac{D(\beta,\phi_n)}{I}. $$ The intuition behind this structure is that there are $J$ observations that are informative to estimate each $\alpha_i$ and $I$ observations that are informative to estimate each $\gamma_j$.
An analytical correction based on Hahn:2004p882 and FW16 can be formed as $$ \widetilde \theta_{\rm ABC} = \widehat \theta - \widehat \mathcal{B}, \ \ \ \ \widehat \mathcal{B} = \mathcal{B}(\widehat \beta, \widehat \phi_n). $$ A split-sample correction based on DJ2015 and FW16 can be formed as $$ \widetilde \theta_{\rm SBC} = 3 \widehat \theta - \bar \theta_{I,J/2} - \bar \theta_{I/2,J}, $$ where $\bar{\theta}_{I,J/2}$ is the average of the estimators in the haft-panels $ \{(i,j) : i = 1,\dots,I; j = 1, \dots, \lceil J/2 \rceil \}$ and $ \{(i,j) : i = 1,\dots,I; j = \lfloor J/2 + 1\rfloor, \dots, J \}$, and $\bar{\theta}_{I/2,J}$ is the average of the estimators in the haft-panels $ \{(i,j) : i = 1,\dots,\lceil I/2 \rceil; j = 1, \dots, J \}$ and $ \{(i,j) : i = \lfloor I/2 + 1\rfloor,\dots,I; j = 1, \dots, J \}$, where $\lceil \cdot \rceil$ and $\lfloor \cdot \rfloor$ are the ceil and floor functions. For network data where $I=J$ and the two dimensions of the data index the same entities, Stata2017 proposed the leave-one-out correction $$ \widetilde \theta_{\rm NBC} = I \widehat \theta - (I-1) \bar \theta_{I-1}, \ \ \bar \theta_{I-1} = I^{-1} \sum_{i=1}^I \widehat \theta_{-i}, $$ where $\widehat \theta_{-i}$ is the estimator in the subpanel $ \{(k,j) : k = 1,\dots,I; j = 1, \dots, I, k \neq i, j \neq i \}$, that is, the original panel leaving out the observations corresponding to the entity $i$ as either sender or receiver.
The discussion of bias correction so far is applicable very generally to network and panel models with two-way fixed effects. We now specialize it to our nonlinear models with interactive fixed effects. For the analytic bias correction and for variance estimation we require consistent estimators for the quantities $\overline B_{\infty} $, $ \overline D_{\infty}$, $\overline{W}_{\infty}$, and $ \overline \Sigma_{\infty}$ defined in Theorem (ref). Let $\widehat B$, $ \widehat D$, $\widehat W$ and $\widehat \Sigma$ be the corresponding sample analogs, obtained by simply dropping expectations and plugging in the fixed effect estimators for the true value of the parameters. For example,
where $ \partial_{z^2} \widehat \ell_{ij} = \partial_{z^2} \ell_{ij} \left( X_{ij}' \widehat \beta + \widehat \alpha_i' \, \widehat \gamma_j \right)$, and $\widehat \Xi_{ij} $ is the $d_x$-vector with elements $ \widehat \Xi_{it,k} = \alpha^{\# \, \prime}_{i,k} \widehat \gamma_j + \widehat \alpha^{ \prime}_i \gamma^{\#}_{t,k} $, with $ \alpha^{\# \, \prime}_{i,k} $ and $\gamma^{\#}_{t,k} $ obtained as the solution to
Once those sample analogs are constructed, then the analytic bias correction of $\widehat \beta$ reads $$\widetilde \beta_{\rm ABC} = \widehat \beta - \frac{I} n \, \widehat W^{-1} \widehat B - \frac{J} n \, \widehat W^{-1} \widehat D. $$ Analogously, we can construct sample analogs for $ \overline B^\delta_{\infty} $, $\overline D^\delta_{\infty} $, $ \overline {(D_{\beta} \Delta)}_{\infty}$, defined in Theorem (ref), in order to construct $ \widetilde \delta_{\rm ABC}$. Also, let $ \widehat V^{\delta}$ be the sample analog of $\overline V_{\infty}^{\delta}$.
Theorem (ref) shows that analytic bias correction can be used to obtain estimators of $\beta^0$ and $\delta^0$ that are asymptotically unbiased. It also shows that the simple plug-in estimators of the asymptotic variances are consistent, thus allowing to perform asymptotically valid hypothesis tests and to construct asymptotically valid confidence intervals for $\beta^0$ and $\delta^0$.
Showing that the Jackknife corrected estimators $\widetilde \beta_{\rm JBC}$ and $\widetilde \delta_{\rm JBC}$ have the same asymptotic distribution as $\widetilde \beta_{\rm ABC}$ and $\widetilde \delta_{\rm ABC}$ requires an additional homogeneity assumption, which guarantees that the unconditional distribution of the data is stationary across $i$ and $j$. This assumption ensures that the terms $B$ and $D$ in the bias expansion of $ \widehat \theta$ are the same as in the bias expansions of the half-panel estimates $\bar \theta_{I,J/2}$ and $\bar \theta_{I/2,J}$, so that forming the Jackknife linear combination $\widetilde \theta_{\rm SBC}$ indeed cancels those bias terms. In other words, the data distribution should not systematically differ across the subsamples used for the Jackknife correction DJ2015, FW16.
The derivation of the asymptotic distribution of the leave-one-out correction $\widetilde \theta_{\rm NBC}$ furthermore requires a third-order bias expansion (i.e., up to terms of order $1/I^2$), because in the expression of $\widetilde \theta_{\rm NBC}$ the estimators $\widehat \theta$ and $\bar \theta_{I-1}$ are multiplied by the factors $I$ and $(I-1)$ that grow with the sample size. We have not worked out those higher-order expansion here, but we refer to SunDhaene2017 for an example of higher-order expansions in nonlinear panel models.
We apply the following EM-type algorithm based on C14 to find the solution to the program (ref):
C14 analyzed the convergence guarantees for this algorithm. She showed that the algorithm converges to a local maximum of the log-likelihood. Since the log-likelihood can have multiple local maxima, we recommend to run the algorithm for several initial values and choose the solution that yields the highest value of the log-likelihood.
The problem of estimating the number of factors $R$ has been extensively discussed for linear factor models without covariates, see for example, BaiNg2002,hallin2007generalized,Onatski2010,alessi2010improved,AhnHorenstein2013. These methods can be extended to linear models with covariates, provided that an appropriate preliminary estimator $\widetilde \beta$ of the regression parameters $\beta$ is available that does not require knowing $R$. In this case the existing methods are applied to the residuals $Y_{ij} - X_{ij}' \widetilde \beta$. If there exists an upper bound for the number of factors, $R_{\max} \geq R$, then the preliminary estimator $\widetilde \beta$ is given by the least squares estimator with $R_{\max}$ factors, see MoonWeidner2015. These methods can also be extended to the nonlinear factor models that we consider. For example, the various information criteria in BaiNg2002 are all based on minimizing the sum of squared residuals plus a penalty function, and can be adapted to the likelihood problem in the spirit of classic model selection criteria (AIC, BIC, etc), see AB16 for an example of this approach.\footnote{kss12 proposed an alternative estimator of the number of factors in linear models specially adapted to i.i.d. errors.} It is less obvious, however, how to extend the eigenvalue ratio (ER) test of AhnHorenstein2013 to nonlinear models. This method is attractive because it does not depend on somewhat arbitrary functional form assumptions or tuning parameters. It only requires to specify $R_{\max}$, but there is no penalty function or any other tuning parameter. Assuming that there exists an upper bound $R_{\max} > R$, we propose adapting this method to nonlinear factor single-index models using the following algorithm:
This algorithm can be seen as a natural generalization of the AhnHorenstein2013 to single-index models. Indeed, if we applied it to the linear model $Y_{ij} = X_{ij}' \beta + \alpha_i' \, \gamma_j + \varepsilon_{ij}$, with $\log f(Y_{ij} \mid X_{ij}' \beta + \alpha_i' \gamma_j )$ replaced by $-(Y_{ij} - X_{ij}' \beta - \alpha_i' \, \gamma_j)^2$, then $$\lambda_r\left( \widetilde \pi \widetilde \pi' \right) = \lambda_r\left[ \left(Y_{ij} - X_{ij}' \widetilde \beta \right) \left(Y_{ij} - X_{ij}' \widetilde \beta \right) ' \right],$$ which corresponds to the eigenvalue ratio criterion of AhnHorenstein2013 applied to the residuals $Y_{ij} - X_{ij}' \widetilde \beta $. Based on this coverage of the linear model, we conjecture that $\widehat R$ is a consistent estimator of $R$ under suitable conditions. To formalize this argument, a key step is to establish the consistency of the preliminary estimator $\widetilde \beta $, extending the results of MoonWeidner2015 from linear to nonlinear models, and the properties of the estimator of the factor structure $\widetilde \pi$. The main technical challenge is to characterize $\widetilde \pi$, which is not even available for the linear model with covariates and $R>R_0$. We leave this analysis to future research. In the rest of the section we show that the method performs well in numerical simulations.
To show how $\widehat R$ performs in small samples, we generate samples from the Poisson model of Example (ref) with additive effects where $z_{ij} = X_{ij} \beta + \alpha_{1i} + \gamma_{1j} + \alpha_{2i}' \gamma_{2j}$, $X_{ij} \sim N(1,1/3)$, $\beta = 0$, $\alpha_{1i} \sim U(0,1)$, $\gamma_{1i} \sim U(0,1)$, $\alpha_{2i}$ is an $R_2$-dimensional standard normal vector with independent components, $\gamma_{2i}$ is an $R_2$-dimensional standard normal vector with independent components, and $X_{ij}$, $\alpha_{1i'}$, $\gamma_{1j'}$, $\alpha_{2i''}$ and $\gamma_{2j''}$ are mutually independent for all $i,i',i'' = 1,\ldots,I$ and $j,j',j'' = 1,\ldots,J$. We generate $1,000$ datasets with $I=J\in\{50,75, 100,150\}$ and $R_2 \in \{1,2,3\}$, and apply Algorithm (ref) with $R_{\max} \in \{4, 5,6 \}$. Table (ref) reports the average of $\widehat R_2$ across simulations and the proportion of simulations where $\widehat R_2 = R_2$. Here, we find that $\widehat R_2$ has little bias and often yields the true $R_2$, specially for the larger sample sizes with $I \geq 75$. Interestingly, the performance of $\widehat R_2$ improves as $R_{\max}$ gets closer to $R_2$. Given this sensitivity, we recommend computing $\widehat R_2$ for several values of $R_{\max}$.
The gravity equation is a fundamental empirical relationship in international economics. We estimate a gravity equation of trade between countries using data from Helpman01052008 on bilateral trade flows and other trade-related variables for 157 countries in 1986.\footnote{The original data set includes 158 countries. We exclude Congo because it did not export to any other country in 1986.} The data set contains a network of trade data where both $i$ and $j$ index countries as senders (exporters) and receivers (importers), such that $I = J = 157$. The outcome $Y_{ij}$ is the volume of trade in thousands of constant 2000 US dollars from country $i$ to country $j$, and the covariates $X_{ij}$ include determinants of bilateral trade flows such as the logarithm of the distance in kilometers between country $i$'s capital and country $j$'s capital and indicators for common colonial ties, currency union, regional free trade area (FTA), border, legal system, language, and religion. Table (ref) reports descriptive statistics of the variables used in the analysis. There are $157 \times 156 = 24,492$ observations corresponding to different pairs of countries. The observations with $i = j$ are missing because we do not observe trade flows from a country to itself. The trade variable in the first row is an indicator of positive volume of trade. There are no trade flows for 55% of the country pairs.
We estimate a Poisson model with the following specification of the intensity $$ \mathbb{E}[Y_{ij} \mid X_{ij}, \alpha_{1i}, \gamma_{1j}, \alpha_{2i}, \gamma_{2j} ] = \exp(X_{ij}'\beta + \alpha_{1i} + \gamma_{1j} + \alpha_{2i}'\gamma_{2j}), $$ where $\alpha_{2i}$ and $\gamma_{2i}$ are $R_2$-dimensional vectors of factors and factor loadings. This model is a special case of Example (ref) with $\alpha_i = (\alpha_{1i}, 1, \alpha_{2i}')'$, $\gamma_j = (1,\gamma_{1j}, \gamma_{2j}')'$, and $R = 2 + R_2$. We explicitly include additive importer and exporter effects to account for scale and multilateral resistance effects following EatonKortum2001 and AndersonWincoop2003. Moreover, we also include interactive country effects to capture possible clustering and homophily induced by latent factors such as country trade partnerships, presence of multinationals or immigrant communities, or differences in natural resources or industrial composition.
Table (ref) reports the estimates and standard errors of the parameter $\beta$.\footnote{We do not report estimates of APEs because in the specification of the Poisson model that we use the parameters can be interpreted as elasticities. } We consider specifications with different number of interactive effects, $R_2$, in addition to the additive effects . The last row of the table reports the maximum value of the average log-likelihood, $L(\widehat \beta, \widehat \phi_n)/n $. We report two sets of standard errors corresponding to the dependence structures of cases (a) and (b) of Assumption (ref)(i). The standard errors in brackets account for possible reciprocity in the data. In this case, the method of Section (ref) selects $R_2 = 3$ factors when $R_{\max} = 4$ and $R_{\max} = 5$. We take $R_2 = 3$ as our preferred specification, but we also note that, relative to the standard errors, the estimates are not very sensitive to the $R_2$ in the range of values that we consider. One possible concern with the use of the Poisson model in the trade application is the excess zeros, i.e. the high probability of zero trade.\footnote{We thank an anonymous referee for raising this issue.} In this case, however, it does not seem to be a problem because the estimated model with $R_2=3$ predicts a probability of zero trade of $0.61$, which is higher than the observed probability of $0.55$.
We find that the sign of most of the effects is robust to the inclusion of latent factors. The only exceptions are the effects of common religion and language, which in the specification with only additive effects have counterintuitive negative signs that turn positive in our preferred specification. Comparing across columns, we observe that the model without factors seems to exaggerate the role of common border, whereas it downplays the effect of distance and colonial links. For example, increasing by 10% the distance reduces by 6.9% the volume of trade and sharing border increases it by 36% according to our preferred specification with $R_2 = 3$, whereas the same effects are 6% and 71% according to the specification with $R_2=0$. Except for language, all the coefficients are individually significant at the 5% level. Overall, increasing the number of factors makes the estimates less precise due to the loss of degrees of freedom. This observation showcases a trade-off in estimation between efficiency and robustness to richer dependence structures in the unobservables. Finally, accounting for reciprocity slightly increases the standard errors, but does not change the statistical significance of the estimates.
We evaluate the finite-sample properties of our estimation and inference methods in a Monte Carlo simulation that mimics the trade application. The design is calibrated to the Poisson model with additive importer and exporter country effects and one factor. We analyze the performance of the estimator of $\beta$ in terms of bias, dispersion and inference accuracy. To speed up computation, we include only one covariate: the log distance. More specifically, we generate $Y_{ij}$ from a Poisson distribution with intensity $\exp(X_{ij} \widehat{\beta} +\widehat{\alpha}_{1i}+\widehat{\gamma}_{1j}+\widehat{\alpha}_{2i}\widehat{\gamma}_{2j})$ independently across $i$ and $j$, where $X_{ij}$ takes the values of log-distance in the trade data set, and $\widehat{\beta}$ and $\{\widehat{\alpha}_{1i}, \widehat{\alpha}_{2i},\widehat{\gamma}_{1i}, \widehat{\gamma}_{2i}\}_{j=1}^{157}$, are equal to the estimates of the parameter, importer effects, exporter effects, factors and factor loadings. We repeat this procedure in $1,000$ simulations for four different sample sizes: $I=50$, $I=75$, $I=100$ and $I=157$ (full sample in the application). For each sample size and simulation, we draw a random sample of $I$ countries both as importers and exporters without replacement, so that the number of observations is $I\times(I-1)$. For each simulated sample, we reestimate the model parameter and standard errors, and construct 95% confidence interval for the model parameter.
Table (ref) reports the bias (Bias), standard deviation (SD), and root mean squared error (RMSE) of the estimator of the parameter $\beta$, together with the ratio of average standard error to the simulation standard deviation (SE/SD), and the empirical coverage in percentage of a confidence interval with 95% nominal value (p;95). We estimate models with four different numbers of factors in addition to the additive effects, $R_2 \in \{1,2,3, R_2^*\}$, where $R_2^*$ is the number of factors selected by the method of Section (ref) with $R_{\max} = 4$, which can vary across simulations. The results for the bias, SD and RMSE are reported in percentage of the true parameter value. We find that the bias is smaller than the standard deviation for every sample size. When we use the true number of factors $R_2=1$, the confidence intervals cover the parameter in more than 95% of the simulations. The excess coverage is due to the overestimation of the dispersion of the estimators by the standard errors. Selecting the number of factors does not introduce bias, but increases the dispersion of the estimator of the parameter. The additional variability yields slight undercoverage of the confidence intervals for small sample sizes. On the other hand, adding unnecessary factors to the specification increases the bias and dispersion of the estimator, but the confidence intervals continue having good coverage properties. This robustness to the inclusion of too many factors is consistent with the theoretical results of MoonWeidner2015 for linear factor models. Overall, the simulations show that the asymptotic theory of Section (ref) provides a good approximation to the finite-sample behavior of the estimator.