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.
69,955 characters · 8 sections · 102 citation commands
Fixed Effects Binary Choice Models: Estimation and Inference with Long Panels
\onehalfspacing
Empirical analyses explaining binary outcomes, such as labour force participation or exporting decisions, are quite common in economics. The increasing number and availability of large and long panel data sets offers several advantages to researchers compared to pure cross-sections or time series (see chapter 1.2 in b2013 and h2014 for a comprehensive list of advantages). Maybe the most important advantage is that they allow to control for different sources of unobserved heterogeneity. In panels it is natural to account for unobserved individual and time specific effects simultaneously, so-called two-way fixed effect models. The corresponding estimators treat the unobserved effects as additional parameters to be estimated and thus allow for unrestricted correlation patterns between the explanatory variables and the unobserved effects. As the researcher does not have to make any distributional assumptions about the unobserved heterogeneity, these models are very flexible and a natural candidate for many empirical applications.
In the early stage of panel data econometrics, panels consisted of relatively few observations per individual. Consequently, when deriving asymptotic properties of estimators, it is very often assumed that the number of individuals ($N$) grows and the number of points in time ($T$) is held fixed. Under this asymptotic framework, non-linear fixed effects estimators are inconsistent, known as the incidental parameter problem (IPP) first mentioned by ns1948. This strand of literature is therefore particularly interested in deriving fixed $T$ consistent estimators. For instance, so-called conditional logit estimators have been proposed for static and dynamic binary choice models with individual fixed effects (see r1960, a1970, c1980, and hk2000). However, it is not possible to derive fixed $T$ consistent fixed effects estimators for all kind of models, e. g. the probit model. Another drawback of all conditional logit estimators is that they preclude the estimation of partial effects, which are of great interest in economics (see ah2007, and fw2018a).
For these reasons, among others, and further motivated by the seminal work of pm1999 and the rising availability of comprehensive panel data, a growing literature now focuses on large $N$ and $T$ asymptotics. The beauty of this asymptotic framework is that the inconsistency problem of the IPP can be turned into an asymptotic bias problem that can be corrected. More precisely, hk2002 were the first to exploit this asymptotic framework and to propose a bias-corrected estimator for dynamic linear panel models with individual fixed effects to address the inference problem induced by the n1981 bias. hm2006 show that the same bias correction is also applicable if the model includes additional time fixed effects. In the meantime, several bias-corrected estimators for non-linear models have been proposed (see among others l2002, w2002, hn2004, c2007, f2009, bh2009, dj2015, fw2016, and ks2016). A remarkable difference compared to estimators for linear models is that the inclusion of time fixed effects leads to an additional bias as shown by fw2016. We refer the interested reader to ah2007 and fw2018a for comprehensive overviews.
Another apparent challenge that discourages researchers from using non-linear fixed effects models is the computational burden associated with the estimation. This problem is especially severe when the model specification leads to high-dimensional fixed effects, which is already the case when $N$ or $T$ are large. If only one of the panel dimensions is large, the algorithm of g2004 can significantly reduce the computational burden. If both dimensions are large algorithms like gp2010 and s2018 can be used. From a practical point of view, however, it is not so obvious how these algorithms can be combined with analytical bias corrections like the one of fw2016.
In this article, we offer new insights that facilitate and validate the usage of different (bias-corrected) estimators for binary choice models with two-way fixed effects in empirical research. First, we show how to address the computational obstacles that often prevent the application of bias corrections. Second, we extend the simulation experiments of fw2016 by several aspects to gain deeper insights into the statistical properties of different estimators. More specifically, we analyse further analytical and split-panel jackknife bias-corrected estimators that were proposed but not studied by the authors. We also consider alternative estimators for average partial effects based on linear fixed effects models. These models are often used in empirical research to avoid the above mentioned pitfalls of non-linear ones. Because many real world panel data sets are unbalanced, our analysis also considers different patterns of randomly missing data. This aspect has so far received little attention in the literature, but is relevant for many empirical applications. Overall, we find that analytical bias corrections are preferable to split-panel jackknife approaches. In general, the latter show higher distortion, lower coverage, and are less robust to different patterns of randomly missing data. Third, we provide an illustrative example using an unbalanced panel from the German Socio-Economic Panel (see wfs2007) to investigate the inter-temporal labour force participation of 6,241 women between 1984 and 2013. Inspired by h1999 we estimate a dynamic fixed effects probit model where bias corrections are required to deal with the inference problem. Finally, we offer the analytical bias correction of fw2016 in our R package alpaca to encourage its application.\footnote{Until now, the analytical bias correction proposed by fw2016 is only provided in a Stata routine by cfw2017, which is not designed for long panels.}
The remainder of this article is organized as follows. In Section (ref) we introduce the model, various bias corrections, and algorithms for handling panels with large $N$ and $T$. In Section (ref) we provide results of extensive simulation experiments. In Section (ref) we apply different bias-corrected estimators to an empirical example from labour economics. Finally, we give some concluding remarks in Section (ref).
Throughout this article, we follow conventional notation: scalars are represented in standard type, vectors and matrices in boldface, and all vectors are column vectors.
The fixed effects binary choice model examined in this article can be derived from a latent variable model with two additive unobserved effects. Let
be the latent variable, where $i$ and $t$ are individual and time specific indexes.\footnote{Without loss of generality, $i$ and $t$ could also be indexes of networked economic activities like trade between countries.} However, instead of the latent variable, we only observe $y_{it} = 1$ if $y_{it}^{\ast} \geq 0$ and $y_{it} = 0$ otherwise. To allow for missing data, we define the following sets: $\mathcal{S}$ is a subset of $\{(i, t) \mid i \in \{1, \ldots, N\} \wedge t \in \{1, \ldots, T\}\}$ containing all observed pairs of indexes and $\mathcal{S}_{t} = \{i \mid (i, t) \in \mathcal{S}\}$ and $\mathcal{S}_{i} = \{t \mid (i, t) \in \mathcal{S}\}$ are subsets of $\mathcal{S}$ containing all indexes of individuals observed for a period $t$ and points in time observed for an individual $i$. $N$ and $T$ are the number of individuals and points in time and $n = \lvert \mathcal{S} \rvert$ is the sample size. Furthermore, $\mathbf{x}_{it}$ is the $it$-th row of the regressor matrix $\mathbf{X}$ and a $J$-dimensional vector of possibly predetermined explanatory variables, $\boldsymbol{\beta}$ are the corresponding structural parameters, $\alpha_{i}$ and $\gamma_{t}$ are incidental parameters that capture unobserved individual and time specific effects, and $e_{it}$ is an idiosyncratic error term assumed to be mean zero and independent of $\mathbf{X}_{i}^{\star} = (\mathbf{x}_{is})_{s \leq t \in \mathcal{S}_{i}}$ and $\boldsymbol{\pi} = (\boldsymbol{\alpha}, \boldsymbol{\gamma})$, where $\boldsymbol{\alpha} = (\alpha_{1}, \ldots, \alpha_{N})$ and $\boldsymbol{\gamma} = (\gamma_{1}, \ldots, \gamma_{T})$. Note that the imposed exogeneity assumption is milder than the often assumed strict exogeneity condition, where $e_{it}$ is independent of $\mathbf{X}_{i} = (\mathbf{x}_{is})_{s \in \mathcal{S}_{i}}$ instead of $\mathbf{X}_{i}^{\star}$. This stricter assumption is often too restrictive in the context of panel data, because it is violated if future realizations of a regressor are affected by the outcome variable, which for instance also precludes lagged dependent variables as regressors. In the presence of missing data, we have to additionally assume that conditional on $\mathbf{X}_{i}^{\star}$ and $\boldsymbol{\pi}$ the observations are missing at random.
Assuming a certain distribution for $e_{it}$ allows us to use the principle of maximum likelihood to derive a parametric estimator for fixed effects binary choice models. Let
be the log-likelihood contribution of individual $i$ at time $t$, where $F_{it}$ is the cumulative distribution function of $e_{it}$ evaluated at the linear index $\eta_{it} = \mathbf{x}_{it}^{\prime} \boldsymbol{\beta} + \alpha_{i} + \gamma_{t}$. Common choices for $F_{it}$ are the standard normal, the logistic, and the complementary log-log distribution. The corresponding maximum likelihood estimator is
Although fw2016 show consistency of $\hat{\boldsymbol{\beta}}$ under asymptotics where $N$ and $T$ grow at the same rate ($\operatorname{\text{plim}\,}_{N, T \rightarrow \infty} \hat{\boldsymbol{\beta}} = \boldsymbol{\beta}$), they also expose an asymptotic bias in the limiting distribution of the estimator with some severe consequences for inference. To get a better understanding of this specific inference problem, we briefly summarize the key findings of fw2016 for binary choice models and combine them with their conjecture about unbalanced panels stated in fw2018a.
Under asymptotic sequences where $N, T \rightarrow \infty$, $N / T \rightarrow \kappa^{2}$, and $0 < \kappa < \infty$ plus certain regularity conditions, like additively separable unobserved effects and concavity of the objective function, an asymptotic approximation to the limiting distribution of $\hat{\boldsymbol{\beta}}$ is given by
where $\overline{N} = n / T$ and $\overline{T} = n / N$ are the average number of individuals and points in time, $\mathbf{B}^{\beta}$ and $\mathbf{C}^{\beta}$ are the leading terms of the asymptotic bias $\mathbf{b}^{\beta} = \overline{T}^{- 1} \mathbf{B}^{\beta} + \overline{N}^{- 1} \mathbf{C}^{\beta}$ stemming from the inclusion of individual and time specific fixed effects, and $\mathbf{V}^{\beta}$ is the asymptotic covariance matrix. Due the asymptotic bias, the approximated limiting distribution is not correctly centred at $\boldsymbol{\beta}$, which means that, even if $n$ is large, confidence intervals constructed around any $\hat{\boldsymbol{\beta}}$ might not cover the true value of the corresponding parameter with probability close to the desired nominal level. However, this inference problem can be corrected by forming suitable estimators for $\mathbf{b}^{\beta}$ that can be subtracted from $\hat{\boldsymbol{\beta}}$.
Most researchers are not directly interested in the structural parameters, but rather in average partial effects. Let
denote the partial effect of a change in $x_{itj}$, where $x_{itj}$ is the $j$-th element in $\mathbf{x}_{it}$, $\partial_{\eta} F_{it}$ is the first-order partial derivative of $F_{it}$ with respect to $\eta_{it}$, and $F_{it}\rvert_{x_{itj} = k}$ indicates that $x_{itj}$ in the linear index is replaced by $k$.\footnote{For simplicity, we ignore the possibility of more complicated functional forms in the linear index, e. g. polynomials.} The average partial effects are then given by $\boldsymbol{\delta} = (\delta_{1}, \ldots, \delta_{J})$, where $\delta_{j} = n^{- 1} \sum_{(i, t) \in \mathcal{S}} \boldsymbol{\Delta}_{itj}$. Under some additional sampling and moment conditions, an asymptotic approximation to the limiting distribution of the estimator of the average partial effects is given by
where $\mathbf{V}^{\delta}$ is the asymptotic covariance matrix. Again, $\mathbf{B}^{\delta}$ and $\mathbf{C}^{\delta}$ are the leading terms of the asymptotic bias $\mathbf{b}^{\delta} = \overline{T}^{- 1} \mathbf{B}^{\delta} + \overline{N}^{- 1} \mathbf{C}^{\delta}$ stemming from the inclusion of individual and time specific fixed effects. Thus, as for $\hat{\boldsymbol{\beta}}$, there is an asymptotic bias problem.
In the next subsection we summarize the different bias corrections proposed by fw2016, but using a slightly modified notation. We change their notation for two reasons: First, we are considering panels that are potentially unbalanced, and second, we want to emphasize the link to recent advances in computational econometrics, which allow us to estimate (bias-corrected) binary choice models even when both panel dimensions are large.
Before we present the different bias corrections proposed by fw2016, we have to introduce some additional notation. Let $\partial_{z^{p}} G_{it}$ denote the $p$-th order partial derivative of an arbitrary function $G_{it}$ with respect to $z_{it}$. Further, let $\partial_{\eta} l_{it} = H_{it} (y_{it} - F_{it})$, $\omega_{it} = H_{it} \partial_{\eta} F_{it}$, $H_{it} = \partial_{\eta} F_{it} / (F_{it} (1 - F_{it}))$, and $\nu_{it} = \partial_{\eta} l_{it} / \omega_{it}$. We use vector notation to indicate that we collect the different quantities for all observations, e. g. $\boldsymbol{\omega} = (\omega_{it})_{(i, t) \in \mathcal{S}}$. Finally we define the residual projection $\operatorname{\mathbb{M}} = \operatorname{\mathbf{1}}_{n} - \operatorname{\mathbb{P}} = \operatorname{\mathbf{1}}_{n} - \mathbf{D}(\mathbf{D}^{\prime} \boldsymbol{\Omega} \mathbf{D})^{+} \mathbf{D}^{\prime} \boldsymbol{\Omega}$, where $\operatorname{\mathbf{1}}_{n}$ is an identity matrix of dimension $(n \times n)$, $\mathbf{D}$ is a sparse indicator matrix of dimension $(n \times N + T)$ arising from dummy encoding of individual and time identifiers, $(\cdot)^{+}$ refers to the Moore-Penrose inverse, and $\boldsymbol{\Omega}$ is a positive definite diagonal weighting matrix with $\operatorname{\text{diag}}(\boldsymbol{\Omega}) = \boldsymbol{\omega}$. The corresponding sample analogues are indicated by a hat. For clarification, we refer to $\hat{\eta}_{it} = \mathbf{x}_{it}^{\prime} \hat{\boldsymbol{\beta}} + \hat{\alpha}_{i} + \hat{\gamma}_{t}$ as the sample analogue of $\eta_{it}$. Table (ref) contains explicit expressions for distributions and the corresponding derivatives of frequently used binary choice models.
fw2016 distinguish between two types of bias corrections: analytical and split-panel jackknife. The latter exploits the relation between sample size and bias to form a non-parametric estimator of the asymptotic bias and is an extension of dj2015, whereas the former relies on explicit expressions derived from asymptotic expansions.\footnote{The idea to reduce bias using jackknife techniques originates from q1949, q1956. In the context of dynamic models they were first mentioned by h2002.} A bias-corrected estimator for $\boldsymbol{\beta}$ is
where $\hat{\mathbf{b}}^{\beta}$ is an estimator of the asymptotic bias such that $\tilde{\boldsymbol{\beta}} \operatorname{\overset{a}{\sim}\,} \operatorname{\mathcal{N}} (\boldsymbol{\beta}, \mathbf{V}^{\beta})$.
We start with an analytical bias correction at the level of the estimator. An explicit expression for an estimator of the asymptotic bias is
where
$L$ is the bandwidth parameter of the truncated spectral density estimator suggested by hk2007, and $\tau_{i}(l) = \lvert \mathcal{S}_{i} \rvert / (\lvert \mathcal{S}_{i} \rvert - l)$ is a finite sample adjustment proposed by fw2016. The corresponding estimator of the asymptotic covariance is $\widehat{\mathbf{V}}^{\beta} = \widehat{\mathbf{W}}^{- 1}$. By making the stronger strict exogeneity assumption, we can set $L = 0$ and drop the second term in $\widehat{\mathbf{B}}^{\beta}$, so that the expressions for $\widehat{\mathbf{B}}^{\beta}$ and $\widehat{\mathbf{C}}^{\beta}$ become identical except for the indexes.\footnote{The second term in $\widehat{\mathbf{B}}^{\beta}$ can be interpreted as an estimator for a n1981-type bias.} However, this stronger assumption is difficult to motivate in practice. Therefore, fw2016, fw2018a recommend to check the sensitivity of the estimates using different values of $L \in \{0, \ldots ,4\}$.
To possibly further improve its finite sample properties, the analytical bias correction can be further iterated. More precisely, we start with an initial $\tilde{\boldsymbol{\beta}}$, afterwards we recompute $\hat{\mathbf{b}}^{\text{abc}}$, update $\tilde{\boldsymbol{\beta}}$, and repeat this procedure a finite number of times. ah2007 refer to this approach as infinitely repeated analytical bias correction.
Next we describe how the split-panel jackknife can be used to form estimators of the asymptotic bias. The idea is to split the panel into sub panels and use them to construct an estimator of the asymptotic bias from different sub panel estimators. We consider two different splitting strategies: the first strategy (SPJ1) is described in fw2016 and the second one (SPJ2) in cfw2017. Let
be estimators of the asymptotic bias, where
$\lceil \cdot \rceil$ and $\lfloor \cdot \rfloor$ are floor and ceiling functions, and the subscript in curly brackets indicates the condition to construct the corresponding sub panel. For clarification, $\{i \leq \lceil N / 2 \rceil \wedge t \leq \lceil T / 2 \rceil\}$ means that the corresponding sub panel only contains the first half of all individuals in the first half of the observation period. In the presence of missing data, we follow the suggestion of fw2018a and ignore the attrition process. Note that, contrary to analytical bias corrections, the split-panel jackknife requires an additional unconditional homogeneity assumption (see assumption 4.3 in fw2016). For instance, this condition rules out trends or structural breaks in the explanatory variables. Further note that both splitting strategies can lead to overlapping sub panels that introduce an additional variance inflation as pointed out by dj2015.\footnote{dj2015 show how to construct non-overlapping sub panels.}
The analytical bias corrections can also be applied at the level of score. The corresponding bias-corrected estimator is the solution to the following system of equations:
fw2016 also suggest a continuously updated score correction where $\widehat{\mathbf{B}}^{\beta}$ and $\widehat{\mathbf{C}}^{\beta}$ are replaced by $\widehat{\mathbf{B}}^{\beta}(\boldsymbol{\beta})$ and $\widehat{\mathbf{C}}^{\beta}(\boldsymbol{\beta})$ such that $\boldsymbol{\beta}$ and the asymptotic bias are estimated simultaneously.
Finally we describe the bias corrections for the average partial effects. A bias-corrected estimator for $\boldsymbol{\delta}$ is
where $\hat{\mathbf{b}}^{\delta}$ is an estimator of the asymptotic bias such that $\tilde{\boldsymbol{\delta}} \operatorname{\overset{a}{\sim}\,} \operatorname{\mathcal{N}} (\boldsymbol{\delta}, \mathbf{V}^{\delta})$. Again, we can either use explicit expressions or we can use the split-panel jackknife described earlier to construct an estimator of the asymptotic bias. Because the application of the split-panel jackknife is generic and already known from the estimation of $\boldsymbol{\beta}$, we omit it for brevity
In the following we assume that $\hat{\boldsymbol{\Delta}}_{it}$ and $\hat{\boldsymbol{\delta}}$ are constructed from $\tilde{\boldsymbol{\beta}}$ and $\tilde{\boldsymbol{\pi}}$, where
An analytical estimator of the asymptotic bias is
where
and $\widehat{\boldsymbol{\Psi}}_{it} = - \partial_{\eta} \hat{\boldsymbol{\Delta}}_{it} / \hat{\omega}_{it}$. Again, by making the stronger strict exogeneity assumption, we can set $L = 0$ and drop the last term in $\widehat{\mathbf{B}}^{\delta}$. Let $\bar{\hat{\boldsymbol{\Delta}}}_{it} = \hat{\boldsymbol{\Delta}}_{it} - \hat{\boldsymbol{\delta}}$, the corresponding estimator of the asymptotic covariance is
where
The first and second term measure the uncertainty caused by the substitution of population by sample means and by the estimation of the structural parameters. The last term is a covariance between both sources of uncertainty that can be dropped if we make the stronger strict exogeneity assumption.\footnote{The estimator of the asymptotic covariance can be adjusted to take into account additional sampling assumptions with respect to the unobserved effects. For instance, if $\{\alpha_{i}\}_{N}$ and $\{\gamma_{t}\}_{T}$ are assumed to be sequences of independent random variables, where $\alpha_{i} \perp \gamma_{t} \, \forall \, i, t$, the estimator of the asymptotic covariance becomes
}
So far we learned how to mitigate the inference problem. In the next subsection we show how the computational costs associated with the estimation and the application of bias corrections can be reduced substantially.
In applications where both panel dimensions are large, estimation with standard software quickly becomes very time-consuming or even infeasible. However, recent advances in computational econometrics embed special solvers in the optimization algorithm of the maximum likelihood estimator to address this problem (see gp2010 and s2018). We first summarize the basic idea of s2018, who proposed an algorithm that can be interpreted as a generalization of g2004 to more than one fixed effect, and then present an extension that can be used to estimate $\boldsymbol{\pi}$ for a given $\tilde{\boldsymbol{\beta}}$.\footnote{To be more precise, s2016 show that the algorithm of g2004 can also be derived using the Frisch-Waugh-Lovell theorem. s2018 combines the resulting projections with h1962's method of alternating projections, leading to a very efficient algorithm for any number of fixed effects.} This extension is necessary e. g. for the bias corrections of the average partial effects. Details about the derivation, the algorithms, and an example code are included in the lementary material.
We start with the estimation of $\boldsymbol{\beta}$. Each step of the optimization algorithm involves solving a weighted least squares problem. More precisely, in iteration $r$ we set
where $\mathbf{Z} = (\mathbf{X}, \mathbf{D})$ and $\mathbf{w}^{[r]} = \boldsymbol{\nu}^{[r]} + \boldsymbol{\eta}^{[r]}$. Remember, $\boldsymbol{\Omega}^{[r]}$ is a diagonal matrix with $\boldsymbol{\omega}^{[r]} = \operatorname{\text{diag}}(\boldsymbol{\Omega}^{[r]})$ and $\boldsymbol{\eta}^{[r]}$ is the collection of linear indexes. Because the rank of $\mathbf{D}$ increases with the sample size, solving the optimization problem quickly becomes infeasible. However we can formulate an alternative weighted least squares problem based on “demeaned” variables so that
Afterwards we can use
to update $\boldsymbol{\omega}$ and $\mathbf{w}$ for the subsequent iteration. Consequently, $\boldsymbol{\beta}$ are the coefficients obtained by a regression of a weighted two-way demeaned $\mathbf{w}$ on a weighted two-way demeaned $\mathbf{X}$ using $\boldsymbol{\omega}$ as weights. $\boldsymbol{\eta}$ are the corresponding fitted values. Thus we can update $\boldsymbol{\beta}$, $\mathbf{w}$, and $\boldsymbol{\omega}$ in each iteration without explicitly updating the incidental parameters. The practical advantage is that $\boldsymbol{\beta}$ and $\boldsymbol{\eta}$ can be very efficiently computed with any software routine developed for weighted least squares problems with high-dimensional fixed effects.\footnote{Some examples available in popular statistical software are reghdfe by c2016 for Stata, lfe by g2013a for R, FixedEffectModels by Matthieu Gomez for Julia, and \textit{pyhdfe} by Jeff Gortmaker for \textit{Python}.} Note that we can also use the same software routines to compute all terms that depend on $\operatorname{\mathbb{M}}$ and $\operatorname{\mathbb{P}}$, as in some expressions of the bias corrections.
Next, we modify the algorithm of s2018 to estimate $\mathbf{D}\boldsymbol{\pi}$ given a fixed $\tilde{\boldsymbol{\beta}}$. Note that we estimate $\mathbf{D}\boldsymbol{\pi}$ instead of $\boldsymbol{\pi}$ as this is simpler and fully sufficient to compute the linear index needed for $F_{it}$ and its derivatives. Again, we start with the weighted least squares problem in iteration $r$:
where $\tilde{\mathbf{w}}^{[r]} = \mathbf{w}^{[r]} - \mathbf{X} \tilde{\boldsymbol{\beta}}$ and $\boldsymbol{\eta}^{[r]} = \mathbf{X}\tilde{\boldsymbol{\beta}} + \mathbf{D}\boldsymbol{\pi}^{[r]}$. Left multiplying by $\mathbf{D}$ yields
and reveals that we only need to demean $\tilde{\mathbf{w}}^{[r]}$ to update the linear index.
Finally, we give a short impression about the capabilities of the algorithms presented in this section. Therefore, we estimate the specification from our empirical illustration using a standard routine and the alternative algorithm we suggest. The former requires more than twelve hours whereas our suggestion only needs three seconds.
In this section, we extend the analysis of fw2016 by two aspects. First, we compare the various analytical and split-panel jackknife bias corrections discussed earlier in this article with respect to their finite sample properties in balanced panels. Second, we analyse whether the improved inference of bias corrections, which is well studied in balanced panels, also shows up in unbalanced panels. For brevity, we restrict ourselves to the analysis of dynamic designs, because fw2016 have already shown that the results with respect to exogenous regressors are very similar for static and dynamic designs in balanced panels.
Besides the asymptotically biased maximum likelihood estimator (MLE), we consider four different analytical bias corrections for the structural parameters. Two of them are applied at the level of the estimator, whereas the others are the solution of modified score equations. ABC1 is the analytical bias correction analysed in fw2016. ABC2 is essentially ABC1, but additionally iterated until convergence. ABC3 and ABC4 are estimators at the level of the score, whereas the latter recomputes the asymptotic bias in each iteration of the non-linear solver. The analytical bias-corrected estimators of the average partial effects are labelled analogously. Further we analyse the two different splitting strategies (\textit{SPJ1--2}) for the split-panel jackknife bias correction and an alternative estimator for the average partial effects (\textit{LPM}) based on the bias-corrected ordinary least squares estimator proposed by hk2002 and the truncated spectral density estimator of hk2007.\footnote{hm2006 show that the bias correction of hk2002 can be applied to dynamic linear models with individual and time specific effects.}
We adapt the dynamic design of fw2016 to allow for the possibility of missing data:
where $\operatorname{\mathbf{1}}[\cdot]$ is an indicator function, $\alpha_{i}, \gamma_{t} \sim \operatorname{\text{iid.}\,} \operatorname{\mathcal{N}}(0, 1 / 16)$, $\epsilon_{it} \sim \operatorname{\text{iid.}\,} \operatorname{\mathcal{N}}(0, 1)$, $\nu_{it} \sim \operatorname{\text{iid.}\,} \operatorname{\mathcal{N}}(0, 0.5)$, and $x_{i0} \sim \operatorname{\text{iid.}\,} \operatorname{\mathcal{N}}(0, 1)$. The corresponding parameters are $\rho = 0.5$ and $\beta = 1$. We consider balanced and unbalanced panels with sample sizes that reflect commonly used panel data sets ($N \gg T$). More specifically, we consider two different patterns of randomly missing data that mimic common situations where some people drop out of a survey and are replaced if necessary. To describe the different patterns of missing data, we distinguish between two types of individuals: type 1 and type 2. The former drop out, whereas the latter are observed over the entire time horizon. To be more precise, let $N_{1}$ and $N_{2}$ be the number of type 1 and type 2 individuals in the unbalanced panel such that $N = N_{1} + N_{2}$. Likewise $T_{1}$ and $T_{2}$ denote the number of consecutive points in time such that $T_{1} < T_{2}$. Both patterns of missing data differ only in the starting point ($s_{i}$) of the time series of each type 1 individual. We set $s_{i} = 1$ in Pattern 1 and sample $s_{i}$ with equal probability from $\{1, \ldots, T_{2} - T_{1} + 1\}$ in \textit{Pattern 2}. For clarification, we set $s_{i} = 1$ for all \textit{type 2} individuals irrespective of the pattern. Figure (ref) provides a graphical illustration of both patterns.
We generate balanced panel data sets with $N = 200$ and $T_{i} = T \in \{15, 20, 25\}$ and unbalanced panel data sets with $(N_{1}, N_{2}) \in \{(300, 100), (150, 150), (60, 180)\}$, $T_{1} = 10$, and $T_{2} = 30$. The different pairs $(N_{1}, N_{2})$ are chosen such that the average number of individuals ($\overline{N}$) and points in time ($\overline{T}$) are $\overline{N} = 200$ and $\overline{T} \in \{15, 20, 25\}$.
To analyse the finite sample performance of the different estimators, we focus on biases relative to the truth and empirical coverage probabilities of 95 % confidence intervals. The latter statistic is especially important, because even if the relative bias is quite small, it might still be large compared to the dispersion of the estimator, with severe consequences for inference. The MLE, ABC1--4, and SPJ1--2 standard errors of the average partial effects are computed using the expression in footnote (ref), which takes the independent sampling of the unobserved effects into account. The LPM standard errors are based on the cluster-robust covariance estimator of cgm2011 to deal with within-individual and within-time correlation of the error terms induced by the probit data generating process. Because there is no obvious way to choose an optimal bandwidth for the estimation of the spectral expectations, we try different choices from a set of values and then report only the results of the choice with the best finite sample performance. We choose $L$ from $\{1, \ldots, 4\}$ for ABC1--4, as suggested by fw2016, fw2018a, and from $\{1, \ldots, \overline{T} - 1\}$ for LPM. All results are based on 1,000 replications.\footnote{We use the \textit{lfe} package of g2013a for the estimation of linear fixed effects models, the non-linear equations solver (\textit{nleqslv}) of h2018 for analytical bias corrections at the level of the scores, and \texttt{R} version 4.0.2 r2020.}
We start by comparing the different analytical bias corrections in a balanced panel. Table (ref) reports the relative biases and coverage probabilities of the estimators for the structural parameters.
As expected, all corrections reduce a larger fraction of the bias and improve coverage as $T$ increases. The difference between the various estimators is in most cases negligible small. Interestingly, ABC2 does not perform better than ABC1, which is remarkable, because often, as for instance in ah2007 and fw2018a, it is noted, that a further iteration of ABC1 could improve the finite sample performance of the estimator.\footnote{j2015 analysed an iterated analytical bias correction, similar to ABC2, for a static design with only individual unobserved effects with similar findings.} ABC3 and ABC4 perform equally well. Our results for the average partial effects are very similar and provided in Table (ref) of the Appendix.
Next we compare the two different split-panel jackknife estimators for the structural parameters in a balanced panel. Table (ref) reports the relative biases and coverage probabilities.
Similar to the analytical bias corrections, we find improved finite sample properties as $T$ increases. One exception is the bias of the estimators in $T = 20$. However, although the bias increases, coverage does not decrease. Overall, we find almost identical properties of both estimators, which is remarkable, because we would expect SPJ2 to have higher dispersion, as significantly smaller sub panels are used to estimate the asymptotic bias. Results for the average partial effects are similar and reported in Table (ref) of the Appendix.
Now we analyse whether the finite sample properties of the different estimators are affected by the two patterns of randomly missing data. Because we already found that the finite sample properties of the different analytical and split-panel jackknife bias corrections barely differ among themselves, we restrict our final analysis to MLE, ABC1, SPJ1, and LPM. Table (ref) and (ref) report the relative biases and coverage probabilities of the estimators for the structural parameters and average partial effects in balanced and unbalanced panels.
We start with the analysis of the finite sample properties in a balanced panel, as these will serve as a benchmark for the properties in unbalanced panels. First, we find that the estimators for the effects of $y_{it - 1}$ are worse than those for $x_{it}$. For instance, we observe larger biases and coverage probabilities that are further away from their nominal level. The distortions in the coefficients are also apparent in the estimates of the average partial effects, which is in contrast to the negligible small biases in the average partial effects of $x_{it}$.\footnote{hn2004, f2009, and fw2016 have similar findings for average partial effects of an exogenous regressor in balanced panels.} In general, the bias corrections perform well in reducing the relative biases and improving coverage. As in fw2016, the properties of SPJ1 are worse than those of ABC1. LPM as an alternative estimator for the average partial effects works well too.\footnote{f2009 has similar findings for linear probability models with only individual fixed effects.} Interestingly, LPM tends to overestimate the average partial effects of $y_{it - 1}$, while the other estimators underestimate them. Further, the optimal bandwidths of LPM for the estimation of the spectral expectations are much larger than for ABC1 and increase rapidly with $T$. This indicates that the temporal dependence induced by the probit data generating process can be very strong, which in turn makes the choice of appropriate bandwidths for \textit{LPM} difficult in practice. Next we compare the finite sample properties of the different estimators in unbalanced panels with our benchmark. First, by comparing the relative biases of \textit{MLE}, we can confirm the conjecture of fw2018a that the magnitude of the asymptotic bias in the limiting distribution depends on $\overline{N}$ and $\overline{T}$. However, the coverage is worse because the sample size of an unbalanced panel is larger than that of a balanced panel (with equal $\overline{N}$ and $\overline{T}$), which in turn implies a smaller standard deviation of the estimators and thus distortions that are larger relative to the variance of the estimator. Second, compared to the benchmark we notice some substantial differences in the performance of \textit{SPJ1}. Whereas the properties of \textit{ABC1} and \textit{LPM} are unaffected by the different patterns of missing data, they are partially significantly worse for \textit{SPJ1}, especially for $T = \overline{T} = 15$. \textit{Pattern 1} stands out in particular, because it clearly shows that the reduction of bias and improvement of coverage are deteriorating. An intuitive explanation is that the splitting strategy leads to sub panels of widely differing sizes. This issue is not so severe in \textit{Pattern 2}, but the performance is still worse than the benchmark.
Finally, we briefly summarize the key findings. We find no differences in the finite sample performance when we compare the various analytically bias-corrected and the different split-panel jackknife estimators among themselves. Although the latter have the advantage that they are relatively easy to implement, this convenience is associated with some performance losses. More precisely, split-panel jackknife estimators have higher distortion and react sensitive to different patterns of randomly missing data, whereas the performance of the analytical bias corrections is unaffected. An alternative estimator for the average partial effects based on the bias-corrected ordinary least squares estimator works well too, but to find an appropriate bandwidth for the required spectral density estimator might be challenging in practice.
In the next section, we demonstrate the usefulness of computational advances and bias corrections with an empirical example from labour economics.
We illustrate one possible area of application by analysing the inter-temporal labour force participation of women using longitudinal micro data (1984--2013) from the German Socio Economic Panel (GSOEP). More specifically, we want to investigate how fertility decisions and non-labour income jointly affect women's labour force participation using an unbalanced panel data set of 6,241 women in labour force observed consecutively for at least ten years. Further details about the sample are provided in the Appendix.
In spirit of h1999, we estimate the following model specification:
where $y_{it}$ is an indicator equal to one if woman $i$ participates in the labour force at time $t$, $\mathbf{x}_{it}$ is a vector of explanatory and further control variables, $\boldsymbol{\beta}$ are the corresponding common parameters, $\alpha_{i}$ and $\gamma_{t}$ are unobserved effects that capture individual specific taste for labour and permanent income as well as control for the business cycle and other time specific shifts in preferences, and $e_{it}$ is an idiosyncratic error term independent of $\mathbf{X}_{i}^{\star}$ and $\boldsymbol{\pi}$ with mean zero. As our data set is unbalanced, we have to additionally assume that conditional on $\mathbf{X}_{i}^{\star}$ and $\boldsymbol{\pi}$ the attrition process is random. We consider the following explanatory variables: number of children in different age groups, various non-labour income classes, and an indicator that is equal to one if a birth occurs in the next year. Further control variables are squared age, marital status, a regional identifier for Eastern Germany, number of children between zero and one in the previous year, and number of other household members.
Table (ref) shows some descriptive statistics of our sample.
The average participation rate is 73 % in the full sample and 67 % for the group of movers who change their labour force participation decision at least once. The group of women who never participate is the smallest and most different from the other groups. On average, this group is older, more likely to be married, and lives in Western Germany. Contrary, women who always participate have less children and live in smaller households. The full sample comprises 97,465 observations and consists of 6,241 women observed for a maximum of 28 years. As some households drop out of the GSOEP and are replaced by new ones, this leads to a pattern of missing data similar to Pattern 2 in the simulation study. On average, we observe each woman for roughly 16 years, or each year we observe about 3,481 women.
We consider the following estimators for the structural parameters and average partial effects: MLE, ABC1, SPJ1, and LPM. The estimators are labelled and bandwidths are chosen as in the simulation study. Table (ref) reports the corresponding estimates.
All results are intuitive and in line with the theoretical model of h1999. We find positive state dependence and negative effects of transitory non-labour income, number of children, and expectations about future fertility. All effects are significant at the 5 % level. Remarkably, most probit estimates of the average partial effects are very close to each other. Exceptions are those with respect to lagged participation which range from 0.23 up to 0.31. LPM estimates are also quite close except for lagged participation and number of children between zero and one. Overall we find evidence for strong state dependence. A woman who has currently a job increases her probability to participate in the future by 23--54 percentage points. Further we find that women respond heterogeneously to changes in transitory non-labour income. Being in the middle class reduces the participation probability by roughly one percentage point compared to a woman in the lower income class. The reduction associated with belonging to the upper income class is significantly stronger with three up to five percentage points. Finally we find that the number of children reduces the likelihood of participation substantially. As expected, the effect is declining in age of children. Each additional child between zero and one reduces the probability to participation 20 up to 30 percentage points. For children older than four, the reduction is only one percentage point. The results are largely consistent with the empirical findings of h1999. However, contrary to him, we find that future birth always negatively affects current participation decision irrespective of the chosen estimator. This might support the author's perfect foresight assumption with respect to life-cycle fertility decisions.
Finally, we check the sensitivity of ABC1 and LPM to different bandwidth choices and conduct a simulation study calibrated to our empirical illustration. The results are reported in Table (ref) and (ref) of the Appendix. With respect to the different bandwidth choices, we find that especially the ABC1 estimates are very robust. We find the largest variation with respect to the effect of lagged participation. The other effects are almost indistinguishable. Most LPM estimates are robust as well. Exceptions are the effects of lagged participation, being in the upper class, and number of children between zero and one. The results of the calibrated simulation study confirm that the finite sample properties of MLE can be improved. However, the improvement is less good compared to the simulation experiments with synthetic data. The performance of SPJ1 and \textit{LPM} is in some cases significantly worse than that of \textit{ABC1}. Given that there have been several labour market reforms with heterogeneous effects on different sub populations that most likely violate the unconditional homogeneity assumption of \textit{SPJ1}, this could partly explain its worse performance.
In this article, we offer new relief and guidance for empirical researchers by showing how popular binary choice estimators benefit from recent advances in econometrics. Especially the analytically bias-corrected estimator of fw2016 convinced by its good performance in all cases.
Although we have focused on panel data binary choice models, we would like to point out that the bias corrections derived by fw2016 are far more general. First, they can also be used to reduce the asymptotic bias of other popular non-linear maximum likelihood estimators e. g. poisson and tobit. The algorithms described in this article can be easily adapted to these problems. Second, the bias corrections can also be applied if we observe cross-sections of networked activities instead of panels e. g. a cross-section of bilateral trade flows. For instance, cfw2017 use the bias corrections of fw2016 to mitigate the asymptotic bias problem in a hmr2008-type model to determine the extensive margin of trade.
Recently, there are also some extensions of fw2016. For instance, wz2020 and hsw2020 have extended these bias corrections to a special three-way error component that is particularly relevant for the estimation of the intensive and extensive margin of trade in a panel of bilateral trade flows. Further, cfw2020 use multiple binary choice regressions to estimate the distribution of non-binary outcomes conditional on strictly exogenous regressors and two unobserved effects. The corresponding inference problem is addressed using an analytical bias correction. These extensions can also benefit from the findings of this article.
Finally, future research could investigate whether bootstrap procedures can further improve inference. In an additional simulation study we find that bootstrapping a bias-corrected estimator in spirit of k2014 slightly improves the coverage of the estimators for the structural parameters but not for the average partial effects, which questions the validity of this bootstrap procedure in our specific setting.
\printbibliography