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.
101,648 characters · 23 sections · 23 citation commands
Fixed Effect Estimation of Large T Panel Data Models
{
{ Keywords:} panel data, fixed effects, incidental parameter problem, bias correction, unobserved heterogeneity, nonlinear models, jackknife. }
One of the main advantages of panel data over cross-sectional or time series data is the possibility of accounting for multiple sources of unobserved heterogeneity. This is accomplished by including individual and time unobserved effects into the model, which control for any unobserved covariate that is either time or cross-sectional invariant. In this article we consider fixed effect approaches to panel data models, where no distributional assumption on the unobserved effects is imposed, therefore allowing the unobserved effects to be arbitrarily related with the observed covariates. Fixed effects contrast with random effects approaches that impose restrictions on the distribution of the unobserved effects conditional on the observed covariates. We refer to \citet*{Arellano2003}, \citet*{Baltagi2008}, \citet*{Hsiao2015}, and \citet*{Wooldridge2010} for modern textbook treatments on the difference between fixed and random effects.
We consider semiparametric models that specify parametrically some characteristic of the distribution of the outcome variable of interest conditional on the observed covariates and unobserved effects. Examples include linear and nonlinear regression models for the conditional expectation, and generalized linear, distribution and quantile regression models for the conditional distribution. Generalized linear models include the most commonly used nonlinear models such as probit, logit, Poisson, negative binomial, proportional hazard, and tobit models. The nonparametric or unspecified part of the model consists of the joint distribution of the unobserved effects, exogenous covariates and initial conditions. We do not review nonparametric panel models such as the nonseparable models considered in \citet*{AltonjiMatzkin2005}, \citet*{Evdokimov2010}, \citet*{GrahamPowell2012}, \citet*{HoderleinWhite2012}, \citet*{CFHN13}, \citet*{CFHHN15}, \citet*{CFN17}, \citet*{Freyberger2017}, and \citet*{Torgovitsky2016}, among others. We refer to \citet*{Matzkin2007} for an excellent survey on nonparametric identification of nonseparable models including results for panel data.
We analyze the properties of fixed effects estimators of model parameters. These estimators treat the unobserved effects as parameters to be estimated. In linear models with strictly exogenous covariates and individual effects, fixed effects is numerically equivalent to the within-group estimator that removes the individual effects by taking differences within each individual. We also analyze fixed effects estimation of average partial effects (APEs), which are averages of functions of the data, parameters and unobserved effects. APEs are often the quantities of interest in nonlinear models because they correspond to marginal effects of the covariates in some characteristic of the distribution of the outcome conditional on covariates and unobserved effects, averaged over the observed and unobserved heterogeneity. The main challenge with fixed effects estimators is to deal with the incidental parameter problem coming from the estimation of the unobserved effects. The estimated unobserved effects, usually called fixed effects, can be very noisy because there are few observations that are informative about them. Section (ref) describes the incidental parameter problem in detail.
Motivated by the first panel data sets available, the econometric theory was initially designed for panels where the number of time periods, $T$, is small compared to the number of cross-sectional units, $N$. This theory analyzed identification, estimation and inference using a fixed-$T$ asymptotic approximation that considers sequences of panels with increasing $N$ and fixed $T$. Consistent approaches under this approximation usually consist of removing the unobserved effects from the model by some smart transformation such as differencing or conditioning on sufficient statistics. Examples include \citet*{Mundlak1978}, \citet*{ArellanoBond1991} and \citet*{Moreira:2009} for linear models; \citet*{Cox1958}, \citet*{Rasch1960}, \citet*{Andersen1970}, \citet*{Chamberlain1980}, \citet*{Manski1987}, \citet*{Horowitz1992}, \citet*{HonoreKyriazidou2000}, and \citet*{ArellanoCarrasco2003} for static and dynamic binary choice models; \citet*{HausmanHallGriliches1984} for static Poisson models; \citet*{Honore1992} for tobit models; \citet*{Kyriazidou1997,Kyriazidou2001} for sample selection models; and \citet*{Muris2016} and \citet*{BotosaruMuris2017} for ordered logit and probit models; \citet*{Bonhomme2012} provided a unifying scheme for many of those methods.
The fixed-$T$ approach has encountered several limitations. First, transformations to remove the unobserved effects are not available for some models such as probit and quantile regression. This issue gets worse in models with unobserved effects in multiple dimensions, where even if sufficient statistics exist conditioning on them can be computationally challenging \citep*{hirjietal87, Charbonneau2014}. Second, differencing over time in dynamic linear models does not work well when the data are very persistent because it removes most of the information about the parameters. Third, removing the unobserved effects from the model precludes the estimation of APEs in nonlinear models. Fourth, under fixed-$T$ asymptotics the incidental parameter problem of fixed effects estimation is a consistency problem, which is difficult to tackle. Fifth, some models are not $\sqrt{N}$-estimable or not even point identified under the fixed-$T$ approximation, see , e.g., \citet*{Chamberlain2010}, \citet*{HonoreTamer2006}, \citet*{CFHN13}, \citet*{ShiShumSong2014}, and \citet*{PakesPorter2013}.
All the limitations of the fixed-$T$ approximation mentioned in the previous paragraph, together with the recent availability of long panel datasets, has led to an alternative asymptotic approximation to panel data that considers sequences of panels where both $N$ and $T$ increase. This large-$T$ approximation, developed by \citet*{Phillips:1999p733} for linear models, deals with some of the shortcomings of the fixed-$T$ approximation. For example, most models are point identified with large-$T$ and the incidental parameter problem of fixed effects estimation becomes an asymptotic bias problem that is easier to tackle. Bias corrections have been developed to provide improved fixed effects estimators of model parameters and APEs. This large-$T$ approximation can also be applied naturally to other types of data with a grouping structure similar to panel data such as network and trade data.
The two alternative asymptotic approximations to panel data should be viewed as complements rather than substitutes. For example, if a fixed $T$ consistent estimator is available for the particular parameter under consideration, then it generally should be used. Otherwise it is likely that the model is not point identified for fixed-$T$. Then, one can either consider bounding the parameter if $T$ is small, or drawing on a large-$T$ consistent estimator if $T$ is moderately large with respect to $N$. \citet*{HonoreTamer2006} and \citet*{CFHN13} showed that in some binary response models the bounds for parameters and APEs can be very tight and shrink very rapidly with $T$. Thus, unless $T$ is very small, one can safely ignore that the model might only be set-identified and use large $T$ consistent estimators. Alternatively, one might entertain (correlated) random effects approaches that achieve fixed-$T$ point identification by imposing restrictions on the distribution of the unobserved effects.
This review focuses solely on fixed effects estimators of semiparametric panel models analyzed under a large-$T$ asymptotic approximation. Available survey articles on panel data estimation under fixed-$T$ approximations include \citet*{Chamberlain1984}, \citet*{ArellanoHonore2001}, \citet*{Honore2002}, \citet*{ArellanoBonhomme2011}, and \citet*{ArellanoBonhomme2017}. We concentrate mainly on the developments since \citet*{ArellanoHahn2007}, which previously reviewed the large-$T$ approach to panel data. This field has grown sufficiently fast in recent years to justify another review. In particular, \citet*{ArellanoHahn2007} considered almost exclusively models with only unobserved individual effects, whereas we consider models with unobserved effects in multiple dimensions such as individual and time effects. We also review the estimation of distributional and quantile effects with panel data, which were not discussed in \citet*{ArellanoHahn2007}. Our review relies heavily upon the important foundational work in this literature by \citet*{NeymanScott1948}, \citet*{Nickell81}, \citet*{coxreid87}, \citet*{kiviet95}, \citet*{Phillips:1999p733}, \citet*{Lancaster2000}, \citet*{Hahn:2002p717}, \citet*{Lancaster:2002p875}, \citet*{Woutersen:2002p3683}, \citet*{arellano03}, \citet*{AlvarezArellano2003}, \citet*{li2003efficiency}, \citet*{Hahn:2004p882}, \citet*{Carro:2007p3601}, \citet*{FernandezVal:2009p3313}, \citet*{HahnKuersteiner2011}, \citet*{gmm13} and \citet*{DhaeneJochmans2015}, among others.
\paragraph*{Outline:} In Section (ref) we describe the bias of maximum likelihood estimators (MLEs) in cross-sectional samples, which allows us to establish some key results in a simplified setting. Those results are then generalized to panel data models in Section (ref), where we discuss the incidental parameter bias and bias correction methods for fixed effected estimators of model parameters and average partial effects. A short heuristic derivation of the main formulas in Section (ref) is provided in the appendix. Section (ref) discusses unbalanced panels, multivariate fixed effects, distributional and quantile effects, and some further extensions of the methods, before we conclude in Section (ref).
\paragraph*{Notation:} In the following we denote by ${\mathrm{E}}$ the expectation with respect to the parametric part of the model. Averages over cross-sectional and time-series samples are denoted by ${\mathbb{E}_n} = n^{-1} \sum_{j=1}^n$, $ \mathbb{E}_N = N^{-1} \sum_{i=1}^N$, $ \mathbb{E}_T = T^{-1} \sum_{t=1}^T$, and $ \mathbb{E}_{NT} = (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T$. For example, we write $\mathbb{E}_n \ell_{j}(\theta) $ for $n^{-1} \sum_{j=1}^n \ell_{j}(\theta) $ in the next section. We also define ${\mathbb{G}_n} = \sqrt{n} \left( {\mathbb{E}_n} - {\mathrm{E}} \right) $, $\mathbb{G}_N = \sqrt{N} \left( {\mathbb{E}_N} - {\mathbb{E}_N} {\mathrm{E}} \right)$, $\mathbb{G}_T = \sqrt{T} \left( {\mathbb{E}_T} - {\mathbb{E}_T} {\mathrm{E}} \right)$, and ${\mathbb{G}_{NT}} = \sqrt{NT} \left( {\mathbb{E}_{NT}} - {\mathbb{E}_{NT}} {\mathrm{E}} \right)$. The notations ${\mathbb{E}_n}$ and ${\mathbb{G}_n}$ are used in Section (ref) for random samples, whereas $ \mathbb{E}_N$, $ \mathbb{E}_T$, $\mathbb{G}_N$, $\mathbb{G}_T$ and ${\mathbb{G}_{NT}}$ are used in Section (ref) for nonrandom panel samples, which explains the additional sample averages before the expectations ${\mathrm{E}}$. For example ${\mathbb{G}_{NT}} \ell_{it} = \sqrt{NT} {\mathbb{E}_{NT}} \left( \ell_{it} - {\mathrm{E}} \ell_{it} \right)$. In Section (ref) we often just write a bar to denote the expectation ${\mathrm{E}}$, e.g. $\overline \ell_{it} = {\mathrm{E}} \ell_{it} $. Vector-valued variables are always arranged as column-vectors in the following. Thus, we write $\theta=(\beta, \alpha, \gamma)$ to form the column-vector $\theta$ as the concatenation of the column-vectors $\beta$, $\alpha$ and $\gamma$. Vector and matrix transposition is denoted by a prime. We use $\overset{a}{\sim}$ to denote asymptotic approximation to the distribution. For example, if $\sqrt{n}(\widehat \beta - \beta_0 - B/a_n) \to_d \mathcal{N}(0,V)$ for some sequence $a_n$ such that $a_n \to \infty$ as $n \to \infty$, then we write $\widehat \beta - \beta_0 \overset{a}{\sim} \mathcal{N}(B/a_n, V/n)$.
Incidental parameters are nuisance parameters whose dimension grows with the sample size. NeymanScott1948 showed that maximum likelihood estimators (MLEs) can be asymptotically biased in models with incidental parameters. Here we provide a derivation based on a second-order asymptotic expansion that the order of the bias of the MLEs in random samples is
where $p$ is the number of parameters and $n$ is the sample size. This order corresponds with the inverse of the number of observations per parameter. In the following sections we will see that the order of the bias of fixed effects estimators of parameters and APEs can also be obtained using this simple heuristic formula in more complex panel data samples.
Given a random sample $\{z_j : 1 \leq j \leq n\}$ of a random variable $z$ that has density $f(z,\theta_0)$ with respect to some dominating measure, let
be the maximum likelihood estimator (MLE) of the value of the $p$-dimensional parameter $\theta$ that is identified by the population program
In the following expansions we assume that $\theta_0$ is uniquely determined by the first order conditions of the program (ref), $\theta \mapsto \ell_j(\theta)$ is a.s.\ differentiable to sufficient order, and other regularity conditions that allow us to bound remainder terms. We also assume that $\widehat \theta$ is consistent, i.e., $\widehat \theta \to_P \theta_0$ under some norm. Establishing this consistency might require delicate arguments in high dimensional settings where we let $p \to \infty$ as $n \to \infty$. We denote derivatives of $\theta \mapsto \ell_j(\theta)$ using superscripts $\theta$, that is, $\ell^{\theta}_{j}(\theta)$ is the $p$-dimensional gradient, $\ell^{\theta \theta}_{j}(\theta)$ is the $p \times p$ Hessian matrix, and, for $k \in \{1,\ldots,p\}$, $\ell^{\theta \theta \theta_k}_{j}(\theta)$ is the $p \times p$ matrix obtained by taking the partial derivative of the Hessian with respect to $\theta_k$. We drop the argument $\theta$ from all the functions whenever they are evaluated at $\theta_0$, e.g. $\ell^{\theta}_{j} := \ell^{\theta}_{j}(\theta_0)$.
By consistency, the asymptotic properties of $\widehat \theta$ are governed by the local properties of $\ell_{j}(\theta)$ around $\theta_0$. A second order Taylor expansion of the first order condition of (ref) around $\theta=\theta_0$ yields
We have $ \mathbb{E}_n \ell^\theta_{j} = n^{-1/2} \mathbb{G}_n \ell^\theta_{j} $, because $ {\mathrm{E}} \ell^\theta_{j} = 0$ by the first order conditions of the population program (ref). Decomposing $ \mathbb{E}_n \ell^{\theta \theta}_{j} = {\mathrm{E}} \ell^{\theta \theta}_{j} + n^{-1/2} \mathbb{G}_n \ell^{\theta \theta}_{j} $ and $ \mathbb{E}_n \ell^{\theta \theta \theta_k}_{j} = {\mathrm{E}} \ell^{\theta \theta \theta_k}_{j} + n^{-1/2} \mathbb{G}_n \ell^{\theta \theta \theta_k}_{j} $, and ignoring the term that includes $\mathbb{G}_n \ell^{\theta \theta \theta_k}_{j}$ because it is of small enough order, gives
The leading two terms of (ref) give the standard {\it first-order} approximation $\widehat \theta - \theta_0 \approx n^{-1/2} \psi_1$, where $\psi_1 := - \left( {\mathrm{E}} \ell^{\theta \theta}_{j} \right)^{-1} \mathbb{G}_n \ell^\theta_{j} $ is the influence function. Plugging this first-order approximation into the third and fourth term of (ref) and dropping low order terms yields the {\it second-order} approximation
where
We now analyze the properties of the components of $\psi_1 = (\psi_{1,1},\ldots, \psi_{1,p})$ and $\psi_2 = (\psi_{2,1},\ldots, \psi_{2,p})$. For each $k \in \{1,\ldots,p\}$, $\psi_{1,k}$ has zero mean, finite variance and asymptotic normal distribution under standard assumptions. Thus, $\psi_1$ is a variance term that does not contribute to the bias of $\widehat \theta$. To determine the mean of the quadratic term $\psi_{2,k}$, we assume for simplicity that ${\mathrm{E}} \ell^{\theta \theta}_{j} $ is a diagonal matrix.\footnote{This simplification is without loss of generality because we can diagonalize the Hessian by reparametrizing the model.} Then,
where we use the information equality $ {\mathrm{E}} \ell^{\theta_l}_{j} \ell^{\theta_m}_{j} = - {\mathrm{E}} \ell^{\theta_l \theta_m}_{j}$ to simplify the expression. This result shows that ${\mathrm{E}} \psi_{2,k}$ is proportional to $p$, because of the sum over $l = 1,\ldots,p$, and therefore $$ {\mathrm{E}} \widehat \theta - \theta_0 \approx \underbrace{ n^{-1/2} {\mathrm{E}} \psi_1 }_{=0} + n^{-1} {\mathrm{E}} \psi_2 \sim n^{-1} p, $$ which verifies the order of the bias in (ref). The term $\psi_2$ also has variance. However, the order of the variance of the term $n^{-1} \psi_2$ is lower than the order of the variance of the term $n^{-1/2} \psi_1$ in the expansion (ref).
The asymptotic bias of the MLE has consequences for inference. Confidence intervals based on the first-order asymptotic distribution might severely undercover even in large samples if the sample size is not sufficiently large relative to the dimension of the parameters, more precisely when $n = O(p^2)$. To see this, let $\rho'\theta$, $\rho \in \mathbb{R}^p$, be the inferential object of interest, and also let ${\mathrm{E}} \psi_2/p = B = (B_1, \ldots, B_p)$. By standard theory for MLE, $\rho'\psi_1 \to_d \mathcal{N}(0,\rho'\Omega\rho)$ with $\Omega = - \left( {\mathrm{E}} \ell^{\theta \theta}_{j} \right)^{-1}$, and $\rho'\psi_2/p \to_p \rho'B$, as $n\to \infty$. Then, as $n \to \infty$, $$ \rho' \left(\widehat \theta - \theta_0 \right) \overset{a}{\sim} \mathcal{N}\left( \frac{p}{ n} \, \rho'B, \; \frac 1 n \, \rho'\Omega\rho \right). $$ Consider $\rho=e_k$, the unit vector with a one in position $k$. The standard asymptotic two-sided $(1-\alpha)$-confidence interval for $\theta_k = e_k'\theta$ is $$ CI_{1-\alpha}(\theta_k) = \widehat \theta_k \pm z_{\alpha/2} \widehat \Omega^{1/2}_{kk}/\sqrt{n}, $$ where $\widehat \Omega_{kk}$ is a consistent estimator of $\Omega_{kk}$, the $(k,k)$ element of $\Omega$, and $z_{\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.
Consider $p/\sqrt{n} \to \kappa$ as $n \to \infty$. Then, the coverage of $CI_{1-\alpha}(\theta_k)$ in large samples is
if $\kappa B_k \neq 0$. This coverage can be much lower than the nominal level $1-\alpha$. For example, if $\kappa B_k/\Omega_{kk}^{1/2} = 1$, then the coverage probability of a 95% interval is less than 83% in large samples.
We observe the panel data set $\{(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 covariates $x_{it}$. The subscripts $i$ and $t$ index individuals and time periods in traditional panels, but they might index other dimensions in more general data structures such as firms and industries or siblings and families. The observations are independent across $i$ and weakly dependent across $t$.\footnote{ Long panel models with cross-sectional dependence are discussed in \citet*{Pakel2014}. } We consider the semiparametric model for each $i = 1,\ldots,N$:
where $x^t_i = (x_{i1}, \ldots, x_{it}),$ $\alpha = (\alpha_1, \ldots, \alpha_N)$, $\gamma = (\gamma_1, \ldots, \gamma_T)$, $f$ is a known density with respect to some dominating measure, and $\beta$ is a $d_{\beta}$-vector of parameters. The variables $\alpha_i$ and $\gamma_t$ are scalar unobserved individual and time effects that in economic applications capture individual heterogeneity and aggregate shocks, respectively.\footnote{We refer to \citet*{gmm13} for GMM panel models defined by moment conditions. We discuss models with multivariate individual and time effects in Section (ref).} The model is semiparametric because it does not specify the distribution of these effects nor their relationship with the covariates. The conditional density $f$ represents the parametric part of the model. The covariates $x_{it}$ are predetermined with respect to $y_{it}$ and might include lags of $y_{it}$ to accommodate dynamic models. If the covariates $x_{it}$ are strictly exogenous with respect to $y_{it}$, then $x^t_i$ can be replaced by $x^T_i=(x_{i1}, \dots, x_{iT})$ in the conditioning set. Strict exogeneity rules out dynamics and more generally any feedback from the outcome to future values of the covariates, which might be restrictive in applications.
\paragraph*{Example (i)} The simplest example of this model is the normal linear model with additive individual and time effects $y_{it} = x_{it}'\delta + \alpha_i + \gamma_t + \varepsilon_{it},$ $\varepsilon_{it} \mid x^t_i, \alpha, \gamma \sim \mathcal{N}(0,\sigma^2),$ where
\paragraph*{Example (ii)} An example of a nonlinear model is the panel binary response single index model with additive unobserved individual and time effects: $$y_{it} = 1(x_{it}'\beta + \alpha_i + \gamma_t \geq \varepsilon_{it}), \quad \varepsilon_{it} \mid x^t_i, \alpha, \gamma \sim F_{\varepsilon},$$ where $F_{\varepsilon}$ is a known CDF such as the standard normal or logistic. In this case,
\paragraph*{Example (iii)} Another example is a panel Poisson count response single index model with additive unobserved individual and time effects where, for $\lambda_{it} = \exp(x_{it}'\beta + \alpha_i + \gamma_t )$,
Let $\beta_0$, $\alpha_0 = (\alpha_{01}, \ldots, \alpha_{0N})$, and $\gamma_0 = (\gamma_{01}, \ldots, \gamma_{0T})$ denote the values of $\beta$, $\alpha$ and $\gamma$ that generate the data. We assume that these values are identified by the population conditional maximum likelihood program
where $\bar \ell_{it}(\beta, \alpha_i, \gamma_t) := {\mathrm{E}} \ell_{it}(\beta, \alpha_i, \gamma_t)$, the expected conditional log-likelihood with respect to the parametric part of the model. This program can have multiple solutions for $\alpha_0$ and $\gamma_0$. For example, in models that are additively separable in $\alpha_i$ and $\gamma_t$ such as the single index models above, $\ell_{it}(\beta, \alpha_i, \gamma_t) = \ell_{it}(\beta, \alpha_i + \gamma_t)$, the location translation $\alpha_i + c$ and $\gamma_t -c$ for a constant $c$ does not change the objective function of the program. In this case we assume that there exists a normalization that selects $\alpha_0$ and $\gamma_0$ (e.g. $\alpha_{01}=0$ or $\gamma_{01}=0$).
The fixed effects (FE) estimator treats the realizations of $\alpha$ and $\gamma$ as parameters to be estimated. It is the solution to the sample conditional maximum likelihood program, for $\widehat \alpha = (\widehat \alpha_{1}, \ldots, \widehat \alpha_{N})$ and $\widehat \gamma = (\widehat \gamma_{1}, \ldots, \widehat \gamma_{T})$,
which is the sample analog of (ref). This program can also have multiple solutions for $\widehat \alpha$ and $\widehat \gamma$. In that case we adopt the same normalization as in the population program to select $\widehat \alpha$ and $\widehat \gamma$. In the binary and count response examples, we can obtain the FE estimator using standard software routines including individual and time indicators for the unobserved effects. We discuss some computational aspects of the program (ref) in Section (ref).
The FE estimator suffers from the incidental parameter problem (NeymanScott1948, see also Lancaster2000 for a review). The incidental parameters are the individual and time fixed effects which have dimensions $N$ and $T$, respectively. Unless $T$ is large, each individual fixed effect is very noisy because only $T$ observations are informative about it. Symmetrically, unless $N$ is large, each time fixed effect is very noisy. The noise in the fixed effects generally contaminates the estimators of the other parameters. Exceptions include linear and Poisson models with strictly exogenous covariates where it is possible to separate the estimation of the fixed effects from the other parameters. The asymptotic consequences of the incidental parameter problem depend on the approximation adopted. It is a consistency problem under a fixed-$T$ or fixed-$N$ approximation, whereas it becomes an asymptotic bias problem under a large-$T$ and large-$N$ approximation, which we adopt here.
The FE program (ref) can be seem as a special case of (ref) with $j = it$, $n=NT$, $\theta=(\beta, \alpha, \gamma)$ and $p = d_{\beta} + N + T$. However, the random sampling assumption of Section (ref) is not plausible for panel data. The presence of the unobserved effects introduces heterogeneity in both dimensions, and assuming independence is often too strong when one of the dimensions is time. Some adjustments to the expansion in (ref) are therefore required. We can still derive a second-order expansion $\widehat \theta - \theta_0 \approx (NT)^{-1/2} \psi_1 + (NT)^{-1} \psi_2$ for $\widehat \theta=(\widehat \beta, \widehat \alpha, \widehat \gamma)$ and $\theta_0=(\beta_0, \alpha_0, \gamma_0)$, but the expressions of $\psi_1$ and $ \psi_2$ need to be modified to account for heterogeneity and weak serial dependence. We provide these modifications in the appendix.
To characterize the asymptotic bias of the component $\widehat \beta$ of $\widehat \theta$, it is convenient to make $\ell_{it}(\beta,\alpha,\gamma)$ information-orthogonal between $\beta$ and the rest of the parameters. This can be achieved by the transformation\footnote{ This transformation corresponds to the reparameterization $\alpha_i^* = \alpha_i - \kappa_i' \beta$ and $\gamma_t^* = \gamma_t - \rho_t' \beta$. The log-likelihood with respect to these parameters is $ \ell_{it}(\beta, \alpha_i^* + \kappa_i' \beta , \gamma_t^* + \rho_t' \beta ) =: \ell^*_{it}(\beta, \alpha_i^* , \gamma_t^* )$, which gives (ref) after renaming $(\alpha_i^* , \gamma_t^*)$ as $(\alpha_i, \gamma_t)$ again. }
where the $d_{\beta}$-vectors $\kappa_i$ and $\rho_t$ are a solution to the system of equations\footnote{ The solution for $\kappa_i$ and $\rho_t$ may not be unique. For example, if $\ell_{it}(\beta,\alpha_i,\gamma_t) =\ell_{it}(\beta,\alpha_i + \gamma_t) $, then the system does not uniquely determine $\kappa_i$ and $\rho_t$, but only $\kappa_i + \rho_t$ for all $i,t$. However, only $\kappa_i + \rho_t$ will appear in our formulas for asymptotic bias and variance of $\widehat \beta$ in that case, so the non-uniqueness of $\kappa_i$ and $\rho_t$ is not important. }
Here, as in Section (ref), we drop the arguments of the partial derivatives of $\ell_{it}$ when they are evaluated at the true values $(\beta_0, \alpha_{0i}, \gamma_{0t})$, and analogously we drop the arguments of partial derivatives of $\ell_{it}^*$ when they are evaluated at the transformed true values $(\beta_0, \alpha_{0i} - \kappa_i' \beta_0, \gamma_{0t} - \rho_t' \beta_0)$. Solving the program (ref) with $\ell^*_{it}(\beta,\alpha_i,\gamma_t)$ in place of $ \ell_{it}(\beta,\alpha_i,\gamma_t)$ does not change the solution for $\beta$. The reparametrized program cannot be used to compute the estimator because $ \ell^*_{it}(\beta,\alpha_i,\gamma_t) $ depends on $\theta_0$, but it is a convenient theoretical device to analyze the properties of $\widehat \beta$ because the Hessian $ \mathbb{E}_{NT} \overline \ell^{* \, \theta \theta}_{it} $ is block-diagonal between $\beta$ and $(\alpha,\gamma)$. We proceed by obtaining the bias from the infeasible log-likelihood $\ell_{it}^*$, and then we can express it in terms of the feasible likelihood $\ell_{it}$ using the one-to-one relationship between $\ell_{it}$ and $\ell_{it}^*$. For example, $\ell_{it}^{* \, \beta} = \ell_{it}^{\beta} + \kappa_i \ell_{it}^{\alpha_i} + \rho_t \ell_{it}^{\gamma_i}$.
Under the orthogonal parametrization, the component of the influence function $ \psi_1 = \allowbreak - ( \mathbb{E}_{NT} \overline \ell^{* \, \theta \theta}_{it})^{-1} \allowbreak {\mathbb{G}_{NT}} \ell^{*\, \theta}_{it} $ corresponding to $\beta$ simplifies to $\psi_{1,\beta} = - ( \mathbb{E}_{NT} \overline \ell^{* \, \beta \beta}_{it} )^{-1} {\mathbb{G}_{NT}} \ell^{*\, \beta}_{it}. $ Under standard conditions for MLE, as $N,T \to \infty$, $$ \psi_{1,\beta} \to_d {\cal N}(0, H^{-1} ), \qquad H = - \operatorname*{plim}_{N,T \rightarrow \infty} \mathbb{E}_{NT} \overline \ell^{*\, \beta \beta}_{it}, $$ where we use that $\psi_{1,\beta}$ is a martingale difference over $t$, and the conditional information equality. As in the cross-sectional case, $\psi_{1,\beta} $ has zero mean and determines the asymptotic variance.
The component of the second term $\psi_2$ corresponding to $\beta$ also simplifies with the reparametrization. After some calculations that are detailed in the appendix,
This second-order term is the primary source of bias. By the law of large numbers for heterogenous weakly dependent sequences as $N,T \to \infty$,
and
The expressions for $B$ and $D$ are almost symmetric with respect to the indices $i$ and $t$, but there is no double sum in $D$ because we are assuming independence of observations across $i$, while we allow for predetermined covariates over $t$. The expressions become symmetric when the covariates are strictly exogenous because the terms of the double sum $ \sum_{s=t+1}^T$ in $B$ drop out. Hence, the order of the bias is $$ \frac{\psi_{2,\beta}}{NT} = \frac{\psi^B_{2,\beta} + \psi^D_{2,\beta}}{NT} \sim \frac{B}{T} + \frac{D }{N} \sim \frac{N+T}{NT}, $$ which corresponds to the prediction from (ref) since $n=NT$ and $p = d_{\beta}+N+T \sim N+T$.
The bias term $B/T$ comes from the individual fixed effects since there are $T$ observations that are informative about each of them. Symmetrically, the bias term $D/N$ comes from the time fixed effects since there are $N$ observations that are informative about each of them. Accordingly, the term $D/N$ drops out in models with only individual effects or becomes negligible relative to $B/T$ when $N \gg T$. \citet*{Hahn:2002p717} and \citet*{AlvarezArellano2003} characterized $B$ in dynamic linear panel models with individual effects. For nonlinear models, \citet*{Hahn:2004p882} and \citet*{HahnKuersteiner2011} derived $B$ in the static and dynamic case, respectively. FernandezValWeidner2016 characterized $B$ and $D$ in static and dynamic nonlinear models with individual and time effects.
The expressions of $B$ and $D$ can be further characterized in specific models. \paragraph*{Example (i)} In the linear model (ref), $$ B = H^{-1} \operatorname*{plim}_{N,T \to \infty} - {\mathbb{E}_{NT}} \left\{ \sum_{s=t+1}^T \varepsilon_{it} \tilde x_{is}\right\}, \quad D = 0, \quad H = \operatorname*{plim}_{N,T \to \infty} {\mathbb{E}_{NT}} \left\{\tilde x_{it} \tilde x_{it}'\right\}, $$ where $\tilde x_{it}$ is the residual of the linear projection of $x_{it}$ on the space spanned by the individual and time effects (i.e. the two-way demeaned $x_{it}$ in the linear model that corresponds to the orthogonal transformation). When the covariates are strictly exogenous, $B=0$, i.e., there is no incidental parameter problem. \citet*{Hahn:2002p717} obtained this large $T$ bias expansion of the \citet*{Nickell81} bias for dynamic linear panel models with only individual effects. \citet*{hahn2006reducing} showed that the bias expression carries over to dynamic linear panel models with individual and time effects because $D=0$.
\paragraph*{Example (ii)} In the binary response model (ref) with standard normal link $F_{\varepsilon}=\Phi$ or probit model, when the covariates are strictly exogenous $$ B = H^{-1} {\mathbb{E}_{NT}} \left\{ \frac{\omega_{it} \tilde x_{it} \tilde x_{it}'}{{\mathbb{E}_T} \omega_{it}}\right\} \beta_0, \ D = H^{-1} {\mathbb{E}_{NT}} \left\{ \frac{\omega_{it} \tilde x_{it} \tilde x_{it}'}{{\mathbb{E}_N} \omega_{it}} \right\} \beta_0, \ H= {\mathbb{E}_{NT}} \left\{\omega_{it} \tilde x_{it} \tilde x_{it}'\right\}, $$ where $\tilde x_{it}$ is the residual of the linear projection of $x_{it}$ on the space spanned by the individual and time effects under a metric weighted by $\omega_{it}$, $\omega_{it} = \phi_{it}^2/[\Phi_{it} (1 - \Phi_{it})]$, and $\phi_{it}$ and $\Phi_{it}$ are the standard normal PDF and CDF evaluated at $x_{it}'\beta_0 + \alpha_{0i} + \gamma_{0t}$. In this case there is bias, which is a positive definite matrix weighted average of $\beta_0$. \citet*{FernandezVal:2009p3313} derived this result in models with individual effects and FernandezValWeidner2016 in models with individual and time effects.
\paragraph*{Example (iii)} In the Poisson model (ref), $$ B = H^{-1} \operatorname*{plim}_{N,T \to \infty} - {\mathbb{E}_{NT}} \left\{\frac{\sum_{s=t+1}^T (y_{it} - \lambda_{it}) \lambda_{is} \tilde x_{is}}{{\mathbb{E}_N} \lambda_{it}}\right\}, \ D =0, \ H = \operatorname*{plim}_{N,T \to \infty} {\mathbb{E}_{NT}} \left\{\lambda_{it} \tilde x_{it} \tilde x_{it}'\right\}, $$ where $\tilde x_{it}$ is the residual of the linear projection of $x_{it}$ on the space spanned by the individual and time effects under a metric weighted by $\lambda_{it}$. When the covariates are strictly exogenous, $B=0$, i.e., there is no incidental parameter problem. This is a well-known result for models with only individual effects (e.g., Palmgren1981),\footnote{ \citealp*{HausmanHallGriliches1984} use the conditional likelihood approach to eliminate the incidental parameters, which for the Poisson models turns out to be equivalent to fixed effect MLE, see \citet*{BlundellGriffithWindmeijer1999,BlundellGriffithWindmeijer2002} and \citet*{Lancaster:2002p875}. } which was extended to models with individual and time effect in FernandezValWeidner2016.
Finally, combining the properties of $\psi_{1,\beta}$ and $\psi_{2,\beta}$, we conclude that as $N,T \to \infty$ such that $N/T \to \kappa$, $0 < \kappa < \infty$,
This asymptotic approximation prescribes that the fixed effects estimator can have significant bias relative to its standard deviation. Moreover, by the argument given in Section (ref), confidence intervals constructed around the fixed effects estimator can severely undercover the components of $\beta_0$ even in large samples. FernandezValWeidner2016 showed that these prescriptions provide a good approximations to the behaviour of the fixed effects estimator for small sample sizes through analytical and simulation examples.
The goal of the bias corrections is to remove the bias from the asymptotic distribution (ref), ideally without increasing the variance. In other words, we want to find a bias corrected estimator $\widehat \beta_{\rm BC}$ such that as $N,T \to \infty$, $$ \widehat \beta_{\rm BC} - \beta_0 \overset{a}{\sim} {\cal N}\left( 0, \frac{H^{-1}} {NT} \right). $$ The corrected estimator therefore has small bias relative to its dispersion and the same variance as the FE estimator in large samples. Moreover, the confidence intervals constructed around the bias corrected estimator should have coverage close to their nominal level. We describe analytical and resampling methods to carry out the bias corrections. These methods rely on the asymptotic distribution (ref), together with consistent estimators of the bias.
The analytically bias corrected (ABC) estimator is $$ \widehat \beta_{\rm ABC} = \widehat \beta - \frac{\widehat B}{T} - \frac{\widehat D}{N}, $$ where $\widehat B$ and $\widehat D$ are consistent estimators of $B$ and $D$, i.e. $\widehat B \to_P B $ and $\widehat D \to_P D $ as $N,T \to \infty$. In this case, as $N,T \to \infty$ such that $N/T \to \kappa$, $0 < \kappa < \infty,$ $$ \widehat \beta_{\rm ABC} - \beta_0 \approx \frac{\psi_{1,\beta}}{\sqrt{NT}} + \frac{1}{T} \left( \frac{\psi^B_{2,\beta}}{N} - \widehat B\right) +\frac{1}{N} \left( \frac{\psi^D_{2,\beta}}{T} - \widehat D\right) \approx \frac{\psi_{1,\beta}}{\sqrt{NT}}, $$ since $\psi^B_{2,\beta}/N - \widehat B \to_P 0$ and $\psi^D_{2,\beta}/T - \widehat D \to_P 0$. Hence, $$ \widehat \beta_{\rm ABC} - \beta_0 \overset{a}{\sim} {\cal N}\left( 0, \frac{H^{-1}} {NT} \right). $$
The estimators $\widehat B$ and $\widehat D$ are constructed from the analytical expressions of $B$ and $D$ given in (ref) and (ref). Let $\widehat \ell_{it} = \ell_{it}(\widehat \beta,\widehat \alpha_i, \widehat \gamma_t)$ denote the log-likelihoods evaluated at the FE estimators, and define their derivatives evaluated at the FE estimators analogously. Let $\widehat \ell^*_{it}(\beta,\alpha_i,\gamma_t) := \ell_{it}(\beta, \alpha_i + \widehat \kappa_i' \beta , \gamma_t + \widehat \rho_t' \beta)$, with $\widehat \kappa_i$ and $\widehat \rho_t$ defined analogously to $\kappa_i$ and $\rho_t$ in (ref), but using the fixed effect estimates instead of the true value of the parameters. All derivatives of $\widehat \ell^*_{it}$ are evaluated at $(\widehat \beta,\widehat \alpha_i - \widehat \kappa_i' \widehat \beta, \widehat \gamma_t - \widehat \rho_t' \widehat \beta)$ in the following. Then, the plug-in estimators of $B$ and $D$ are
and
where $\widehat H = - \mathbb{E}_{NT} \widehat \ell^{*\, \beta \beta}_{it}$, and $M$ is a trimming parameter such that $M/T \to 0$ and $M\to\infty$ as $T \to \infty$ \citep*{HK2007}. We can set $M=0$ when the covariates are strictly exogenous. These estimators use $\widehat \beta$, i.e. $\widehat B = \widehat B(\widehat \beta)$ and $\widehat D = \widehat D(\widehat \beta)$. It is possible to iterate the correction by (i) estimating $B$ and $D$ using $\widehat \beta_{\rm ABC}$ to obtain $\widehat B_2 = \widehat B(\widehat \beta_{\rm ABC})$, $\widehat D_2 = \widehat D(\widehat \beta_{\rm ABC})$ and $\widehat \beta_{\rm ABC2} = \widehat \beta - \widehat B_2/T - \widehat D_2/N$; (ii) estimating $B$ and $D$ using $\widehat \beta_{\rm ABC2}$ to obtain $\widehat B_3 = \widehat B(\widehat \beta_{\rm ABC2})$, $\widehat D_3 = \widehat D(\widehat \beta_{\rm ABC2})$ and $\widehat \beta_{\rm ABC3} = \widehat \beta - \widehat B_3/T - \widehat D_3/N$; and so on. The iteration does not affect the asymptotic distribution of the correction but might improve its small sample properties.
\citet*{HahnKuersteiner2011} developed the ABC for general dynamic nonlinear models with unobserved individual effects, building on the analysis of \citet*{Hahn:2002p717} for dynamic linear models and \citet*{Hahn:2004p882} for static nonlinear models. FernandezValWeidner2016 extended the ABC to models with unobserved individual and time effects.
This correction is based on the Jackknife method introduced by \citet*{quenouille1956notes} and \citet*{tukey58} for cross-sectional data. We start by giving some intuition on how this method works using the example of Section (ref).\footnote{We refer to ShaoTu1995 for a rigourous analysis of the properties of the Jackknife.} Let $\widehat \theta_{(-j)}$ denote the estimator of $\theta$ that leaves out the $j^{th}$ observation. The leave-one-out bias corrected estimator is $$ \widehat \theta_{\rm JC} = n \widehat \theta - (n-1) \bar{\theta}_{n-1}, \ \ \bar{\theta}_{n-1} = \frac{1}{n} \sum_{j=1}^n \widehat \theta_{(-j)}. $$ To understand how the correction works, assume the second-order expansion for the bias ${\mathrm{E}} \widehat \theta - \theta_0 = B_1/n + B_2/n^{2} + o(n^{-2})$. Then, under identical distribution ${\mathrm{E}} \widehat \theta_{(-j)} - \theta_0 = B_1/(n-1) + B_2/(n-1)^{2} + o(n^{-2})$ for all $j$, so that $$ {\mathrm{E}} \widehat \theta_{\rm JC} - \theta_0 = B_1 + \frac{B_2}{n} - B_1 - \frac{B_2}{n-1} + o(n^{-1}) = o(n^{-1}). $$ In other words $(n-1)(\bar{\theta}_{n-1} - \widehat \theta)$ is an estimator of the first-order bias since $ (n-1){\mathrm{E}}(\bar{\theta}_{n-1} - \widehat \theta) = B_1/n + o(n^{-1}). $
\citet*{Hahn:2004p882} introduced the Jackknife to panel models with individual effects, and FernandezValWeidner2016 extended it to panel models with individual and time effects. To describe how to apply the Jackknife to panel data, it is convenient to introduce some notation. Let ${\bf N} = \{ 1,\ldots, N\}$ and ${\bf T} = \{ 1,\ldots, T\}$ be the sets of indexes for the two dimensions of the panel. For the subsets of indexes $A \subseteq {\bf N}$ and $C \subseteq {\bf T}$, let $\widehat \beta_{A,C}$ denote the FE estimator of $\beta$ in the subpanel with indexes $(i,t) \in A \times C = \{(i,t): i \in A, t \in C\}$, that is
where $|A|$ denotes the cardinality of the set $A$. With this notation the FE estimator of $\beta$ is $\widehat \beta = \widehat \beta_{{\bf N},{\bf T}}$. Define the average leave-one-out estimators in each dimension as
The panel Jackknife estimator for models with individual effects is $T \widehat \beta_{{\bf N},{\bf T}} - (T-1)\overline \beta_{N,T-1}$. The correction $(T-1)(\overline \beta_{N,T-1} - \widehat \beta_{{\bf N},{\bf T}})$ removes the bias term $B/T$, analogously to the cross-sectional case discussed above. A panel jackknife bias corrected (JBC) estimator that removes both bias terms can be formed as $$ \widehat \beta_{\rm JBC} = (N+T-1) \widehat \beta_{{\bf N},{\bf T}} - (N-1) \overline \beta_{N-1,T} - (T-1) \overline \beta_{N,T-1}, $$ This correction removes the bias in large samples because $(N-1){\mathrm{E}}(\overline{\beta}_{N-1,T} - \widehat \beta_{{\bf N},{\bf T}} ) \approx D/N$ and $(T-1){\mathrm{E}}(\overline{\beta}_{N,T-1} - \widehat \beta_{{\bf N},{\bf T}} ) \approx B/T$. Moreover, it can be shown that under suitable conditions $$ \widehat \beta_{\rm JBC} - \beta_0 \overset{a}{\sim} {\cal N}\left( 0, \frac{H^{-1}} {NT} \right). $$ An important limitation in the application of the leave-one-out Jackknife to panel data is that it requires independence in both dimensions of all the variables. This is a very restrictive condition for panel data, specially when one of the dimensions is time, as it rules out lagged-dependent variables or serially correlated variables as covariates.
This correction is based on the split-sample Jackknife method introduced by \citet*{quenouille1949approximate} for time series data. We again start by giving an intuitive description of this method in the context of the simple model of Section (ref). Assume that the sample size $n$ is even and split the sample in two halves: $\{z_j : 1 \leq j \leq n/2\}$ and $\{z_j : n/2 + 1 \leq j \leq n\}$. Let $\widehat \theta^1$ and $\widehat \theta^2$ denote the estimator of $\theta$ in each of the half-samples. The split-sample jackknife bias corrected estimator is $$ \widehat \theta_{\rm SC} = 2 \widehat \theta - \widetilde{\theta}_{1/2} = \widehat \theta - (\widetilde{\theta}_{1/2} - \widehat \theta), \qquad \widetilde{\theta}_{1/2} = (\widehat \theta^1 + \widehat \theta^2)/2. $$ Assume the first-order expansion for the bias ${\mathrm{E}} \widehat \theta - \theta_0 = B/n + o(n^{-1})$. The correction works asymptotically because under stationarity ${\mathrm{E}} \widehat \theta^1 = {\mathrm{E}} \widehat \theta^2 = \theta_0 + 2B/n + o(n^{-1}) $, so that $$ {\mathrm{E}} \widehat \theta_{\rm SC} - \theta_0 = \frac{2B}{n} - \frac{2B}{n} + o(n^{-1}) = o(n^{-1}). $$ In other words $\widetilde{\theta}_{1/2} - \widehat \theta$ is an estimator of the bias since $ {\mathrm{E}}(\widetilde{\theta}_{1/2} - \widehat \theta) = B/n + o(n^{-1}). $
To define the panel version of this Jackknife correction we introduce
where $\lceil . \rceil$ and $\lfloor . \rfloor$ are the ceiling and floor functions, and we use the same notation as in the previous section. \citet*{DhaeneJochmans2015} introduced the split-sample jackknife to panel models with individual effects. In this case the split-sample bias corrected (SBC) estimator is $2 \widehat \beta_{{\bf N},{\bf T}} - \widetilde \beta_{N,T/2} $. The correction $\widetilde \beta_{N,T/2} - \widehat \beta_{{\bf N},{\bf T}}$ removes the bias term $B/T$, analogously to the cross-sectional case. FernandezValWeidner2016 extended the SBC estimator to panel models with individual and time effects, in which case the correct linear combination is $$ \widehat \beta_{\rm SBC} = 3 \widehat \beta_{{\bf N},{\bf T}} - \widetilde \beta_{N/2,T} - \widetilde \beta_{N,T/2}. $$ This correction removes both bias terms in large samples because ${\mathrm{E}}(\widetilde{\beta}_{N/2,T} - \widehat \beta_{{\bf N},{\bf T}} ) \approx D/N$ and ${\mathrm{E}}(\widetilde{\beta}_{N,T/2} - \widehat \beta_{{\bf N},{\bf T}} ) \approx B/T$. The intuition here is simple. The estimator $\widetilde{\beta}_{N,T/2}$ has double the bias than $\widehat \beta_{{\bf N},{\bf T}}$ coming from the estimation of the individual effects, because there are only half of observations, $T/2,$ informative about each of them. However, $\widetilde{\beta}_{N,T/2}$ has the same bias as $\widehat \beta_{{\bf N},{\bf T}}$ coming from the estimation of the time effects, because there are the same number of observations, $N$, informative about each of them. A similar argument shows that the bias of $ \widetilde \beta_{N/2,T} $ is $2D/N + B/T$. Under suitable conditions, $$ \widehat \beta_{\rm SBC} - \beta_0 \overset{a}{\sim} {\cal N}\left( 0, \frac{H^{-1}} {NT} \right). $$ \citet*{quenouille1956notes} already noted that Jackknife correction can be extended to higher order, and \citet*{DhaeneJochmans2015} discussed high-order Jackknife corrections, that is, elimination not only of the leading bias of order $T^{-1}$, but also of the next order bias term of order $T^{-2}$, which can again be achieved by an appropriate linear combination of split-panel fixed-effect estimates. A formal discussion of those higher-order bias terms is given in \citet*{BunKiviet2003} for dynamic linear models panel models and in \citet*{SunDhaene2017} for non-linear models.
When the data are independent across one of the dimensions of the panel, it is not clear how to carry out the sample split in this dimension. One possibility is to make multiple splits and average the resulting corrected estimators for each split. Another possibility is a hybrid method that uses leave-one-out along the independent dimension and split-sample along the weakly dependent dimension. For example, if we assume independence over the cross-sectional dimension, then a hybrid jacknife estimator is $$ \widehat {\beta}_{\rm HBC} = (N+1)\widehat{\beta}_{{\bf N},{\bf T}} - (N-1)\overline{\beta}_{N-1,T} - \widetilde{\beta}_{N,T/2}. $$ This correction removes the bias in large samples because $(N-1){\mathrm{E}}(\overline{\beta}_{N-1,T} - \widehat \beta_{{\bf N},{\bf T}} ) \approx D/N$ and ${\mathrm{E}}(\widetilde{\beta}_{N,T/2} - \widehat \beta_{{\bf N},{\bf T}} ) \approx B/T$. Again, under suitable conditions $$ \widehat \beta_{\rm HBC} - \beta_0 \overset{a}{\sim} {\cal N}\left( 0, \frac{H^{-1}} {NT} \right). $$ We refer to Stata2015 for other examples of hybrid corrections.
All the corrections produce estimators with the same asymptotic distribution, but they rely on different sampling conditions and might involve the choice of tuning parameters. Compared to the JBC; ABC, SBC and HBC work under more general sampling conditions because they do not require independence of the data along both dimensions. However, SBC and HBC usually increase dispersion in small samples because they estimate the parameters in smaller panels than JBC. Compared to ABC, SBC and HBC work under less general sampling conditions because they rely on homogeneity in both dimensions in order for the bias to be the same in all the subpanels. ABC does not use homogeneity, but requires coding the estimators of the bias, which involves a delicate choice for the trimming parameter $M$ when the covariates are not strictly exogenous. In practice, we recommend to always carry out a sensitivity analysis reporting results from both jackknife and analytical corrections and to try several values of $M$ in the analytical correction starting with $M=0,1,2,3,4$.
Some of the differences can be appreciated in a simple cross sectional example where we can characterize the moments of all the estimators. Consider the normal model $z \sim \mathcal{N}(\mu,\sigma^2)$, $\sigma^2 > 0$. The MLE of $\sigma^2$ is the sample variance $\widehat \sigma^2 = n^{-1} \sum_{j=1}^n (z_j - \bar z)^2,$ where $\bar z = n^{-1} \sum_{j=1}^n z_j$. As it is well-known, $\widehat \sigma^2$ is biased with ${\rm Bias}(\widehat \sigma^2) = - \sigma^2/n$ and $\mathrm{Var}(\widehat \sigma^2 ) = 2 \sigma^4 (n-1)/n^2$. The ABC is $\widetilde \sigma^2_{\rm ABC} = (n+1) \widehat \sigma^2/n$, which has ${\rm Bias}( \widetilde \sigma^2_{\rm ABC})= - \sigma^2/n^2$ and $\mathrm{Var}(\widetilde \sigma^2_{\rm ABC} ) = (n+1)^2\mathrm{Var}(\widehat \sigma^2 )/n^2$. Iterating $k$ times yields $\widetilde \sigma^2_{\rm ABCk} = (\sum_{r=0}^{k} n^{-r}) \widehat \sigma^2$, which has ${\rm Bias}(\widetilde \sigma^2_{\rm ABCk})= - \sigma^2/n^{k+1}$ and $\mathrm{Var}(\widetilde \sigma^2_{\rm ABCk} ) = (\sum_{r=0}^{k} n^{-r})^2 \mathrm{Var}(\widehat \sigma^2 )$. The JBC can be shown to be $\widetilde \sigma^2_{\rm JBC} = n \widehat \sigma^2/(n-1),$ the degrees-of-freedom corrected estimator, which has ${\rm Bias}( \widetilde \sigma^2_{\rm JBC})= 0$ and $\mathrm{Var}(\widetilde \sigma^2_{\rm JBC} ) = n^2 \mathrm{Var}(\widehat \sigma^2 )/(n-1)^2$. Let $\bar z_1$ and $\bar z_2$ be the sample means of $z$ in the first and second halves of the sample. Simple algebra shows that the SBC is $\widetilde \sigma^2_{\rm SBC} = \widehat \sigma^2 + \bar z^2 - \bar z_1 \bar z_2,$ which has ${\rm Bias}( \widetilde \sigma^2_{\rm SBC})= 0$ and $\mathrm{Var}(\widetilde \sigma^2_{\rm SBC} ) = (n+2) \mathrm{Var}(\widehat \sigma^2 )/n$. A comparison of biases and variances reveals that in this case ${\rm Bias}(\widehat \sigma^2) > {\rm Bias}( \widetilde \sigma^2_{\rm ABC}) > {\rm Bias}( \widetilde \sigma^2_{\rm ABCk}) > {\rm Bias}( \widetilde \sigma^2_{\rm JBC}) = {\rm Bias}( \widetilde \sigma^2_{\rm SBC})$, while $\mathrm{Var}(\widehat \sigma^2 ) < \mathrm{Var}(\widetilde \sigma^2_{\rm ABC} ) < \mathrm{Var}(\widetilde \sigma^2_{\rm ABCk} ) < \mathrm{Var}(\widetilde \sigma^2_{\rm JBC} ) < \mathrm{Var}(\widetilde \sigma^2_{\rm SBC} )$.
Table (ref) reports the results of a numerical simulation for a probit model with strictly exogenous covariates and unobserved individual and time effects. The design is calibrated to the female labor force participation (LFP) application in \citet*{FernandezVal:2009p3313}. Thus, we draw 500 panels of size $N=664$ and $T=9$, with $$ y_{it} = 1(\beta'x_{it} + \alpha_i + \gamma_t \geq \varepsilon_{it}), \ \ i = 1,\ldots, N, \ \ t = 1, \ldots, T, $$ where $x_{it}$ includes three fertility variables (the numbers of children aged 0-2, 3-5, and 6-17), the logarithm of the husband's earnings in 1995 thousands of dollars, and a quadratic function of age in years divided by 10, whose values are taken for the PSID 1980--1988; and $\varepsilon_{it}$ are independent draws from the standard normal distribution.\footnote{The original PSID sample includes $1,461$ women, but only $664$ of them have variation in the LFP variable over the years of the panel.} The parameters $(\beta, \alpha_i, \gamma_t)$ are calibrated to the FE probit estimates in the PSID 1980--1988 with the observed LFP as the dependent variable. The table reports biases, standard deviations, root mean square errors, and empirical coverage probabilities of confidence intervals with nominal level of 95% for the FE, ABC, SBC, HBC, and JBC estimators of the coefficients of the three fertility variables. We consider two versions of the SBC, one based on one partition over the cross sectional dimension following the ordering in the data set, and the other based on the average over 50 random partitions. All the results, except for the coverage probabilities, are in percentage of the true value of the parameter. We find that the ABC and JBC drastically reduce bias, dispersion and rmse, and have coverage probabilities close to their nominal level. The evidence for the SBC and HBC is more mixed. While always reducing bias, they increase dispersion resulting in higher rmse and lower coverage than the FE for the kids6_17 fertility variable. Increasing the number of partitions in the SBC does not improve the finite sample performance of this correction.
The methods discussed so far correct the estimator. The same methods can also be applied to correct the first order conditions of the FE estimator of $\beta$. Namely, let $L(\beta) = \max_{(\alpha, \gamma) \in \mathbb{R}^{N+T}} \, \mathbb{E}_{NT} \ell_{it}(\beta, \alpha_i,\gamma_t)$ be the profile objective function of (ref), where all the fixed effects have been concentrated out. Then, a second-order expansion similar to (ref) yields
where $b=HB$ and $d=HD$. This expression was derived in \citet*{BesterHansen2009} for models with individual effects and \citet*{Sun2016} and \citet*{JochmansOtsu2016} for models with individual and time effects. Thus, the score of $L(\beta)$ at the true parameter value is not centered at zero, which is the source of the bias in $\widehat \beta$. A “profile-score bias corrected estimator” is the solution to $$d L(\widehat \beta_{\rm PSBC}) / d \beta + \widehat b/T + \widehat d/N= 0,$$ where $\widehat b= \widehat b(\widehat \beta)$ and $\widehat d = \widehat d(\widehat \beta)$ are consistent estimators for $b$ and $d$, respectively. An alternative score correction that estimates simultaneously the parameter and bias terms can be formed as the solution to $$ d L(\widehat \beta_{\rm PSBC2}) / d \beta + \widehat b(\widehat \beta_{\rm PSBC2})/T + \widehat d(\widehat \beta_{\rm PSBC2})/N= 0. $$ The analytical properties of this corrected profile-score are model specific, since they depend on the functional form of $ \widehat b(\beta) $ and $ \widehat d(\beta) $. In particular, they may have multiple or no solution, even in the case of concave log-likelihood, where the solution for the previous equation for $\widehat \beta_{\rm PSBC}$ is unique. \citet*{DhaeneJochmans2016likelihood} discussed this issue in detail for linear autoregressive panel models, and dhaene2015profile explored the behaviour of these profile-score adjustments for other panel data models with individual effects. \citet*{DhaeneJochmans2015} developed jackknife methods to correct the scores in models with individual effects.
Finally, \citet*{gonccalves2015bootstrap} show that bootstrap methods can correct for the Nickell bias in dynamic linear panel models, and \citet*{KimSun2016} show how to use the bootstrap to construct bias-corrected estimators in non-linear panel models.
The objects of interest in panel models are often ceteris paribus or partial effects, i.e. effects in the outcome of changing each covariate while holding the rest of covariates and unobserved effects fixed \citep*{Chamberlain1984}. In linear models the parameters correspond to partial effects. For example, in the model (ref), the components of $\beta$ measure the partial effects of each of the covariates. This effect is the same for all the individuals and time periods by linearity. In nonlinear models the partial effects depend not only on the parameters but also on the covariates and unobserved effects. For example in the nonlinear model for the conditional expectation ${\mathrm{E}}[y_{it} \mid x_{it},\alpha_i,\gamma_t] = m(x_{it}, \alpha_i, \gamma_t, \beta)$, the partial effect of changing the covariates from $x_{it}^0$ to $x_{it}^1$ is $$ \delta_{it}(\alpha_{i},\gamma_{t}, \beta) = m(x^1_{it}, \alpha_{i}, \gamma_t,\beta) - m(x^0_{it}, \alpha_i, \gamma_t,\beta), $$ where $x^1_{it}$ and $x^0_{it}$ might depend on $x_{it}$. If $x_{it}$ is continuous and $x \mapsto m(x, \alpha_i, \gamma_t, \beta)$ a.s. differentiable, the partial effect of a marginal change in the covariates is $$ \delta_{it}(\alpha_i,\gamma_t, \beta) = \partial m(x, \alpha_i, \gamma_t)/\partial x \big|_{x = x_{it}}. $$ In both cases the partial effects are heterogenous across $i$ and $t$.
\paragraph*{Example (ii)} In the binary response model, if the $k$th covariate $x_{k,it}$ is binary, then the partial effect of changing $x_{k,it}$ from $0$ to $1$ on the probability of $y_{it}=1$ conditional on the rest of the covariates and unobserved effects is $$ \delta_{it}(\alpha_i,\gamma_t, \beta) = F_{\varepsilon}(\beta_k + x_{it,-k}'\beta_{-k} + \alpha_i + \gamma_t) - F_{\varepsilon}(x_{it,-k}'\beta_{-k} + \alpha_i + \gamma_t) , $$ 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.
\paragraph*{Example (iii)} In the count response example, if the $k$th covariate $x_{k,it}$ is continuous, then the partial effect of a marginal change in $x_{k,it}$ on the average of $y_{it}$ conditional on the rest of the covariates and unobserved effects is $$ \delta_{it}(\alpha_i,\gamma_t, \beta) = \beta_k \exp(x_{it}'\beta + \alpha_i + \gamma_t), $$ where $\beta_k$ is the $k$th element of $\beta$.
One way to summarize the heterogeneity in the partial effects is to average them across the individuals and time periods in the panel. This yields the in-sample average partial effect (APE) $$ \delta_{NT} = {\mathbb{E}_{NT}} \delta_{it}(\alpha_i,\gamma_t, \beta). $$ Another possibility is to assume that there is a population of individuals and time periods from where the observed panel is drawn and consider the APE in this population. Under weak conditions this in-population APE corresponds to $$ \delta = \operatorname*{plim}_{N,T\to\infty} \delta_{NT}. $$ Yet another possibility is to average the partial effects of the individuals in the population over the time periods in the sample $$ \delta_T = \operatorname*{plim}_{N\to\infty} \delta_{NT}. $$ The choice of the relevant APE is application-specific.
The fixed effect estimator is the same for all the previous APEs. Thus, applying the plug-in principle, $$ \widehat \delta_{NT} = {\mathbb{E}_{NT}} \delta_{it}(\widehat \alpha_i, \widehat \gamma_t, \widehat \beta). $$ However, the asymptotic distribution of $\widehat \delta_{NT} $ depends on the APE of interest. The asymptotic distribution around $\delta_{NT}$ is $$ \widehat \delta_{NT} - \delta_{NT} \overset{a}{\sim} \mathcal{N}\left( \frac{E}{T} + \frac{F}{N}, \frac{\Sigma}{NT}\right), $$ where the expressions for the bias terms $E$, $F$, and the variance $\Sigma$ are given in FernandezValWeidner2016. The bias and variance here come from the estimation of the parameters and unobserved effects. The distribution around $\delta$ becomes
where the additional variance term, $ \Omega = \operatorname*{plim}_{N,T \to \infty} \left({\mathbb{E}_{NT}} [\delta_{it}(\alpha_i,\gamma_t,\beta) - \delta]\right)^2/a_{NT}, $ comes from the estimation of the population mean $\delta$ using the sample mean $\delta_{NT}$, and $a_{NT} = N \vee T$ because of the correlation of $\delta_{it}(\alpha_i,\gamma_t,\beta)$ with $\delta_{jt}(\alpha_j,\gamma_t,\beta)$ and $\delta_{is}(\alpha_i,\gamma_s,\beta)$ induced by the individual and time effects. Similarly, the distribution around $\delta_T$ is $$ \widehat \delta_{NT} - \delta_{T} \overset{a}{\sim} \mathcal{N}\left( \frac{E}{T} + \frac{F}{N}, \frac{\Sigma + T \Omega_{T}}{NT} \right), $$ where $ \Omega_{T} = \operatorname*{plim}_{N \to \infty} \left({\mathbb{E}_{NT}} [\delta_{it} (\alpha_i,\gamma_t,\beta)- \delta_T]\right)^2/T. $ Analogously to the parameters, the bias in the asymptotic distribution of the estimators of the APEs can be removed using analytical and jackknife corrections \citep*{FernandezValWeidner2016}.
The rate of convergence of $\widehat \delta_{NT}$ differs depending on the APE of interest. Thus, $\delta_{NT}$ can be estimated at the rate $\sqrt{NT}$, the same rate as the parameter $\beta$, whereas $\delta$ can be estimated at the rate $\sqrt{NT/a_{NT}} = \sqrt{N} \wedge \sqrt{T}$ and $\delta_T$ at the rate $\sqrt{T}$. As a consequence of the slower rate of convergence, the asymptotic distribution (ref) simplifies to
after dropping terms of lower order.\footnote{Similarly the asymptotic distribution around $\delta_T$ simplifies to $\widehat \delta_{NT}-\delta_{T} \overset{a}{\sim} \mathcal{N}(0 , \Omega_{T}/N)$.} In other words the estimation of the parameters does not have first-order effect on the distribution around $\delta$. The reason is that the bias and standard deviation introduced by the parameter estimation are of lower order compared to $\sqrt{NT/a_{NT}} $. Despite this possible asymptotic simplification, we recommend using the higher-order distribution in (ref) to perform inference on $\delta$, because it provides a more accurate approximation to the finite-sample distribution of $\widehat \delta_{NT}$ than the first-order distribution (ref).
The FE program (ref) might look like computationally challenging due to the high dimension of the fixed effects. However, there are two aspects of the program that greatly facilitate the computation. First, the objective function is concave and smooth in the most commonly used cases such as the linear, Poisson, logit, probit, ordered probit, and tobit models. Second, the design matrix is sparse allowing the use of sparse algebra methods to speed up the computation. The SBC correction preserves the computational properties of the FE estimator as it only involves solving the program (ref) in a small number of subsamples. The JBC and HBC can be computationally more intensive than the SBC when the cross-sectional dimension is large. The ABC also preserves the computational properties of the FE estimator, but requires to code the estimators of the bias that are model specific.
The bias corrections are available in the statistical software packages \verb"Stata" and \verb"R" for some models. Stata2015 provided two \verb"Stata" commands that implement the analytical and several jackknife corrections for estimators of parameters and APEs in logit and probit models with individual and time effects. These commands also report standard errors for the estimates of the parameter and APEs constructed using the asymptotic distribution. \citet*{Sun2016stata} implemented in \verb"Stata" the split-sample jackknife estimator and score corrections for parameters in linear, logit and probit models with individual effects. The command also has functionality for other single index user-supplied models and reports standard errors based on the asymptotic distribution. All the previous \verb"Stata" commands provide functionality for unbalanced panels. \citet*{StammannHeissMacFadden2016} implemented the analytical correction in \verb"R" for estimators of parameters and APEs in logit and probit models with exogenous covariates and individual effects.
In the previous section, we assumed that the observations for all the combinations of the two indexes $i$ and $t$ were observed, i.e. the panel was balanced. In empirical applications, however, it is common to have unbalanced panels where some of the observations are missing due to sample attrition. This does not introduce special theoretical complications provided that the source of the missing observations is random. Perhaps due to this reason, we are not aware of any work in the panel literature that extends the bias corrections explicitly to unbalanced panels. To fill this void, we provide without proof the asymptotic distribution of the FE for unbalanced panels.
Let ${\cal D}$ be the set of all observed pairs $(i,t)$, and $n = |{\cal D}|$ be the sample size. Define the operators $\mathbb{E}_n := n^{-1} \sum_{(i,t) \in {\cal D}}$, $ \mathbb{E}_{T,i} = |{\cal D}_i|^{-1} \sum_{t \in {\cal D}_i}$, and $ \mathbb{E}_{N,t} = |{\cal D}_t|^{-1} \sum_{i \in {\cal D}_t}$, where ${\cal D}_i = \{ t : (i,t) \in {\cal D}\}$ and ${\cal D}_t = \{ i : (i,t) \in {\cal D}\}$. Let $a_{it}$ be an attrition indicator for the observation $(i,t)$. The key regularity conditions that are needed to derive the bias expressions are (i) $|{\cal D}_i|/T \to c_i >0$ as $T \to \infty$ for all $i$; (ii) $|{\cal D}_t|/N \to c_t >0$ as $N \to \infty$ for all $t$; and (iii) $y_{it}$ is independent of $a_{it}$ conditional on $(x^t_i, \alpha, \gamma)$. The first two conditions guarantee that the number of observations for each unobserved effects increases with the sample size. The third condition imposes conditional missing at random in the attrition process. The asymptotic distribution of $\widehat \beta$ then becomes
where $H_u = - \operatorname*{plim}_{n \rightarrow \infty} \mathbb{E}_n \overline \ell^{*\, \beta \beta}_{it}$, $\overline T = n / N$ is the average number of observations available for each cross-sectional unit,
$\overline N = n / T$ is the average number of observations available for each time period, and
This result agrees with formula (ref), because the order of the bias is $\overline T^{-1} B + \overline N^{-1} D \sim (N+T)/n$. The magnitude of the bias depends on the averages $\overline T$ and $\overline N$, so having some units $i$ and $t$ with a very small number of observations is not necessarily a serious problem from the perspective of incidental parameter bias.
Regarding bias correction, the ABC for the case of unbalanced panels can be constructed using the empirical analog of the bias expressions evaluated at the FE estimator. Jackknife corrections can be formed by partitioning the panel in the same way as in the balanced case, i.e. without taking into account the attrition. We refer to DhaeneJochmans2015, Stata2015, ChudikPesaranYang2016 and Sun2016stata for details.
We now briefly discuss the case where the unobserved effects are vector-valued, i.e. $\alpha_i \in \mathbb{R}^{d_{\alpha}}$ and $\gamma_t \in \mathbb{R}^{d_{\gamma}}$ with $d_{\alpha} \geq 1$ and $d_{\gamma} \geq 1$. We start by providing some motivating examples.
\paragraph*{Example (iv)} An important case is the normal linear model with interactive effects or factor structure, $y_{it} = x_{it}'\delta + \alpha_{i}'\gamma_{t} + \varepsilon_{it},$ $\varepsilon_{it} \mid x^t_i, \alpha, \gamma, \beta \sim \mathcal{N}(0,\sigma^2),$ where
Here, $d_{\alpha} = d_{\gamma}$ is the number of interactive effects. \citet*{bai2009panel} showed that this model does not suffer from the incidental parameter problem when the covariates are strictly exogenous and the semiparametric model is correctly specified. However, he found that the FE least squares estimator has an asymptotic bias structure of the form $B/T+D/N$ under conditional heteroskedasticity and derived the expressions for $B$ and $D$ that can be used for analytic bias correction. \citet*{MoonWeidner2017} showed that the FE least squares estimator of $\delta$ suffers from a \citet*{Nickell81} type incidental parameter bias when the covariates are predetermined, even when the model is correctly specified. We refer to \citet*{bai2016econometric} for a recent review, and to \citet*{Pesaran2006} and \citet*{ChudikPesaran2015} for alternative estimation methods of factor models.
\paragraph*{Example (v)} Another example is a nonlinear single index factor model where
where $f$ is a known function such as $f(y,u) = F_{\varepsilon}(u)^y(1-F_{\varepsilon}(u))^{1-y}1(y \in \{0,1\})$ with a CDF $F_{\varepsilon}$, for a binary response model, or $f(y,u) = [\exp(u)^y \exp(-\exp u)/y!] 1(y \in \{0,1, ... \})$ for a count response model. \citet*{CFW2014} characterized the bias of the FE estimator for models where the function $u \mapsto f(\cdot, u)$ is log-concave including the probit, logit, ordered probit and Poisson. See also \citet*{Chen2016estimation}, \citet*{BonevaLinton2017} and \citet*{Wang2017} for other articles discussing this model.
The analysis of Section (ref) carries over to the multivariate case with some minor adjustments. For the second-order expansion of $\widehat \beta$ it is still convenient to define $\ell^*_{it}(\beta,\alpha_i,\gamma_t)$ as in (ref), but now $ \kappa_i$ is a $d_{\beta} \times d_{\alpha}$ matrix and $\rho_t$ is a $d_{\beta} \times d_{\gamma}$ matrix, which are solutions to
Here we need to be more careful with the column and row dimensions of matrices in the notation. For example, $\overline \ell_{it}^{\, \beta \alpha_i'} := \partial(\partial \overline \ell_{it}/ \partial \beta)/ \partial \alpha_i'$ is the $d_\beta \times d_{\alpha}$ matrix of second derivatives of $\overline \ell_{it}$. A second-order asymptotic expansion similar to Section (ref) gives
where $H = - \operatorname*{plim}_{N,T \rightarrow \infty} \mathbb{E}_{NT} \overline \ell^{*\, \beta \beta'}$, $\Omega = \operatorname*{plim}_{N,T \rightarrow \infty} \mathbb{E}_{NT} {\mathrm{E}}\left( \overline \ell^{*\, \beta} \overline \ell^{*\, \beta'} \right)$, $B= H^{-1} b$ and $D = H^{-1} d$, and $b$ and $d$ are $d_\beta$-vectors with elements
and
In contrast to Section (ref), the expressions of the asymptotic bias and variance here do not use the information equality to simplify terms, e.g., $H^{-1} \Omega H^{-1}= H^{-1}$. These formulas remain valid in conditional moment models where only a characteristic of the conditional distribution such as the expectation is correctly specified. This covers for example the linear model with interactive effects under heteroskedasticity, that is, a version of Example (iv) that only imposes $$ {\mathrm{E}}[y_{it} \mid x_{it},\alpha_i, \gamma_t; \beta] = x_{it}'\beta + \alpha_{i}'\gamma_{t}. $$ The above formulas for asymptotic bias and variance thus contain those in \citet*{bai2009panel} and \citet*{MoonWeidner2017} as special cases.\footnote{See \citet*{Sun2016} for a related derivation.} However, these results rely on an auxiliary consistency proof, which is usually model specific and therefore can be delicate. For example, it is difficult to show consistency in models with interactive effects because the log-likelihood $(\beta,\alpha,\gamma) \mapsto {\mathbb{E}_{NT}} \ell_{it}(\beta,\alpha_i,\gamma_t)$ is not concave.
Analytical and jackknife corrections can be formed analogously to Section (ref). For details we refer to \citet*{ArellanoHahn2016} for models with individual effects, \citet*{bai2009panel} and \citet*{MoonWeidner2017} for linear models with interactive fixed effects, and to \citet*{CFW2014} for single-index factor models. \citet*{twostep11} derived bias corrections for two-step estimators of selection and other control variable models with unobserved individual effects in both steps. \citet*{chudik2015common} and \citet*{Westerlund2018} developed bias corrections for the common correlated effects estimator of \citet*{Pesaran2006} for linear factor models.
\citet*{CFW17} used nonlinear panel data methods to estimate distributional and quantile effects. The idea is to model the distribution of the outcome conditional on the covariates and unobserved effects via distribution regression, that is $$ \Pr(y_{it} \leq y \mid x_{it}, \alpha_i, \gamma_t) = F_{y}(x_{it}'\delta(y) + \alpha_i'\pi(y) + \gamma_t' \xi(y) ), \ \ \beta(y) = (\delta(y),\pi(y),\xi(y)), $$ where $F_{y}$ is a CDF such as the normal or logistic, $y \mapsto \beta(y)$ is a function-valued parameter, and the dimension of $\alpha_i$ and $\gamma_t$ is unrestricted. Both $F_y$ and $\beta(y)$ can vary with $y$ to accommodate heterogeneity accross the distribution. At each $y$, the parameters are estimated by a binary response regression of the outcome $1(y_{it} \leq y)$ on the covariates and fixed effects with the parametrization $\alpha_i(y) = \alpha_i'\pi(y)$ and $\gamma_t(y) = \gamma_t' \xi(y)$. The distributional and quantile effects are functionals of APEs. For example, the marginal distribution of the potential outcome corresponding to $x_{it}= x^0_{it}$ is the APE $$ F^0(y) = (NT)^{-1} \sum_{i=1}^N\sum_{t=1}^T F_{y}(\delta(y)'x^0_{it} + \alpha_i(y) + \gamma_t(y) ), $$ and the corresponding $\tau$-quantile is the functional $$ Q^0(\tau) = \inf\{y \in \mathbb{R} : F^0(y) \geq \tau\}. $$ The $\tau$-quantile effect of a change from $x_{it}= x^0_{it}$ to $x_{it}= x^1_{it}$ is
where $Q^1(\tau)$ is constructed using the same procedure as $Q^0(\tau)$ but with $x_{it}^1$ in place of $x_{it}^0$. \citet*{CFW17} applied bias corrections to estimate and make inference on quantile effects uniformly over quantile indexes. These methods applied to continuous, discrete and mixed continuous-discrete outcomes.
The effect in (ref) is a marginal quantile effect that has causal interpretation under standard conditional ignorability conditions. Alternative conditional effects can be estimated by modelling the distribution of the outcome conditional on the covariates and unobserved effects via quantile regression. The FE quantile regression estimator, introduced by \citet*{Koenker2004}, in general suffers from the incidental parameter problem. \citet*{KatoGalvaoMontes-Rojas2012} showed asymptotic normality and unbiasedness of the FE quantile regression estimator in models with only individual effects under sequences where $N/T\rightarrow 0$ as $N,T \to \infty$. \citet*{GalvaoKato2016smooth} derived the asymptotic distribution of the estimator in the same models, including the leading order $1/T$ bias, for a smoothed version of the quantile regression objective function under sequences for which $N$ and $T$ grow at the same rate. \citet*{ArellanoWeidner2017} showed that the bias of the FE estimator without smoothing is the same as in \citet*{GalvaoKato2016smooth} in models with strictly exogenous regressors. Thus, in terms of leading order bias and variance the FE quantile regression estimator behaves analogous to any other non-linear panel data model. However, the analysis in Section (ref) does not directly carry over to this case, because the quantile regression objective function is not sufficiently smooth. We refer to \citet*{GalvaoKato2018} for a recent review of quantile regression methods for panel data.
Another natural extension of the panel semiparametric model of Section (ref) is a model with three-way effects for data with a multidimensional structure. Let $\ell_{ijt}( \beta, \alpha_i , \gamma_j, \delta_t) = \log f(y_{ijt} \mid x_{ijt}, \alpha_i , \gamma_j, \delta_t; \beta)$ be the log-likelihood of an outcome $y_{ijt}$ conditional on $x_{ijt}$, with common parameter $\beta$ and unobserved effects $\alpha_i$, $\gamma_j$ and $\delta_t$. Here, we have three panel indices $i=1,\ldots,I$, $j=1,\ldots,J$ and $t=1,\ldots,T$. For example in an international trade application, \citet*{helpman2008estimating} are interested in analyzing the volume of trade from country $i$ to country $j$ at year $t$. The corresponding fixed effect estimator is
\paragraph*{Example (vi)} As an extension of Example (ii) above, consider a binary response single index model $$y_{ijt} = 1(x_{ijt}'\beta + \alpha_i + \gamma_j + \delta_t \geq \varepsilon_{it}), \quad \varepsilon_{it} \mid x^t_i, \alpha, \gamma, \delta \sim F_{\varepsilon},$$ where $F_{\varepsilon}$ is a known CDF. For example, $y_{ijt}$ can be an indicator of trade from country $i$ to country $j$ at year $t$, $x_{ijt}$ the level of tariffs from country $i$ to country $j$ at year $t$, $\alpha_i$ is an exporter country effect, $\gamma_j$ is an importer country effect, and $\delta_t$ is a year effect.
Similar to the problem of unbalanced panels above, we are not aware of any formal discussion of three-way fixed effect models in the econometrics literature on large panels, but such estimation problems certainly occur in applications. We provide a brief discussion here. In particular, we highlight the applicability of the formula (ref) to this case. According to that formula, the expected order of the incidental parameter bias in $\widehat \beta$ is $(I+J+T)/ (IJT)$, that is, we expect a leading order bias structure $B_1/(IJ) + B_2/(IT) + B_3/(JT) $. Compared to the order of the standard deviation of $\widehat \beta$, $1/\sqrt{IJT}$, the bias is small, as long as all three panel dimensions are sufficiently large, in which case the incidental parameter problem can be ignored. Thus, if $I,J,T \rightarrow \infty$ such that $I/JT \to 0$, $J/IT \to 0$ and $T/IJ \to 0$, then $\widehat \beta$ is asymptotically normal and unbiased, and standard MLE inference results are applicable to $\widehat \beta$, without requiring any bias correction.
So far we have discussed a relatively benign three-way fixed effect model. A more difficult case is $\ell_{ijt}( \beta, \alpha_{ij} , \gamma_{it}, \delta_{jt}) = \log f(y_{ijt} \mid x_{ijt}, \alpha_{ij} , \gamma_{it}, \delta_{jt}; \beta)$ and
Here, the unobserved effects $\alpha_{ij}$, $\gamma_{it}$ and $\delta_{jt}$ are pair-specific, only constant along one of the panel dimensions. It is still possible to estimate them consistently if $I,J,T \rightarrow \infty$. Applying the formula (ref), we expect the leading order incidental parameter bias to be of the form $B_1/I + B_2/J + B_3/T$ since $p=IJ+IT+JT$ and $n=IJT$. The order of this bias is never negligible compared to the order of the standard deviation of $\widehat \beta$, $1/\sqrt{IJT}$. For example, if $I = J = T$, then the order of the bias is $I^{-1}$, greater than the order of the standard deviation $I^{-3/2}$. It is therefore particularly useful to develop bias corrections for models with fixed effects in multiple dimensions.
We have reviewed developments for large panel data models since \citet*{ArellanoHahn2007}. A key insight is that the order of the incidental parameter bias of the fixed effects estimators can be obtained from the simple formula (ref). This formula is useful to get an initial idea of the relevance of the bias in a given application, but the exact magnitude of the bias is still very much model and data generating process dependent. In particular, it is impossible to give any general answer to the question of how large $T$ needs to be for the large-$T$ methods to perform well, see for example \citet*{Greene2004} for simulation results on various non-linear panel data models. We have discussed bias correction methods that estimate the actual size of the leading order bias. If this leading order bias is small (i.e. if the bias corrected estimator is close to the uncorrected estimator), then it is plausible to assume that the higher order bias terms are also small and can be neglected. If the leading order bias is large (i.e. if the bias corrected estimator is very different from the uncorrected estimator), then it is useful to conduct a Monte Carlo simulation calibrated to the application of interest, in order to verify that the remaining bias is small and that the inference procedure works well.
Fixed effect estimation methods for long panels have recently been applied to other data structures such as network and trade data. Examples include \citet*{Harrigan1996}, \citet*{AndersonWincoop2003}, \citet*{SantosSilvaTenreyro2006} and \citet*{helpman2008estimating}, which estimate gravity equations with unobserved importer and exporter country effects; and \citet*{Graham2017}, \citet*{Dzemski2017}, \citet*{Jochmans2017semiparametric}, \citet*{jochmans2017two}, \citet*{shi2016structural} and \citet*{Candelaria2016}, which apply large-$T$ panel methods to network data with unobserved sender and receiver effects. The large-$T$ panel methods are particularly well-suited for this type of data because typically $N=T$. Thus, trade/network usually correspond to square panel data sets where the two dimensions index the same set of countries/individuals as senders and receivers.
There are alternative semiparametric methods to FE that are not reviewed here. An interesting recent approach is the grouped fixed effects (GFE) of \citet*{hahn2010panel} and \citet*{bonhomme2015grouped} for models with individual effects. Compared to FE, GFE is less affected by the incidental parameter problem, because it restricts the distribution of the individual effects to be discrete. Large $T$ is still required to consistently estimate the group membership of individuals, but if the true heterogeneity structure is indeed discrete, then no asymptotic bias correction is required for the parameters of interest, even if $N$ is much larger than $T$. The GFE approach imposes restrictions on the distribution of the individual effects, but is more flexible in other dimensions than FE. For example, it permits to include individual effects that change over time. \citet*{bonhomme2017discretizing} characterized the approximation properties of GFE when the distribution of the individual effects is continuous. They showed that GFE is an effective dimension reduction device in this case, but suffers from the incidental parameter problem.