Individual and Time Effects in Nonlinear Panel Models with Large N, T

Individual and Time Effects in Nonlinear Panel Models with Large $N$, $T$

\abstract{

}

{

{ Keywords:} Panel data, nonlinear model, dynamic model, asymptotic bias correction, fixed effects, time effects.\\ { JEL:} C13, C23. }

Introduction

Fixed effects estimators of nonlinear panel data models can be severely biased because of the incidental parameter problem Neyman:1948p881. A growing literature, surveyed in Arellano and Hahn ArellanoHahn2007, shows that the leading term of an asymptotic expansion of the bias as both the cross-sectional dimension $N$ and time series dimension $T$ of the panel grow, can be characterized and corrected for. In models with individual effects, the leading bias term is of order $1/T$ and comes from the estimation of the individual effects. This result, however, does not apply to models with individual and time effects, where both of these effects are treated as parameters to be estimated. In this paper we show that the estimation of the time effects causes an additional incidental parameter bias of order $1/N$. Thus, if $N$ and $T$ are similarly large, the bias produced by the estimation of the time effects is of similar order of magnitude to the bias produced by the estimation of the individual effects, and both biases need to be corrected. We provide the corresponding analytical and jackknife bias corrections.

The asymptotic approximation to the fixed effects estimators that lets the two dimensions of the panel grow with the sample size is motivated by the recent availability of long panels and other large pseudo-panel data structures where the indexes might not correspond to individuals and time periods. Examples of these datasets include traditional microeconomic panel surveys with a long history of data such as the PSID and NLSY, international cross-country panels such as the Penn World Table, U.S. state level panels over time such as the CPS, and square pseudo-panels of trade flows across countries such as the Feenstra's World Trade Flows and CEPII, where the indexes correspond to the same countries indexed as importers and exporters.

We focus on semi-parametric models with log-likelihood functions that are concave in all parameters, and where each individual effect $\alpha_i$ and time effect $\gamma_t$ enter the log-likelihood for observation $(i,t)$ additively as $\alpha_i + \gamma_t$. This is the most common specification for the individual and time effects in linear models and is also a natural specification in the nonlinear models that we consider. Imposing concavity of the log-likelihood function greatly facilitates showing consistency in our setting where the dimension of the parameter space grows with the sample size. The most popular limited dependent variable models, including logit, probit, ordered probit, Tobit and Poisson models have concave log-likelihood functions, possibly after reparametrization (Olsen Olsen:1978p3375, and Pratt Pratt:1981p654). We note here that the general expansion that we derive in Appendix B do not impose additivity and concavity, but we use these restrictions to apply the expansion to fixed effects estimators. The models that we consider are semi-parametric because the joint distribution of the explanatory variables and the unobserved effects is left unspecified. The explanatory variables can be either strictly exogenous or predetermined.

We derive bias expansions and corrections for fixed effects estimators of common parameters $\beta$ and average partial effects (APEs). The vector $\beta$ includes all the unknown parameters that enter the log-likelihood function other than the individual and time effects, such as index coefficients in a probit model. The APEs are functions of the data, the common parameters, and the individual and time effects in nonlinear models. We find that the properties of the fixed effects estimators of $\beta$ and the APEs are different. For $\beta$, the order of the bias is $1/T + 1/N$, which is of the same as the rate of convergence $1/\sqrt{NT}$ under sequences where $N/T$ converge to a constant. For the APEs, we uncover that the incidental parameter problem is negligible asymptotically because the order of the bias, $1/N + 1/T$, is smaller than the rate of convergence, which is $1/\sqrt{N} + 1/\sqrt{T}$, slower than for model parameters. To the best of our knowledge, this rate result is new for fixed effects estimators of average partial effects in nonlinear panel models with individual and time effects.\footnote{Galvao and Kato GalvaoKato:2013 also found slow rates of convergence for fixed effects estimators in linear models with individual effects under misspecification. Fernandez-Val and Lee FL13 pointed out this issue in nonlinear models with only individual effects.} In numerical examples we find that the bias corrections, while not necessary to center the asymptotic distribution of APE estimators, do improve their finite-sample properties, specially in dynamic models.

The bias correction eliminates the bias terms of orders $1/T$ and $1/N$ from the fixed effects estimators. We considerer two methods to implement the correction: an analytical bias correction similar to Hahn and Newey Hahn:2004p882 and Hahn and Kuersteiner HahnKuersteiner2011, and a suitable modification of the split panel jackknife of Dhaene and Jochmans DhaeneJochmans2015.\footnote{A similar split panel jackknife bias correction method was outlined in Hu Hu2002.} However, the theory of the previous papers does not cover the models that we consider, because, in addition to not allowing for time effects, it assumes either identical distribution or stationarity over time for the processes of the observed variables, conditional on the unobserved effects. These assumptions are violated in our models due to the presence of the time effects, so we need to adjust the asymptotic theory accordingly. The individual and time effects introduce strong correlation in both dimensions of the panel. Conditional on the unobserved effects, we impose cross-sectional independence and weak time-serial dependence, and we allow for heterogeneity in both dimensions.

Simulation evidence indicates that our corrections improve the estimation and inference performance of the fixed effects estimators of parameters and average effects. The analytical corrections dominate the jackknife corrections in a probit model for sample sizes that are relevant for empirical practice. In the online supplement, Fern{\'a}ndez-Val and Weidner Supp2015, we illustrate the corrections with an empirical application on the relationship between competition and innovation using a panel of U.K. industries, following Aghion, Bloom, Blundell, Griffith and Howitt AghionBloomBlundellGriffithHowitt2005. We find that the inverted-U pattern relationship found by Aghion et al is robust to relaxing the strict exogeneity assumption of competition with respect to the innovation process and to the inclusion of innovation dynamics. We also uncover substantial state dependence in the innovation process.

\paragraph{Literature review.} The Neyman and Scott incidental parameter problem has been extensively discussed in the econometric literature; see, for example, Heckman Heckman:1981p2940, Lancaster Lancaster:2000p879, and Greene Greene:2004p3125. There is also a vast literature that shows how to tackle the problem in specific models under asymptotic sequences where $T$ is fixed and $N$ grows to infinity. However, there are results, e.g. from Honor{\'e} and Tamer HonoreTamer2006, Chamberlain Chamberlain2010, and Chernozhukov, Fern{\'a}ndez-Val, Hahn and Newey CFHN13, showing that model parameters and APEs are not point identified in important nonlinear panel data models under fixed-$T$ asymptotic sequences, implying that no fixed-$T$ consistent point estimators exist in these models.

A recent response to the incidental parameter problem is to adopt an alternative asymptotic approximation where both $N$ and $T$ grow with the sample size. Under these large-$T$ sequences, the fixed effects estimator is consistent but has bias in the asymptotic distribution. This asymptotic bias is the large-$T$ version of the incidental parameter problem and has motivated the development of bias corrections. Examples of papers that use this approximation include Phillips and Moon Phillips:1999p733, Hahn and Kuersteiner Hahn:2002p717, Lancaster Lancaster:2002p875, Woutersen Woutersen:2002p3683, Alvarez and Arellano AlvarezArellano2003, Hahn and Newey Hahn:2004p882, Carro Carro:2007p3601, Arellano and Bonhomme ArellanoBonhomme2009, Fernandez-Val FernandezVal:2009p3313, Hahn and Kuersteiner HahnKuersteiner2011, Fernandez-Val and Vella FernandezValVella2011, and Kato, Galvao and Montes-Rojas KatoGalvaoMontes-Rojas2012. This previous work, however, does not cover models with time effects.\footnote{An exception is Woutersen Woutersen:2002p3683, which considers a special type of grouped time effects whose number is fixed with $T$. We instead consider an unrestricted set of $T$ time effects, one for each time period.} Our contribution to this literature is to extend the large-$T$ bias corrections to models with two-way unobserved effects such as the individual and time effects commonly included in linear models.

The large-$T$ panel literature on models with both individual and time effects is sparse. Pesaran Pesaran2006, Bai Bai:2009p3321, and Moon and Weidner MoonWeidner2015a,MoonWeidner2015b study linear regression models with interactive individual and time fixed effects. The fixed effects estimators in these models also have asymptotic bias of order $1/T + 1/N$, but the methods used to derive this bias rely on linearity and therefore cannot be applied to the nonlinear models that we consider. Hahn and Moon HahnMoon2006 consider bias corrected fixed effects estimators in panel linear autoregressive models with additive individual and time effects. Regarding non-linear models, there is independent and contemporaneous work by Charbonneau Charbonneau2011,Charbonneau2014, which extends the conditional fixed effects estimators to logit and Poisson models with individual and time effects. She differences out the individual and time effects by conditioning on sufficient statistics. The conditional approach completely eliminates the asymptotic bias coming from the estimation of the incidental parameters, but it does not permit estimation of average partial effects and has not been developed for models with predetermined regressors. We instead consider estimators of model parameters and average partial effects in nonlinear models with predetermined regressors. The two approaches can therefore be considered as complementary.

\paragraph{Outline of the paper.} The rest of the paper is organized as follows. Section (ref) introduces the model and fixed effects estimators. Section (ref) describes the bias corrections to deal with the incidental parameters problem and illustrates how the bias corrections work through an example. Section (ref) provides the asymptotic theory. Section (ref) presents Monte Carlo results. The Appendix collects the proofs of the main results, and an online supplement to the paper contains additional technical derivations, numerical examples, and an empirical application Supp2015.

Model and Estimators

Model

The data consist of $N \times T$ observations $\{(Y_{it}, X'_{it})': 1 \leq i \leq N, 1 \leq t \leq T \},$ for a scalar outcome variable of interest $Y_{it}$ and a vector of explanatory variables $X_{it}$. We assume that the outcome for individual $i$ at time $t$ is generated by the sequential process:

where $X^t_i = (X_{i1}, \ldots, X_{it}),$ $\alpha = (\alpha_1, \ldots, \alpha_N)$, $\gamma = (\gamma_1, \ldots, \gamma_T)$, $f_{Y}$ is a known probability function, and $\beta$ is a finite dimensional parameter vector. The variables $\alpha_i$ and $\gamma_t$ are unobserved individual and time effects that in economic applications capture individual heterogeneity and aggregate shocks, respectively. The model is semiparametric because we do not specify the distribution of these effects nor their relationship with the explanatory variables. The conditional distribution $f_{Y}$ represents the parametric part of the model. The vector $X_{it}$ contains predetermined variables with respect to $Y_{it}$. Note that $X_{it}$ can include lags of $Y_{it}$ to accommodate dynamic models.

We consider two running examples throughout the analysis:

For estimation, we adopt a fixed effects approach, treating the realization of the unobserved individual and time effects as parameters to be estimated. We collect all these effects in the vector $\phi_{NT} = (\alpha_1, ..., \alpha_N, \gamma_1, ..., \gamma_T)'$. The model parameter $\beta$ usually includes regression coefficients of interest, while the vector $\phi_{NT}$ is treated as a nuisance parameter. The true values of the parameters, denoted by $\beta^0$ and $\phi_{NT}^0 = ({\alpha^{0}_1}, ..., {\alpha^{0}_N}, {\gamma^{0}_1}, ..., {\gamma^{0}_T})'$, are the solution to the population conditional maximum likelihood problem

for every $N,T$, where $\mathbb{E}_{\phi}$ denotes the expectation with respect to the distribution of the data conditional on the unobserved effects and initial conditions including strictly exogenous variables, $b > 0$ is an arbitrary constant, $v_{NT} = ( 1_N' ,- 1_T' )'$, and $1_N$ and $1_T$ denote vectors of ones with dimensions $N$ and $T$. Existence and uniqueness of the solution to the population problem will be guaranteed by our assumptions in Section (ref) below, including concavity of the objective function in all parameters. The second term of $\mathcal{L}_{NT}$ is a penalty that imposes a normalization needed to identify $\phi_{NT}$ in models with scalar individual and time effects that enter additively into the log-likelihood function as $\alpha_i + \gamma_t$.\footnote{ In Appendix (ref) we derive asymptotic expansions that apply to general models with multiple unobseved effects. In order to use these expansions to obtain the asymptotic distribution of the panel fixed effects estimators, we need to derive the properties of the expected Hessian of the incidental parameters, a matrix with increasing dimension, and to show the consistency of the estimator of the incidental parameter vector. The additive specification $\alpha_i + \gamma_t$ is useful to characterize the Hessian and we impose strict concavity of the objective function to show the consistency. } In this case, adding a constant to all $\alpha_i$, while subtracting it from all $\gamma_t$, does not change $\alpha_{i} + \gamma_t$. To eliminate this ambiguity, we normalize $\phi^0_{NT}$ to satisfy $v_{NT}' \phi^0_{NT} = 0$, i.e. $\sum_i \alpha_i^0 = \sum_t \gamma_t^0$. The penalty produces a maximizer of ${\cal L}_{NT}$ that is automatically normalized. We could equivalently impose $v_{NT}' \phi_{NT} = 0$ as a constraint, but for technical reasons we prefer to work with an unconstrained optimization problem. There are other possible normalizations for $\phi_{NT}$, such as $\alpha_1 = 0$. The model parameter $\beta$ is invariant to the choice of normalization, that is, our asymptotic results on the estimator for $\beta$ are independent of this choice of normalization. Our choice is convenient for certain intermediate results that involve the incidental parameter $\phi_{NT}$, its score vector and its Hessian matrix. The pre-factor $(NT)^{-1/2}$ in $\mathcal{L}_{NT}(\beta,\phi_{NT})$ is just a rescaling.

Other quantities of interest involve averages over the data and unobserved effects

where $\mathbb{E}$ denotes the expectation with respect to the joint distribution of the data and the unobserved effects, provided that the expectation exists. $\delta^0_{NT} $ is indexed by $N$ and $T$ because the marginal distribution of $\{(X_{it}, \alpha_i, \gamma_t) : 1 \leq i \leq N, 1 \leq t \leq T\}$ can be heterogeneous across $i$ and/or $t$; see Section (ref). These averages include average partial effects (APEs), which are often the ultimate quantities of interest in nonlinear models. The APEs are invariant to the choice of normalization for $\phi_{NT}$ if $\alpha_i$ and $\gamma_t$ enter $\Delta(X_{it}, \beta, \alpha_i, \gamma_t)$ as $\alpha_i + \gamma_t$. Some examples of partial effects that satisfy this condition are the following:

Example (ref) (Binary response model). If $X_{it,k}$, the $k$th element of $X_{it}$, is binary, its partial effect on the conditional probability of $Y_{it}$ is

where $\beta_k$ is the $k$th element of $\beta$, and $X_{it,-k}$ and $\beta_{-k}$ include all elements of $X_{it}$ and $\beta$ except for the $k$th element. If $X_{it,k}$ is continuous and $F$ is differentiable, the partial effect of $X_{it,k}$ on the conditional probability of $Y_{it}$ is

where $\partial F$ is the derivative of $F$. }

Example (ref) (Poisson model). If $X_{it}$ includes $Z_{it}$ and some known transformation $H(Z_{it})$ with coefficients $\beta_k$ and $\beta_j$, the partial effect of $Z_{it}$ on the conditional expectation of $Y_{it}$ is

}

Fixed effects estimators

We estimate the parameters by solving the sample analog of problem (ref), i.e.

As in the population case, we shall impose conditions guaranteeing that the solution to this maximization problem exists and is unique with probability approaching one as $N$ and $T$ become large. For computational purposes, we note that the solution to the program (ref) for $\beta$ is the same as the solution to the program that imposes $v_{NT}' \phi_{NT} = 0$ directly as a constraint in the optimization, and is invariant to the normalization. In our numerical examples we impose either $\alpha_1 = 0$ or $\gamma_1 = 0$ directly by dropping the first individual or time effect. This constrained program has good computational properties because its objective function is concave and smooth in all the parameters. We have developed the commands probitfe and logitfe in Stata to implement the methods of the paper for probit and logit models Stata2015.\footnote{We refer to this companion work for computational details.} When $N$ and $T$ are large, e.g., $N > 2,000$ and $T > 50$, we recommend the use of optimization routines that exploit the sparsity of the design matrix of the model to speed up computation such as the package Speedglm in R Speedglm2012. For a probit model with $N=2,000$ and $T = 52$, Speedglm computes the fixed effects estimator in less than 2 minutes with a 2 x 2.66 GHz 6-Core Intel Xeon processor, more than 7.5 times faster than our \texttt{Stata} command \texttt{probitfe} and more than 30 times faster than the \texttt{R} command \texttt{glm}.\footnote{Additional comparisons of computational times are available from the authors upon request.}

To analyze the statistical properties of the estimator of $\beta$ it is convenient to first concentrate out the nuisance parameter $\phi_{NT}$. For given $\beta$, we define the optimal $\widehat \phi_{NT}(\beta)$ as

The fixed effects estimators of $\beta^0$ and $\phi_{NT}^0$ are

Estimators of APEs can be formed by plugging-in the estimators of the model parameters in the sample version of (ref), i.e.

Again, $\widehat \delta_{NT}$ is invariant to the normalization chosen for $\phi_{NT}$ if $\alpha_i$ and $\gamma_t$ enter $\Delta(X_{it}, \beta, \alpha_i, \gamma_t)$ as $\alpha_i + \gamma_t$.

Incidental parameter problem and bias corrections

In this section we give a heuristic discussion of the main results, leaving the technical details to Section (ref). We illustrate the analysis with numerical calculations based on a variation of the classical Neyman and Scott Neyman:1948p881 variance example.

Incidental parameter problem

Fixed effects estimators in nonlinear models suffer from the incidental parameter problem Neyman:1948p881. The source of the problem is that the dimension of the nuisance parameter $\phi_{NT}$ increases with the sample size under asymptotic approximations where either $N$ or $T$ pass to infinity. To describe the problem let

The fixed effects estimator is inconsistent under the traditional Neyman and Scott asymptotic sequences where $N \to \infty$ and $T$ is fixed, i.e., $\operatorname*{plim}_{N \to \infty} \overline{\beta}_{NT} \neq \beta^0$. Similarly, the fixed effects estimator is inconsistent under asymptotic sequences where $T \to \infty$ and $N$ is fixed, i.e., $\operatorname*{plim}_{T \to \infty} \overline{\beta}_{NT} \neq \beta^0$. Note that $\overline \beta_{NT} = \beta^0$ if $\widehat \phi_{NT}(\beta)$ is replaced by $\phi_{NT}(\beta) = \operatorname*{argmax}_{\phi_{NT} \in \mathbb{R}^{\dim \phi_{NT}}} \, \mathbb{E}_{\phi}[ {{\cal L}}_{NT}(\beta,\,\phi_{NT})]$. Under asymptotic approximations where either $N$ or $T$ are fixed, there is only a fixed number of observations to estimate some of the components of $\phi_{NT}$, $T$ for each individual effect or $N$ for each time effect, rendering the estimator $\widehat \phi_{NT}(\beta)$ inconsistent for $\phi_{NT}(\beta)$. The nonlinearity of the model propagates the inconsistency to the estimator of $\beta$.

A key insight of the large-$T$ panel data literature is that the incidental parameter problem becomes an asymptotic bias problem under an asymptotic approximation where $N \to \infty$ and $T \to \infty$ (e.g., Arellano and Hahn, ArellanoHahn2007). For models with only individual effects, this literature derived the expansion $\overline{\beta}_{NT} = \beta^0 + B/T + o_P(T^{-1})$ as $N,T \to \infty$, for some constant $B$. The fixed effects estimator is consistent because $\operatorname*{plim}_{N,T \to \infty} \overline{\beta}_{NT} = \beta^0$, but has bias in the asymptotic distribution if $B/T$ is not negligible relative to $1/\sqrt{NT}$, the order of the standard deviation of the estimator. This asymptotic bias problem, however, is easier to correct than the inconsistency problem that arises under the traditional Neyman and Scott asymptotic approximation. We show that the same insight still applies to models with individual and time effects, but with a different expansion for $\overline{\beta}_{NT}$. We characterize the expansion and develop bias corrections.

Bias Expansions and Bias Corrections

Some expansions can be used to explain our corrections. For smooth likelihoods and under appropriate regularity conditions, as $N,T \to \infty$,

for some $\overline B_{\infty}^{\beta}$ and $\overline D_{\infty}^{\beta}$ that we characterize in Theorem (ref) and explain in Remark (ref), where $a \vee b := \max(a,b)$. Unlike in nonlinear models without incidental parameters, the order of the bias is higher than the inverse of the sample size $(NT)^{-1}$ due to the slow rate of convergence of $\widehat \phi_{NT}$. Note also that by the properties of the maximum likelihood estimator $$ \sqrt{NT}(\widehat \beta_{NT} - \overline{\beta}_{NT}) \to_d \mathcal{N}(0, \overline V_{\infty}), $$ for some $\overline V_{\infty}$ that we also characterize in Theorem (ref). Under asymptotic sequences where $N/T \to \kappa^2$ as $N,T \to \infty$, the fixed effects estimator is asymptotically biased because

Relative to fixed effects estimators with only individual effects, the presence of time effects introduces additional asymptotic bias through $\overline D_{\infty}^{\beta}$. This asymptotic result predicts that the fixed effects estimator can have significant bias relative to its dispersion. Moreover, confidence intervals constructed around the fixed effects estimator can severely undercover the true value of the parameter even in large samples. We show that these predictions provide a good approximations to the finite sample behavior of the fixed effects estimator through analytical and simulation examples in Sections (ref) and (ref).

The analytical bias correction consists of subtracting estimates of the leading terms of the bias from the fixed effect estimator of $\beta^0$. Let $\widehat B_{NT}^{\beta}$ and $\widehat D_{NT}^{\beta}$ be estimators of $\overline B_{\infty}^{\beta}$ and $\overline D_{\infty}^{\beta}$ as defined in (ref). The bias corrected estimator can be formed as

If $N/T \to \kappa^2$, $\widehat B_{NT}^{\beta} \to_P \overline B_{\infty}^{\beta}$, and $\widehat D_{NT}^{\beta} \to_P \overline D_{\infty}^{\beta},$ then

The analytical correction therefore centers the asymptotic distribution at the true value of the parameter, without increasing asymptotic variance. This asymptotic result predicts that in large samples the corrected estimator has small bias relative to dispersion, the correction does not increase dispersion, and the confidence intervals constructed around the corrected estimator have coverage probabilities close to the nominal levels. We show that these predictions provide a good approximations to the behavior of the corrections in Sections (ref) and (ref) even in small panels with $N < 60$ and $T < 15$.

We also consider a jackknife bias correction method that does not require explicit estimation of the bias. This method is based on the split panel jackknife (SPJ) of Dhaene and Jochmans DhaeneJochmans2015 applied to the time and cross-section dimension of the panel. Alternative jackknife corrections based on the leave-one-observation-out panel jackknife (PJ) of Hahn and Newey Hahn:2004p882 and combinations of PJ and SPJ are also possible. We do not consider corrections based on PJ because they are theoretically justified by second-order expansions of $\overline{\beta}_{NT}$ that are beyond the scope of this paper.

To describe our generalization of the SPJ, define the fixed effects estimator of $\beta$ in the subpanel with cross sectional indexes $A$ and time series indexes $B$ as $$ \widehat \beta_{A,B} \in \operatorname*{argmax}_{\beta \in \mathbb{R}^{\dim \beta} } \; \max_{\alpha(A) \in \mathbb{R}^{|A|}} \; \max_{\gamma(B) \in \mathbb{R}^{|B|}} \; \sum_{i,t} d_{it}(A,B) \, \log f_{Y}(Y_{it} \mid X_{it}, \alpha_i, \gamma_t, \beta) , $$ where $\alpha(A)=( \alpha_i : i \in A )$, $\gamma(B)=( \gamma_t : t \in B )$, and $d_{it}(A,B) = 1(i \in A) \times 1(t \in B)$. Let $\widetilde{\beta}_{N,T/2}$ be the average of the 2 split jackknife estimators in the subpanels with $A = \{1,2,\ldots,N\}$, and $B = \{1,2,\ldots,T/2\}$ or $B = \{T/2+1, T/2+2, \ldots, T\}$, i.e. including all the individuals and leaving out the first and second halves of the time periods. Let $\widetilde{\beta}_{N/2,T}$ be the average of the 2 split jackknife estimators in the subpanels with $B = \{1,2,\ldots,T\}$, and $A = \{1,2,\ldots,N/2\}$ or $A = \{N/2+1,N/2+2,\ldots,N\}$, i.e. including all the time periods and leaving out half of the individuals of the panel.\footnote{When $T$ is odd we define $\widetilde{\beta}_{N,T/2}$ as the average of the 2 split jackknife estimators that use overlapping subpanels with $B = \{1,2,\ldots,(T+1)/2\}$ and $B = \{(T+1)/2, (T+1)/2+1, \ldots,T\}$. We define $\widetilde{\beta}_{N/2,T}$ similarly when $N$ is odd.} In choosing the cross sectional indexing of the panel, one might want to take into account individual clustering structures and other dependencies to preserve them in the SPJ. For example, all the individuals belonging to the same cluster should be indexed such that they remain in the same subpanel after the cross sectional split. If there are no cross sectional dependencies, the indexing of the individuals is unrestricted. We recommend to construct $\widetilde{\beta}_{N/2,T}$ as the average of the estimators obtained from all possible partitions of $N/2$ individuals to avoid ambiguity and arbitrariness in the choice of the division.\footnote{There are $P = {N \choose N/2}$ different cross sectional partitions with $N/2$ individuals. When $N$ is large, we can approximate the average over all possible partitions by the average over $S \ll P$ randomly chosen partitions to speed up computation.} The bias corrected estimator is

To give some intuition about how the corrections works, note that

where $ \widetilde{\beta}_{N,T/2} - \widehat \beta_{NT} = \overline B_{\infty}^{\beta}/T + o_P(T^{-1} \vee N^{-1})$ and $\widetilde{\beta}_{N/2,T} - \widehat \beta_{NT} = \overline D_{\infty}^{\beta}/N + o_P(T^{-1} \vee N^{-1}).$ Relative to $\widehat \beta_{NT}$, $ \widetilde{\beta}_{N,T/2} $ has double the bias coming from the estimation of the individual effects because it is based on subpanels with half of the time periods, and $ \widetilde{\beta}_{N/2,T} $ has double the bias coming from the estimation of the time effects because it is based on subpanels with half of the individuals. The time series split removes the bias term $\overline B_{\infty}^{\beta}$ and the cross sectional split removes the bias term $\overline D_{\infty}^{\beta}.$

Illustrative Example

To illustrate how the bias corrections work in finite samples, we consider a simple model where the solution to the population program (ref) has closed form. This model corresponds to a variation of the classical Neyman and Scott Neyman:1948p881 variance example that includes both individual and time effects, $Y_{it} \mid \alpha, \gamma, \beta \sim \mathcal{N}(\alpha_i + \gamma_t, \beta)$. It is well-know that in this case $$ \widehat \beta_{NT} = (NT)^{-1} \sum_{i,t} \left( Y_{it} - \bar Y_{i.} - \bar Y_{.t} + \bar Y_{..} \right)^2, $$ where $\bar Y_{i.} = T^{-1} \sum_t Y_{it}$, $\bar Y_{.t} = N^{-1} \sum_i Y_{it},$ and $\bar Y_{..} = (NT)^{-1} \sum_{i,t} Y_{it}.$ Moreover, from the well-known results on the degrees of freedom adjustment of the estimated variance $$ \overline{\beta}_{NT} = \mathbb{E}_{\phi} [\widehat \beta_{NT}] = \beta^0 \frac{(N-1)(T-1)}{NT} = \beta^0 \left(1 - \frac{1}{ T} - \frac{1}{ N} + \frac{1}{NT} \right), $$ so that $\overline{B}_{\infty}^{\beta} = -\beta^0$ and $\overline{D}_{\infty}^{\beta} = -\beta^0$.\footnote{Okui Okui2013 derived the bias of fixed effects estimators of autocovariances and autocorrelations in this model.}

To form the analytical bias correction we can set $\widehat{B}_{NT}^{\beta} = - \widehat \beta_{NT}$ and $\widehat{D}_{NT}^{\beta} = - \widehat \beta_{NT}$. This yields $\widetilde \beta^A_{NT} = \widehat \beta_{NT} (1 + 1/T + 1/N)$ with $$ \overline{\beta}^A_{NT} = \mathbb{E}_{\phi}[\widetilde \beta^A_{NT}] = \beta^0 \left( 1 - \frac{1}{T^2} - \frac{1}{N^2} - \frac{1}{NT} + \frac{1}{NT^2} + \frac{1}{N^2T}\right). $$ This correction reduces the order of the bias from $(T^{-1} \vee N^{-1})$ to $(T^{-2} \vee N^{-2}),$ and introduces additional higher order terms. The analytical correction increases finite-sample variance because the factor $(1 + 1/T + 1/N) > 1$. We compare the biases and standard deviations of the fixed effects estimator and the corrected estimator in a numerical example below.

For the Jackknife correction, straightforward calculations give $$ \overline{\beta}^{J}_{NT} = \mathbb{E}_{\phi}[\widetilde \beta^{J}_{NT}] = 3 \overline{\beta}_{NT} - \overline{\beta}_{N,T/2} - \overline{\beta}_{N/2,T} = \beta^0 \left( 1 - \frac{1}{NT} \right).$$ The correction therefore reduces the order of the bias from $(T^{-1} \vee N^{-1})$ to $(TN)^{-1}.$\footnote{ In this example it is possible to develop higher-order jackknife corrections that completely eliminate the bias because we know the entire expansion of $\overline{\beta}_{NT}$. For example, $\mathbb{E}_{\phi}[ 4 \widehat \beta_{NT} - 2 \widetilde \beta_{N,T/2} - 2 \widetilde \beta_{N/2,T} + \widetilde \beta_{N/2,T/2}] = \beta^0,$ where $\widetilde \beta_{N/2,T/2}$ is the average of the four split jackknife estimators that leave out half of the individuals and the first or the second halves of the time periods. See Dhaene and Jochmans DhaeneJochmans2015 for a discussion on higher-order bias corrections of panel fixed effects estimators.}

Table (ref) presents numerical results for the bias and standard deviations of the fixed effects and bias corrected estimators in finite samples. We consider panels with $N,T \in \{10,25, 50\},$ and only report the results for $T \leq N$ since all the expressions are symmetric in $N$ and $T$. All the numbers in the table are in percentage of the true parameter value, so we do not need to specify the value of $\beta^0$. We find that the analytical and jackknife corrections offer substantial improvements over the fixed effects estimator in terms of bias. The first and fourth row of the table show that the bias of the fixed effects estimator is of the same order of magnitude as the standard deviation, where $\overline{V}_{NT} = {\rm Var}[\widehat \beta_{NT}] = 2 (N-1)(T-1) (\beta^0)^2 / (NT)^2$ under independence of $Y_{it}$ over $i$ and $t$ conditional on the unobserved effects. The fifth row shows the increase in standard deviation due to analytical bias correction is small compared to the bias reduction, where $\overline{V}_{NT}^A = {\rm Var}[\widetilde \beta_{NT}^A] = (1+1/N+1/T)^2 \overline{V}_{NT}$. The last row shows that the jackknife yields less precise estimates than the analytical correction when $T=10$.

Table (ref) illustrates the effect of the bias on the inference based on the asymptotic distribution. It shows the coverage probabilities of 95% asymptotic confidence intervals for $\beta^0$ constructed in the usual way as $$\text{CI}_{.95}(\widehat \beta) = \widehat \beta \pm 1.96 \widehat V_{NT}^{1/2} = \widehat \beta (1 \pm 1.96 \sqrt{2/(NT)}),$$ where $\widehat \beta = \{\widehat \beta_{NT}, \widetilde \beta_{NT}^{A}, \widetilde \beta_{NT}^{J}\}$ and $\widehat V_{NT} = 2 \widehat \beta^2/(NT)$ is an estimator of the asymptotic variance $\overline{V}_{\infty}/(NT) = 2 (\beta^0)^2/(NT)$. To find the coverage probabilities, we use that $NT \widehat \beta_{NT}/ \beta^0 \sim \chi^2_{(N-1)(T-1)}$ and $\widetilde \beta_{NT}^{A} = (1 + 1/N + 1/T) \widehat \beta_{NT}$. These probabilities do not depend on the value of $\beta^0$ because the limits of the intervals are proportional to $\widehat \beta$. For the Jackknife we compute the probabilities numerically by simulation with $\beta^0=1$. As a benchmark of comparison, we also consider confidence intervals constructed from the unbiased estimator $\widetilde \beta_{NT} = NT \widehat \beta_{NT}/[(N-1)(T-1)]$. Here we find that the confidence intervals based on the fixed effect estimator display severe undercoverage for all the sample sizes. The confidence intervals based on the corrected estimators have high coverage probabilities, which approach the nominal level as the sample size grows. Moreover, the bias corrected estimators produce confidence intervals with very similar coverage probabilities to the ones from the unbiased estimator.

Asymptotic Theory for Bias Corrections

In nonlinear panel data models the population problem ((ref)) generally does not have closed form solution, so we need to rely on asymptotic arguments to characterize the terms in the expansion of the bias ((ref)) and to justify the validity of the corrections.

Asymptotic distribution of model parameters

We consider panel models with scalar individual and time effects that enter the likelihood function additively through $\pi_{it} = \alpha_i + \gamma_t$. In these models the dimension of the incidental parameters is $\dim \phi_{NT} = N+T$. The leading cases are single index models, where the dependence of the likelihood function on the parameters is through an index $ X_{it}'\beta + \alpha_i + \gamma_t$. These models cover the probit and Poisson specifications of Examples (ref) and (ref). The additive structure only applies to the unobserved effects, so we can allow for scale parameters to cover the Tobit and negative binomial models. We focus on these additive models for computational tractability and because we can establish the consistency of the fixed effects estimators under a concavity assumption in the log-likelihood function with respect to all the parameters.

The parametric part of our panel models takes the form

We denote the derivatives of the log-likelihood function $\ell_{it}$ by $\partial_\beta \ell_{it}(\beta, \pi) := \partial \ell_{it}(\beta,\pi)/\partial \beta$, $\partial_{\beta \beta'} \ell_{it}(\beta, \pi) := \partial^2 \ell_{it}(\beta,\pi)/(\partial \beta \partial \beta')$, $\partial_{\pi^q} \ell_{it}(\beta, \pi) := \partial^q \ell_{it}(\beta,\pi)/\partial \pi^q$, $q = 1,2,3$, etc. We drop the arguments $\beta$ and $\pi$ when the derivatives are evaluated at the true parameters $\beta^0$ and $\pi^0_{it} := \alpha_i^0 + \gamma_t^0$, e.g. $\partial_{\pi^q} \ell_{it} := \partial_{\pi^q}\ell_{it}(\beta^0,\pi^0_{it})$. We also drop the dependence on $NT$ from all the sequences of functions and parameters, e.g. we use $\mathcal{L}$ for $\mathcal{L}_{NT}$ and $\phi$ for $\phi_{NT}$.

We make the following assumptions:

To describe the asymptotic distribution of the fixed effects estimator $\widehat \beta,$ it is convenient to introduce some additional notation. Let $\overline{\cal H}$ be the $(N+T) \times (N+T)$ expected Hessian matrix of the log-likelihood with respect to the nuisance parameters evaluated at the true parameters, i.e.

where $\overline{\mathcal{H}}_{(\alpha\alpha)}^* = \text{diag}( \sum_{t} \mathbb{E}_{\phi}[- \partial_{\pi^2} \ell_{it}])/\sqrt{NT}$, $\overline{\mathcal{H}}_{(\alpha\gamma)it}^* = \mathbb{E}_{\phi}[- \partial_{\pi^2} \ell_{it}]/\sqrt{NT}$, and $\overline{\mathcal{H}}_{(\gamma\gamma)} ^*= \linebreak \text{diag}( \sum_{i} \mathbb{E}_{\phi}[- \partial_{\pi^2} \ell_{it}])/\sqrt{NT}$. Furthermore, let $\overline{\cal H}^{-1}_{(\alpha\alpha)}$, $\overline{\cal H}^{-1}_{(\alpha\gamma)}$, $\overline{\cal H}^{-1}_{(\gamma\alpha)}$, and $\overline{\cal H}^{-1}_{(\gamma\gamma)}$ denote the $N\times N$, $N\times T$, $T\times N$ and $T\times T$ blocks of the inverse $\overline{\cal H}^{-1}$ of $\overline{\cal H}$. We define the $\dim \beta$-vector $\Xi_{it}$ and the operator $D_{\beta \pi^q}$ as

with $q=0,1,2$. The $k$-th component of $\Xi_{it}$ corresponds to the population least squares projection of $\mathbb{E}_{\phi}(\partial_{\beta_k \pi} \ell_{it})/\mathbb{E}_{\phi}( \partial_{\pi^2} \ell_{it} )$ on the space spanned by the incidental parameters under a metric given by $ \mathbb{E}_{\phi}( - \partial_{\pi^2} \ell_{it})$, i.e.

The operator $D_{\beta \pi^q}$ partials out individual and time effects in nonlinear models. It corresponds to individual and time differencing when the model is linear. To see this, consider the normal linear model $Y_{it} \mid X_{i}^t, \alpha_i,\gamma_t \sim {\cal N}(X_{it}'\beta + \alpha_i + \gamma_t, 1).$ Then, $\Xi_{it} = T^{-1} \sum_{t=1}^T \mathbb{E}_{\phi}[X_{it}] + N^{-1} \sum_{i=1}^N \mathbb{E}_{\phi}[X_{it}] - (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T \mathbb{E}_{\phi}[X_{it}]$, $D_{\beta} \ell_{it} = - \tilde X_{it} \varepsilon_{it}, $ $D_{\beta \pi} \ell_{it} = - \tilde X_{it},$ and $D_{\beta \pi^2} \ell_{it} = 0$, where $\varepsilon_{it} = Y_{it} - X_{it}'\beta - \alpha_i - \gamma_t$ and $\tilde X_{it} = X_{it} - \Xi_{it}$ is the individual and time demeaned explanatory variables.

The following theorem establishes the asymptotic distribution of the fixed effects estimator $\widehat \beta.$

It is instructive to evaluate the expressions of the bias in our running examples.

Example (ref) (Binary response model). In this case $$\ell_{it}(\beta,\pi) = Y_{it} \log F(X_{it}'\beta + \pi) + (1 - Y_{it}) \log [1 - F(X_{it}'\beta + \pi)],$$ so that $\partial_{\pi} \ell_{it} = H_{it} (Y_{it} - F_{it}),$ $\partial_{\beta} \ell_{it} = \partial_{\pi} \ell_{it} X_{it},$ $ \partial_{{\pi}^2} \ell_{it} = - H_{it} \partial F_{it} + \partial H_{it} (Y_{it} - F_{it})$, $\partial_{\beta \beta'} \ell_{it} = \partial_{\pi^2} \ell_{it} X_{it} X_{it}'$, $\partial_{\beta \pi} \ell_{it} = \partial_{\pi^2} \ell_{it} X_{it},$ $\partial_{{\pi}^3} \ell_{it} = - H_{it} \partial^2 F_{it} - 2 \partial H_{it} \partial F_{it} + \partial^2 H_{it} (Y_{it} - F_{it})$, and $\partial_{\beta {\pi}^2} \ell_{it} = \partial_{{\pi}^3} \ell_{it} X_{it}$, where $H_{it} = \partial F_{it} / [F_{it}(1-F_{it})], and $ $\partial^{j} G_{it} := \partial^{j} G(Z)|_{Z = X_{it}'\beta^0 + \pi_{it}^0}$ for any function $G$ and $j = 0,1,2$. Substituting these values in the expressions of the bias of Theorem (ref) yields

where $\omega_{it} = H_{it} \partial F_{it}$ and $\tilde X_{it}$ is the residual of the population projection of $X_{it}$ on the space spanned by the incidental parameters under a metric weighted by $\mathbb{E}_{\phi}(\omega_{it})$. For the probit model where all the components of $X_{it}$ are strictly exogenous, $$ \overline B_{\infty} = \overline{\mathbb{E}} \left[ \frac {1} {2N} \sum_{i=1}^{N} \frac{ \sum_{t=1}^T \mathbb{E}_{\phi}[\omega_{it} \tilde{X}_{it} \tilde{X}_{it}'] } { \sum_{t=1}^T \mathbb{E}_{\phi}\left( \omega_{it} \right) }\right] \beta^0 , \ \ \overline D_{\infty} = \overline{\mathbb{E}} \left[ \frac {1} {2T} \sum_{t=1}^{T} \frac{ \sum_{i=1}^N \mathbb{E}_{\phi}[\omega_{it} \tilde{X}_{it} \tilde{X}_{it}'] } { \sum_{i=1}^N \mathbb{E}_{\phi}\left( \omega_{it} \right) } \right] \beta^0. $$ 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 only individual effects FernandezVal:2009p3313.}

Example (ref) (Poisson model). In this case $$\ell_{it}(\beta,\pi) = (X_{it}'\beta + \pi) Y_{it} - \exp(X_{it}'\beta + \pi) - \log Y_{it}!,$$ so that $\partial_{\pi} \ell_{it} = Y_{it} - \omega_{it},$ $\partial_{\beta} \ell_{it} = \partial_{\pi} \ell_{it} X_{it},$ $\partial_{{\pi}^2} \ell_{it} = \partial_{{\pi}^3} \ell_{it}= - \omega_{it} $, $\partial_{\beta \beta'} \ell_{it} = \partial_{{\pi}^2} \ell_{it} X_{it} X_{it}',$ and $\partial_{\beta {\pi}} \ell_{it} = \partial_{\beta {\pi}^2} \ell_{it} = \partial_{{\pi}^3} \ell_{it} X_{it},$ where $\omega_{it} = \exp(X_{it}'\beta^0 + \pi_{it}^0)$. Substituting these values in the expressions of the bias of Theorem (ref) yields

and $\overline D_{\infty} = 0$, where $\tilde X_{it}$ is the residual of the population projection of $X_{it}$ on the space spanned by the incidental parameters under a metric weighted by $\mathbb{E}_{\phi}(\omega_{it})$. If in addition all the components of $X_{it}$ are strictly exogenous, then we get the no asymptotic bias result $\overline B_{\infty} = \overline D_{\infty} = 0$.}

Asymptotic distribution of APEs

In nonlinear models we are often interested in APEs, in addition to model parameters. These effects are averages of the data, parameters and unobserved effects; see expression (ref). For the panel models of Assumption (ref) we specify the partial effects as $\Delta(X_{it}, \beta, \alpha_i, \gamma_t) = \Delta_{it}(\beta, \pi_{it})$. The restriction that the partial effects depend on $\alpha_i$ and $\gamma_t$ through $\pi_{it}$ is natural in our panel models since $$ \mathbb{E} [Y_{it} \mid X^t_i, \alpha_i, \gamma_t, \beta] = \int y \exp [\ell_{it}(\beta, \, \pi_{it} )] dy, $$ and the partial effects are usually defined as differences or derivatives of this conditional expectation with respect to the components of $X_{it}$. For example, the partial effects for the binary response and Poisson models described in Section (ref) satisfy this restriction.

The distribution of the unobserved individual and time effects is not ancillary for the APEs, unlike for model parameters. We therefore need to make assumptions on this distribution to define and interpret the APEs, and to derive the asymptotic distribution of their estimators. We control the heterogeneity of the partial effects assuming that the individual effects and explanatory variables are identically distributed cross sectionally and/or stationary over time. If $(X_{it}, \alpha_i, \gamma_t)$ is identically distributed over $i$ and can be heterogeneously distributed over $t$, $\mathbb{E}[\Delta_{it}] = \delta_t^0$ and $\delta_{NT}^0 = T^{-1} \sum_{t=1}^T \delta_t^0$ changes only with $T$. If $(X_{it}, \alpha_i, \gamma_t)$ is stationary over $t$ and can be heterogeneously distributed over $i$, $\mathbb{E}[\Delta_{it}] = \delta_i^0$ and $\delta_{NT}^0 = N^{-1} \sum_{i=1}^N \delta_i^0$ changes only with $N$. Finally, if $(X_{it}, \alpha_i, \gamma_t)$ is identically distributed over $i$ and stationary over $t$, $\mathbb{E}[\Delta_{it}] = \delta_{NT}^0$ and $\delta_{NT}^0 = \delta^0$ does not change with $N$ and $T.$ We also impose smoothness and moment conditions on the function $\Delta$ that defines the partial effects. We use these conditions to derive higher-order stochastic expansions for the fixed effect estimator of the APEs and to bound the remainder terms in these expansions. Let $\{\alpha_i\}_N := \{\alpha_i : 1 \leq i \leq N\}$, $\{\gamma_t\}_T := \{\gamma_t : 1 \leq t \leq T\},$ and $\{X_{it}, \alpha_i, \gamma_t \}_{NT} := \{(X_{it}, \alpha_i, \gamma_t) : 1 \leq i \leq N, 1 \leq t \leq T\}.$

Analogous to $\Xi_{it}$ and $D_{\beta \pi^q} \ell_{it} $ in equation (ref) we define

for $q \in \{1,2\}$. Here, $\Psi_{it}$ is the population projection of $ \partial_{\pi} \Delta_{it} / \mathbb{E}_{\phi}[ \partial_{\pi^2} \ell_{it} ]$ on the space spanned by the incidental parameters under the metric given by $ \mathbb{E}_{\phi}[ - \partial_{\pi^2} \ell_{it} ]$. We use analogous notation to the previous section for the derivatives with respect to $\beta$ and higher order derivatives with respect to $\pi$.

Let $\delta_{NT}^0$ and $\widehat \delta$ be the APE and its fixed effects estimator, defined as in equations (ref) and (ref) with $\Delta(X_{it}, \beta, \alpha_i, \gamma_t) = \Delta_{it}(\beta, \alpha_i + \gamma_t).$\footnote{We keep the dependence of $\delta_{NT}^0$ on $NT$ to distinguish $\delta_{NT}^0$ from $\delta^0 = \lim_{N,T \to \infty} \delta_{NT}^0$.} The following theorem establishes the asymptotic distribution of $\widehat \delta.$

In the following examples we assume that the APEs are constructed from asymptotically unbiased estimators of the model parameters.

Example (ref) (Binary response model). Consider the partial effects defined in ((ref)) and ((ref)) with $$ \Delta_{it}(\beta , \pi) = F(\beta_k + X_{it,-k}'\beta_{-k} + \pi) - F(X_{it,-k}'\beta_{-k} + \pi) \text{ and } \Delta_{it}(\beta, \pi) = \beta_k \partial F(X_{it}'\beta + \pi). $$ Using the notation previously introduced for this example, the components of the asymptotic bias of $\widetilde \delta$ are

where $\tilde \Psi_{it}$ is the residual of the population regression of $ - \partial_{\pi} \Delta_{it} / \mathbb{E}_{\phi}[ \omega_{it}]$ on the space spanned by the incidental parameters under the metric given by $\mathbb{E}_{\phi}[\omega_{it}]$. If all the components of $X_{it}$ are strictly exogenous, the first term of $ \overline B_{\infty}^{\delta}$ is zero.}

Example (ref) (Poisson model). Consider the partial effect $$ \Delta_{it}(\beta, \pi) = g_{it}(\beta) \exp(X_{it}'\beta + \pi), $$ where $g_{it}$ does not depend on $\pi$. For example, $g_{it}(\beta) = \beta_k + \beta_j h(Z_{it})$ in ((ref)). Using the notation previously introduced for this example, the components of the asymptotic bias are

and $\overline D_{\infty}^{\delta} = 0$, where $\tilde g_{it}$ is the residual of the population projection of $g_{it}$ on the space spanned by the incidental parameters under a metric weighted by $\mathbb{E}_{\phi}[\omega_{it}]$. The asymptotic bias is zero if all the components of $X_{it}$ are strictly exogenous or $g_{it}(\beta)$ is constant. The latter arises in the leading case of the partial effect of the $k$-th component of $X_{it}$ since $g_{it}(\beta) = \beta_k$. This no asymptotic bias result applies to any type of regressor, strictly exogenous or predetermined.}

Bias corrected estimators

The results of the previous sections show that the asymptotic distributions of the fixed effects estimators of the model parameters and APEs can have biases of the same order as the variances under sequences where $T$ grows at the same rate as $N$. This is the large-$T$ version of the incidental parameters problem that invalidates any inference based on the fixed effect estimators even in large samples. In this section we describe how to construct analytical and jackknife bias corrections for the fixed effect estimators and give conditions for the asymptotic validity of these corrections.

The jackknife correction for the model parameter $\beta$ in equation (ref) is generic and applies to the panel model. For the APEs, the jackknife correction is formed similarly as

where $\widetilde{\delta}_{N,T/2}$ is the average of the 2 split jackknife estimators of the APE that use all the individuals and leave out the first and second halves of the time periods, and $\widetilde{\delta}_{N/2,T}$ is the average of the 2 split jackknife estimators of the APE that use all the time periods and leave out half of the individuals.

The analytical corrections are constructed using sample analogs of the expressions in Theorems (ref) and (ref), replacing the true values of $\beta$ and $\phi$ by the fixed effects estimators. To describe these corrections, we introduce some additional notation. For any function of the data, unobserved effects and parameters $g_{itj}(\beta,\alpha_i + \gamma_t,\alpha_i + \gamma_{t-j})$ with $0 \leq j < t$, let $\widehat g_{itj} = g_{it}(\widehat \beta, \widehat \alpha_i + \widehat \gamma_t, \widehat \alpha_i + \widehat \gamma_{t-j})$ denote the fixed effects estimator, e.g., $\widehat{\mathbb{E}_{\phi}[\partial_{\pi^2} \ell_{it}]}$ denotes the fixed effects estimator of $\mathbb{E}_{\phi}[\partial_{\pi^2} \ell_{it}].$ Let $\widehat{\cal H}^{-1}_{(\alpha\alpha)}$, $\widehat{\cal H}^{-1}_{(\alpha\gamma)}$, $\widehat{\cal H}^{-1}_{(\gamma\alpha)}$, and $\widehat{\cal H}^{-1}_{(\gamma\gamma)}$ denote the blocks of the matrix $\widehat{\cal H}^{-1}$, where $$ \widehat{\cal H} = \left(

\right) + \frac{b} {\sqrt{NT}} \, vv' ,$$ $\widehat{\mathcal{H}}_{(\alpha\alpha)}^* = diag( - \sum_{t} \widehat{\mathbb{E}_{\phi}[\partial_{\pi^2} \ell_{it}]})/\sqrt{NT}$, $\widehat{\mathcal{H}}_{(\alpha\alpha)} ^*= diag( - \sum_{i} \widehat{\mathbb{E}_{\phi}[\partial_{\pi^2} \ell_{it}]})/\sqrt{NT}$, and $\widehat{\mathcal{H}}_{(\alpha\gamma)it}^* = - \widehat{\mathbb{E}_{\phi}[\partial_{\pi^2} \ell_{it}]}/\sqrt{NT}$. Let

The $k$-th component of $\widehat \Xi_{it}$ corresponds to a least squares regression of $ \widehat{\mathbb{E}_{\phi} \left( \partial_{\beta_k \pi} \ell_{it} \right)}/\widehat{\mathbb{E}_{\phi}( \partial_{\pi^2} \ell_{it} )}$ on the space spanned by the incidental parameters weighted by $ \widehat{\mathbb{E}_{\phi}( - \partial_{\pi^2} \ell_{it})}.$

The analytical bias corrected estimator of $\beta^0$ is

where $\widehat B_{NT}^{\beta} = \widehat W^{-1} \widehat B$, $\widehat D_{NT}^{\beta} = \widehat W^{-1} \widehat D$,

and $L$ is a trimming parameter for estimation of spectral expectations such that $L \to \infty$ and $L/T \to 0$ HahnKuersteiner2011. Here we use truncation instead of kernel smoothing in the estimation of spectral expectations following Hahn and Kuersteiner HK2007. Note that, unlike for variance estimation, a kernel is not needed to ensure that the bias estimator be positive. Instead of choosing a value of $L$, our recommendation for practice is to conduct a sensitivity analysis by reporting estimates for multiple values of $L$ starting from $L=1$. From our experience based on extensive Monte Carlo simulations, we do not recommend values of $L$ greater than $4$, because the finite-sample dispersion of the estimator quickly increases with $L$. We refer to Section (ref) for an example of sensitivity analysis with respect to $L$. The factor $T/(T-j)$ is a degrees of freedom adjustment that rescales the time series averages $T^{-1} \sum_{t=j+1}^T$ by the number of observations instead of by $T$. Similar corrections for conditional mean models can be formed using the sample analogs of the expressions of $\overline{B}_{\infty}$ and $\overline{D}_{\infty}$ in Remark (ref). We do not spell out these estimators for the sake of brevity.

Asymptotic $(1-p)$--confidence intervals for the components of $\beta^0$ can be formed as $$ \widetilde \beta_k^A \pm z_{1-p} \sqrt{\widehat W_{kk}^{-1} / (NT)}, \ \ k = \{1, ..., \dim \beta^0\}, $$ where $z_{1-p}$ is the $(1-p)$--quantile of the standard normal distribution, and $\widehat W_{kk}^{-1} $ is the $(k,k)$-element of the matrix $\widehat W^{-1}$. In conditional moment models we replace $\widehat W_{kk}$ by the $(k,k)$-element of the matrix $\widehat W^{-1} \widehat \Omega \widehat W^{-1}$, where $$ \widehat \Omega = \frac 1 {NT} \sum_{i=1}^N \sum_{t=1}^T \sum_{\tau=1}^T \widehat{\mathbb{E}_{\phi} \left[ D_{\beta} \ell_{it} (D_{\beta} \ell_{i\tau})' \right]}. $$

We have implemented the analytical correction at the level of the estimator. Alternatively, we can implement the correction at the level of the score or first order conditions by solving

for $\beta$. Global concavity of the objective function guarantees that the solution to (ref) is unique. Other possible extensions such as continuously updated score corrections where $\overline B_{\infty}$ and $\overline D_{\infty}$ are estimated together with $\beta$, corrections at the level of the objective function, or iterative corrections are left to future research.

The analytical bias corrected estimator of $\delta^0_{NT}$ is $$ \widetilde \delta^A = \widehat \delta - \widehat B^\delta/T - \widehat D^\delta/N, $$ where $\widetilde \delta$ is the APE constructed from a bias corrected estimator of $\beta$. Let

The fixed effects estimators of the components of the asymptotic bias are

The estimator of the asymptotic variance depends on the sampling properties of the unobserved effects. Under the independence assumption of Remark (ref) with all the components of $X_{it}$ strictly exogenous,

where $ \widehat{\tilde \Delta}_{it} = \widehat \Delta_{it} - N^{-1} \sum_{i=1}^N \widehat \Delta_{it}$ under identical distribution over $i$, $ \widehat{\tilde \Delta}_{it} = \widehat \Delta_{it} - T^{-1} \sum_{t=1}^T \widehat \Delta_{it}$ under stationarity over $t$, and $ \widehat{\tilde \Delta}_{it} = \widehat \Delta_{it} - \widehat \delta$ under both. Note that we do not need to specify the convergence rate $r_{NT}$ to make inference because the standard errors $\sqrt{\widehat{V}^{\delta}}/r_{NT}$ do not depend on $r_{NT}$. Bias corrected estimators and confidence intervals can be constructed in the same fashion as for the model parameter.

We use the following homogeneity assumption to show the validity of the jackknife corrections for the model parameters and APEs. It implies that $ \widetilde{\beta}_{N,T/2} - \widehat \beta_{NT} = \overline B_{\infty}^{\beta}/T + o_P(T^{-1} \vee N^{-1})$ and $\widetilde{\beta}_{N/2,T} - \widehat \beta_{NT} = \overline D_{\infty}^{\beta}/N + o_P(T^{-1} \vee N^{-1})$, which are weaker but higher level sufficient conditions for the validity of the jackknife for the model parameter. For APEs, Assumption (ref) also ensures that these effects do not change with $T$ and $N$, i.e. $\delta_{NT}^0 = \delta^0$. The analytical corrections {\it do not} require this assumption.

This assumption might seem restrictive for dynamic models where $X_{it}$ includes lags of the dependent variable because in this case it restricts the unconditional distribution of the initial conditions of $Y_{it}$. Note, however, that Assumption (ref) allows the initial conditions to depend on the unobserved effects. In other words, it does not impose that the initial conditions are generated from the stationary distribution of $Y_{it}$ conditional on $X_{it}$ and $\phi$. Assumption (ref) rules out time trends and structural breaks in the processes for the unobserved effects and observed variables.

The following theorems are the main result of this section. They show that the analytical and jackknife bias corrections eliminate the bias from the asymptotic distribution of the fixed effects estimators of the model parameters and APEs without increasing variance, and that the estimators of the asymptotic variances are consistent.

Monte Carlo Experiments

This section reports evidence on the finite sample behavior of fixed effects estimators of model parameters and APEs in static models with strictly exogenous regressors and dynamic models with predetermined regressors such as lags of the dependent variable. We analyze the performance of uncorrected and bias-corrected fixed effects estimators in terms of bias and inference accuracy of their asymptotic distribution. In particular we compute the biases, standard deviations, and root mean squared errors of the estimators, the ratio of average standard errors to the simulation standard deviations (SE/SD); and the empirical coverages of confidence intervals with 95% nominal value (p; .95).\footnote{The standard errors are computed using the expressions ((ref)) and ((ref)) with $ \widehat{\tilde \Delta}_{it} = \widehat \Delta_{it} - \widehat \delta$, evaluated at uncorrected estimates of the parameters. We find little difference in performance of constructing standard errors based on corrected estimates.} Overall, we find that the analytically corrected estimators dominate the uncorrected and jackknife corrected estimators.\footnote{Kristensen and Salani{\'e} KS2013 also found that analytical corrections dominate jackknife corrections to reduce the bias of approximate estimators.} A possible explanation for the better finite-sample performance of the analytical over the jackknife corrections is that the jackknife increases dispersion because the components of the bias are estimated from subsamples that include half of the observations of the panel. We observe this variance increase in all our numerical examples, specially in short panels. The jackknife corrections are also more sensitive than the analytical corrections to Assumption (ref). All the results are based on 500 replications. The designs correspond to static and dynamic probit models. As in the analytical example of Section (ref), we find that our large $T$ asymptotic approximations capture well the behavior of the fixed effects estimator and the bias corrections in moderately long panels with $N=56$ and $T=14$.

Static probit model

The data generating process is

where $\alpha_{i} \sim \mathcal{N}(0,1/16)$, $\gamma_{t} \sim \mathcal{N}(0,1/16)$, $\varepsilon_{it} \sim \mathcal{N}(0,1)$, and $\beta = 1$. We consider two alternative designs for $X_{it}$: autoregressive process and linear trend process both with individual and time effects. In the first design, $X_{it} = X_{i,t-1} / 2 + \alpha_{i} + \gamma_{t} + \upsilon_{it}$, $\upsilon_{it} \sim \mathcal{N}(0,1/2)$, and $X_{i0} \sim \mathcal{N}(0,1)$. In the second design, $X_{it} = 2 t / T + \alpha_{i} + \gamma_{t} + \upsilon_{it}$, $\upsilon_{it} \sim \mathcal{N}(0,3/4)$, which violates Assumption (ref). In both designs $X_{it}$ is strictly exogenous with respect to $\varepsilon_{it}$ conditional on the individual and time effects. The variables $\alpha_i$, $\gamma_t$, $\varepsilon_{it}$, $\upsilon_{it}$, and $X_{i0}$ are independent and $i.i.d.$ across individuals and time periods. We generate panel data sets with $N=56$ individuals and three different numbers of time periods $T$: 14, 28 and 56.\footnote{Following a suggestion from an anonymous referee, we obtained results for panel data sets with $T=56$ and $N$ in $\{14, 28, 56\}$. These results are similar to the results reported and are available from the authors upon request.}

Table 3 reports the results for the probit coefficient $\beta$, and the APE of $X_{it}$. We compute the APE using ((ref)). Throughout the table, MLE-FETE corresponds to the probit maximum likelihood estimator with individual and time fixed effects, Analytical is the bias corrected estimator that uses the analytical correction, and Jackknife is the bias corrected estimator that uses SPJ in both the individual and time dimensions. The cross-sectional division in the jackknife follows the order of the observations. All the results are reported in percentage of the true parameter value.

We find that the bias is of the same order of magnitude as the standard deviation for the uncorrected estimator of the probit coefficient causing severe undercoverage of the confidence intervals. This result holds for both designs and all the sample sizes considered. The bias corrections, specially Analytical, remove the bias without increasing dispersion, and produce substantial improvements in rmse and coverage probabilities. For example, Analytical reduces rmse by 50% and increases coverage by 26% in the first design with $T=14$. As in Hahn and Newey Hahn:2004p882 and Fernandez-Val FernandezVal:2009p3313, we find very little bias in the uncorrected estimates of the APE, despite the large bias in the probit coefficients. Jackknife performs relatively worse in the second design that does not satisfy Assumption (ref).

Dynamic probit model

The data generating process is

where $\alpha_{i} \sim \mathcal{N}(0,1/16)$, $\gamma_{t} \sim \mathcal{N}(0,1/16)$, $\varepsilon_{it} \sim \mathcal{N}(0,1)$, $\beta_Y = 0.5$, and $\beta_Z = 1$. We consider two alternative designs for $Z_{it}$: autoregressive process and linear trend process both with individual and time effects. In the first design, $Z_{it} = Z_{i,t-1} / 2 + \alpha_{i} + \gamma_{t} + \upsilon_{it}$, $\upsilon_{it} \sim \mathcal{N}(0,1/2)$, and $Z_{i0} \sim \mathcal{N}(0,1)$. In the second design, $Z_{it} = 1.5 t / T + \alpha_{i} + \gamma_{t} + \upsilon_{it}$, $\upsilon_{it} \sim \mathcal{N}(0,3/4)$, which violates Assumption (ref). The variables $\alpha_i$, $\gamma_t$, $\varepsilon_{it}$, $\upsilon_{it}$, and $Z_{i0}$ are independent and $i.i.d.$ across individuals and time periods. We generate panel data sets with $N=56$ individuals and three different numbers of time periods $T$: 14, 28 and 56.

Table 4 reports the simulation results for the probit coefficient $\beta_Y$ and the APE of $Y_{i,t-1}$. We compute the partial effect of $Y_{i,t-1}$ using the expression in equation ((ref)) with $X_{it,k} = Y_{i,t-1}$. This effect is commonly reported as a measure of state dependence for dynamic binary processes. Table 5 reports the simulation results for the estimators of the probit coefficient $\beta_Z$ and the APE of $Z_{it}$. We compute the partial effect using ((ref)) with $X_{it,k} = Z_{it}$. Throughout the tables, we compare the same estimators as for the static model. For the analytical correction we consider two versions, Analytical (L=1) sets the trimming parameter to estimate spectral expectations $L$ to one, whereas Analytical (L=2) sets $L $ to two.\footnote{In results not reported for brevity, we find little difference in performance of increasing the trimming parameters to $L=3$ and $L=4$. These results are available from the authors upon request.} Again, all the results in the tables are reported in percentage of the true parameter value.

The results in table 4 show important biases toward zero for both the probit coefficient and the APE of $Y_{i,t-1}$ in the two designs. This bias can indeed be substantially larger than the corresponding standard deviation for short panels yielding coverage probabilities below 70% for $T=14$. The analytical corrections significantly reduce biases and rmse, bring coverage probabilities close to their nominal level, and have little sensitivity to the trimming parameter $L$. The jackknife corrections reduce bias but increase dispersion, producing less drastic improvements in rmse and coverage than the analytical corrections. The results for the APE of $Z_{it}$ in table 5 are similar to the static probit model. There are significant bias and undercoverage of confidence intervals for the coefficient $\beta_Z$, which are removed by the corrections, whereas there are little bias and undercoverage in the APE. As in the static model, Jackknife performs relatively worse in the second design.

Concluding remarks

In this paper we develop analytical and jackknife corrections for fixed effects estimators of model parameters and APEs in semiparametric nonlinear panel models with additive individual and time effects. Our analysis applies to conditional maximum likelihood estimators with concave log-likelihood functions, and therefore covers logit, probit, ordered probit, ordered logit, Poisson, negative binomial, and Tobit estimators, which are the most popular nonlinear estimators in empirical economics.

We are currently developing similar corrections for nonlinear models with interactive individual and time effects (Chen, Fern\'andez-Val, and Weidner CFW2014). Another interesting avenue of future research is to derive higher-order expansions for fixed effects estimators with individual and time effects. These expansions are needed to justify theoretically the validity of alternative corrections based on the leave-one-observation-out panel jackknife method of Hahn and Newey Hahn:2004p882.

\setcounter{section}{0} \setcounter{footnote}{0} \setcounter{page}{1}

\abstract{ This supplemental material contains five appendices. Appendix (ref) presents the results of an empirical application and a Monte Carlo simulation calibrated to the application. Following Aghion et al. AghionBloomBlundellGriffithHowitt2005, we use a panel of U.K. industries to estimate Poisson models with industry and time effects for the relationship between innovation and competition. Appendix (ref) gives the proofs of Theorems (ref) and (ref). Appendices (ref), (ref), and (ref) contain the proofs of Appendices (ref), (ref), and (ref), respectively. Appendix (ref) collects some useful intermediate results that are used in the proofs of the main results.}

Relationship between Innovation and Competition

Empirical Example

To illustrate the bias corrections with real data, we revisit the empirical application of Aghion, Bloom, Blundell, Griffith and Howitt AghionBloomBlundellGriffithHowitt2005 (ABBGH) that estimated a count data model to analyze the relationship between innovation and competition. They used an unbalanced panel of seventeen U.K. industries followed over the 22 years between 1973 and 1994.\footnote{We assume that the observations are missing at random conditional on the explanatory variables and unobserved effects and apply the corrections without change since the level of attrition is low in this application.} The dependent variable, $Y_{it},$ is innovation as measured by a citation-weighted number of patents, and the explanatory variable of interest, $Z_{it},$ is competition as measured by one minus the Lerner index in the industry-year.

Following ABBGH we consider a quadratic static Poisson model with industry and year effects where $$ Y_{it} \mid Z_i^T, \alpha_i,\gamma_t \sim \mathcal{P}(\exp[\beta_1 Z_{it} + \beta_2 Z_{it}^2 + \alpha_i + \gamma_t]), $$ for $(i = 1,...,17; t = 1973, ..., 1994),$ and extend the analysis to a dynamic Poisson model with industry and year effects where $$ Y_{it} \mid Y_i^{t-1}, Z_i^t, \alpha_i, \gamma^t \sim \mathcal{P}(\exp[\beta_{Y} \log(1 + Y_{i,t-1}) + \beta_1 Z_{it} + \beta_2 Z_{it}^2 + \alpha_i + \gamma_t]), $$ for $(i = 1,...,17; t = 1974, ..., 1994).$ In the dynamic model we use the year 1973 as the initial condition for $Y_{it}$.

Table S1 reports the results of the analysis. Columns (2) and (3) for the static model replicate the empirical results of Table I in ABBGH (p. 708), adding estimates of the APEs. Columns (4) and (5) report estimates of the analytical corrections that do not assume that competition is strictly exogenous with $L=1$ and $L=2$, and column (6) reports estimates of the jackknife bias corrections described in equation ((ref)) of the paper. Note that we do not need to report separate standard errors for the corrected estimators, because the standard errors of the uncorrected estimators are consistent for the corrected estimators under the asymptotic approximation that we consider.\footnote{In numerical examples, we find very little gains in terms of the ratio SE/SD and coverage probabilities when we reestimate the standard errors using bias corrected estimates.} Overall, the corrected estimates, while numerically different from the uncorrected estimates in column (3), agree with the inverted-U pattern in the relationship between innovation and competition found by ABBGH. The close similarity between the uncorrected and bias corrected estimates gives some evidence in favor of the strict exogeneity of competition with respect to the innovation process.

\setcounter{table}{1}

The results for the dynamic model show substantial positive state dependence in the innovation process that is not explained by industry heterogeneity. Uncorrected fixed effects underestimates the coefficient and APE of lag patents relative to the bias corrections, specially relative to the jackknife. The pattern of the differences between the estimates is consistent with the biases that we find in the numerical example in Table S4. Accounting for state dependence does not change the inverted-U pattern, but flattens the relationship between innovation and competition.

Table (ref) implements Chow-type homogeneity tests for the validity of the jackknife corrections. These tests compare the uncorrected fixed effects estimators of the common parameters within the elements of the cross section and time series partitions of the panel. Under time homogeneity, the probability limit of these estimators is the same, so that a standard Wald test can be applied based on the difference of the estimators in the sub panels within the partition. For the static model, the test is rejected at the 1% level in both the cross section and time series partitions. Since the cross sectional partition is arbitrary, these rejection might be a signal of model misspecification. For the dynamic model, the test is rejected at the 1% level in the time series partition, but it cannot be rejected at conventional levels in the cross section partition. The rejection of the time homogeneity might explain the difference between the jackknife and analytical corrections in the dynamic model.

Calibrated Monte Carlo Simulations

We conduct a simulation that mimics the empirical example. The designs correspond to static and dynamic Poisson models with additive individual and time effects. We calibrate all the parameters and exogenous variables using the dataset from ABBGH.

Static Poisson model

The data generating process is

where $\mathcal{P}$ denotes the Poisson distribution. The variable $Z_{it}$ is fixed to the values of the competition variable in the dataset and all the parameters are set to the fixed effect estimates of the model. We generate unbalanced panel data sets with $T=22$ years and three different numbers of industries $N$: 17, 34, and 51. In the second (third) case, we double (triple) the cross-sectional size by merging two (three) independent realizations of the panel.

Table S3 reports the simulation results for the coefficients $\beta_1$ and $\beta_2$, and the APE of $Z_{it}$. We compute the APE using the expression ((ref)) with $H(Z_{it}) = Z_{it}^2$. Throughout the table, MLE corresponds to the pooled Poisson maximum likelihood estimator (without individual and time effects), MLE-TE corresponds to the Poisson estimator with only time effects, MLE-FETE corresponds to the Poisson maximum likelihood estimator with individual and time fixed effects, Analytical (L=l) is the bias corrected estimator that uses the analytical correction with $L=l$, and Jackknife is the bias corrected estimator that uses SPJ in both the individual and time dimensions. The analytical corrections are different from the uncorrected estimator because they do not use that the regressor $Z_{it}$ is strictly exogenous. The cross-sectional division in the jackknife follows the order of the observations. The choice of these estimators is motivated by the empirical analysis of ABBGH. All the results in the table are reported in percentage of the true parameter value.

The results of the table agree with the no asymptotic bias result for the Poisson model with exogenous regressors. Thus, the bias of MLE-FETE for the coefficients and APE is negligible relative to the standard deviation and the coverage probabilities get close to the nominal level as $N$ grows. The analytical corrections preserve the performance of the estimators and have very little sensitivity to the trimming parameter. The jackknife correction increases dispersion and rmse, specially for the small cross-sectional size of the application. The estimators that do not control for individual effects are clearly biased.

Dynamic Poisson model

The data generating process is

The competition variable $Z_{it}$ and the initial condition for the number of patents $Y_{i0}$ are fixed to the values in the dataset and all the parameters are set to the fixed effect estimates of the model. To generate panels, we first impute values to the missing observations of $Z_{it}$ using forward and backward predictions from a panel AR(1) linear model with individual and time effects. We then draw panel data sets with $T=21$ years and three different numbers of industries $N$: 17, 34, and 51. As in the static model, we double (triple) the cross-sectional size by merging two (three) independent realizations of the panel. We make the generated panels unbalanced by dropping the values corresponding to the missing observations in the original dataset.

Table S4 reports the simulation results for the coefficient $\beta^0_Y$ and the APE of $Y_{i,t-1}$. The estimators considered are the same as for the static Poisson model above. We compute the partial effect of $Y_{i,t-1}$ using ((ref)) with $Z_{it} = Y_{i,t-1}$, $H(Z_{it}) = \log (1+Z_{it}),$ and dropping the linear term. Table S5 reports the simulation results for the coefficients $\beta_1^0$ and $\beta_2^0$, and the APE of $Z_{it}$. We compute the partial effect using ((ref)) with $H(Z_{it}) = Z_{it}^2$. Again, all the results in the tables are reported in percentage of the true parameter value.

The results in table S4 show biases of the same order of magnitude as the standard deviation for the fixed effects estimators of the coefficient and APE of $Y_{i,t-1}$, which cause severe undercoverage of confidence intervals. Note that in this case the rate of convergence for the estimator of the APE is $r_{NT} = \sqrt{NT}$, because the individual and time effects are hold fixed across the simulations. The analytical corrections reduce bias by more than half without increasing dispersion, substantially reducing rmse and bringing coverage probabilities closer to their nominal levels. The jackknife corrections reduce bias and increase dispersion leading to lower improvements in rmse and coverage probability than the analytical corrections. The results for the coefficient of $Z_{it}$ in table 8 are similar to the static model. The results for the APE of $Z_{it}$ are imprecise, because the true value of the effect is close to zero.

Proofs of Theorems (ref) and (ref)

We start with a lemma that shows the consistency of the fixed effects estimators of averages of the data and parameters. We will use this result to show the validity of the analytical bias corrections and the consistency of the variance estimators.