EconBase
← Back to paper

Grouped fixed effects regularization for binary choice models

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.

72,729 characters · 12 sections · 74 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Grouped fixed effects regularization for binary choice models

\thispagestyle{empty}

abstractWe study the application of the grouped fixed effects approach to binary choice models for panel data in presence of severe complete separation. Through data loss, complete separation may lead to biased estimates of Average Partial Effects and imprecise inference. Moreover, forecasts are not available for units without variability in the response configuration. The grouped fixed effects approach discretizes unobserved heterogeneity via k-means clustering, thus reducing the number of fixed effects to estimate. This regularization reduces complete separation, since it relies on within-cluster rather than within-subject response transitions. Drawing from asymptotic theory for the APEs, we propose choosing a number of groups such that clustering delivers a good approximation of the latent trait while keeping the incidental parameters problem under control. The simulation results show that the proposed approach delivers unbiased estimates and reliable inference for the APEs. Two empirical applications illustrate the sensitivity of the results to the choice of the number of groups and how nontrivial forecasts for a much larger number of units can be obtained.

\vskip 3mm {\bf Keywords:} {\sc Average Partial Effects, Dynamic models, Grouped fixed effects,\\ Rare events, Regularization } \vskip 3mm {\bf JEL Classification: C13,C23, C25} {\sc }

\setcounter{page}{0}

Introduction

Fixed-effects (FE) binary choice models are prominently used in applied econometrics and popular examples arise from a wide range of applications.\footnote{Noteworthy examples come from labor market participation heckman1980does with a focus on fertility choices for female married workers hyslop1999state, self-reported health status contoyannis2004dynamics, transitions in income dynamics cappellari2004modelling, household finance alessie2004ownership, and drivers of unionization choices wooldridge2005simple. Recent applications can be found in studies on firms' behavior in accessing credit pigini2016state, migrants’ remitting choices bettin2018dynamic, link formation models dzemski2019, energy poverty drescher2021determinants,persistence of innovation in firms arroyabe2022estimation.} Estimation of model parameters in this context, where one or more sets of FE are included, is usually carried out by Maximum Likelihood (ML).

With binary panel data, it may happen that the dependent variable does not exhibit within-subject variation when outcomes describe highly persistent phenomena or extremely rare events, such as employment status and the occurrence of financial crises. Practitioners estimating FE binary choice models based on these data are often unable to recover finite estimates of all the individual intercepts, an issue known in the literature as the complete separation (CS) problem albert1984existence: for a given subject, if their outcome configuration does not vary over time, the log-likelihood will be monotone in their intercept, leading to a non-finite ML estimate of their FE. Although the true population intercept is finite, the observed time series might not be long enough to observe a time-varying configuration.\footnote{To give the dimension of the CS problem in typical settings, kunz2021predicting describes an application on health care utilization where 29% to 45% of subjects do not exhibit outcome variation over time. In the study on labor market participation in dj2015 and fernandez2009fixed, revisited here in Section (ref), 60% of the observations are in CS.}

Statistical software typically remove subjects in CS, which, in absence of cross-sectional dependence, has no direct effects on the regression parameter estimates. The sub-sampling, however, impacts other quantities of interest in three main respects: i) the Average Partial Effects (APEs) are overestimated, as the subjects dropped from the dataset are likely to have small individual population partial effects, due to their large index functions; ii) as APEs converge at the rate $1/\sqrt{N}$, with $N$ being the number of subjects, the reduced sample size leads to an imprecise large sample approximation of the APE sampling distribution, resulting in poor finite-sample coverage; iii) forecasts for discarded units are trivial, as their predicted probability would always be zero or one, in and out of sample.

This paper motivates the application of the Group Fixed Effects (GFE) approach, put forward by blm2022, in settings where the use of FE leads to pervasive CS. The GFE approach is based on a two-step procedure: in the first step individual, possibly continuous, unobserved heterogeneity (UH) is discretized by {\em k-means} clustering based on the model covariates; in the second step, group-membership indicators enter the main specification as cluster-specific intercepts. The intrinsic regularization introduced by GFE, which limits the number of FE to be estimated, reduces the instances of CS. This happens because the existence of finite estimates for the group-specific intercepts relies on the within-cluster, as opposed to within-subject, variability in outcome configuration. Therefore, subjects without outcome variation are retained if they end up in a cluster together with individuals who exhibit time variability in their response configuration.

We show that the GFE regularization effectively overcomes the finite-sample issues entailed by the CS-related sample size reduction: i) APEs computed using GFE estimates account for the systematically smaller marginal effects of subjects otherwise dropped, providing a more precise quantification of the population APEs; ii) the larger sample size actually used yields more accurate coverage for the APEs; iii) the GFE approach allows one to make non-trivial predictions for units without variation in the response variable, as long as these are clustered in groups where outcome variation at cluster-level is observed.

Ways of dealing with CS are the subject of the stream of literature that relies on shrinkage to obtain finite ML estimates. These approaches are inspired by the modified score correction introduced by firth1993bias for the logit model, applied to handle CS in cross-section data by Heinze2002 and Heinze2006, then generalized by Kosmidis2009 to nonlinear models of the exponential family. Modified versions of this approach have later been used to shrink FE estimates in binary choice models by kunz2021predicting and pigini2021penalized, who focus on forecasts, and Cook2018, who suggest FE shrinkage to reliably quantify population APEs by means of a plug-in estimator. Despite the conjecture put forward by Cook2018, a thorough study of finite-sample properties nor complete asymptotic theory for APEs with a Firth-type shrinkage is available. For instance, it is well known in the literature that plug-in APE estimators still suffer from the typical incidental parameters problem, which might not be negligible when the individual time series is short, thus requiring a bias correction dj2015.

Further to providing evidence of better coverage of the GFE plug-in APE estimator, we show that the cluster regularization employed by the proposed approach can be used to limit the effects of the incidental parameters problem on the APE estimator in finite samples. Relying on the asymptotic properties of the proposed estimator, we provide the practitioner with a guideline to choose a number of groups that simultaneously makes the incidental parameters bias and the approximation error entailed by discretization both negligible in finite samples. Therefore, no further bias reduction is required. Finally, it is worth to stress that the GFE approach can directly be applied to dynamic binary choice models, differently from the shrinkage-type estimators that would require a modification of the score correction term.

The simulation study analyzes the finite-sample properties of the GFE plug-in estimator of the APE, for both static and dynamic logit models, in presence of moderate to severe degrees of CS. The results show that the GFE approach mitigates the APEs overestimation, which would otherwise result from dropping subjects in CS, as witnessed by the performance of the infeasible estimator. The performance of the proposed approach is also compared to the APEs plug-in estimators obtained using ML and to the analytical and jackknife bias-corrected APE estimators hahn2011bias,dj2015. By discarding significantly fewer observations, the GFE APE estimator exhibits minimal bias and better empirical coverage. Moreover, choosing the number of groups approximating individual UH according to the proposed rule makes the incidental parameters bias negligible in finite samples, signaling that further bias reduction can be avoided.\footnote{We also explore the performance of the plug-in estimator based on the Firth-type score correction firth1993bias. While this approach does not lead to loss of observations, the shrinkage of the FE estimates does not seem to be effective as a bias reduction device.}

commentprovide imprecise estimates for APE, with misleading inference based on poor coverage. In particular, in the static case we show how CS leads to the overestimation of APE.\footnote{Despite being present also in the dynamic setting, the overestimation of APE does not emerge in simulation due to the opposite effect of Nickell’s bias, which leads to a compensation, see Section (ref) for details.}

We present the results of two real-data applications. The first revisits the empirical application on the participation of young working women in the labor market proposed, among others, by dj2015 and fernandez2009fixed. In this setting, CS involves around 60% of the original sample due to the strong intertemporal correlation of employment status, a phenomenon often observed in labor market studies. We show that, as in the simulation study, the GFE approach retains a larger portion of the dataset and leads to a quantification of APEs that coherently lies between the pooled and the ML-based bias-corrected estimators. The second application presents a forecast exercise based on rare events. We use the panel data on financial crises issued by laeven2018systemic, where the dependent variable is equal to one if a country in a particular year witnessed financial turmoil. We show that the GFE approach manages to offer non-zero predicted probabilities for a higher number of countries with respect to ML alternatives and has a good forecasting performance.

commentOverall, the simulation study and the empirical application suggest the use of GFE approach when the degree of CS is particularly high.

The rest of the paper is organized as follows: Section (ref) outlines the effects of CS in ML estimation of FE binary choice models and motivates the use of the GFE approach; Section (ref) presents the simulation study; Section (ref) illustrates the two empirical applications. Finally, Section (ref) concludes.

Econometric methods

Background on fixed-effects binary choice models

For $i=1,\ldots,N$ and $t=1,\ldots,T$, we study the model

equation[equation omitted — 110 chars of source]

where $ \mathbbm{1}(\cdot)$ is the indicator function, $x_{it}$ denotes a set of $J$ individual-specific covariates associated with a conformable vector of unknown parameters $\beta_0$ and may include $y_{i,t-1}$; $\alpha_{i0}$ parameterizes the UH as time-invariant individual effects, while $u_{it}$ is the i.i.d error, whose distribution is either standard logistic or normal.

The structural and nuisance parameters in model (ref) can be jointly estimated using ML, leading to $(\hat{\beta}', \hat{\alpha}_1, \ldots,\hat{\alpha}_N)'$. As is well known, the ML estimator suffers from the so-called incidental parameters problem (IPP), which is due to the estimation noise introduced by the nuisance parameters entering the profile likelihood for the structural ones NS1948. The IPP leads to an asymptotic bias in the limiting distribution, even if both $N$ and $T$ $\to \infty$, but in a fixed proportion to each other.\footnote{This framework is referred to as rectangular array asymptotics Li2003, where $N,T \to \infty$ with $N/T \to \rho$, $0 <\rho < \infty$.} Bias reduction techniques for the ML estimator are available, in the form of both analytical fernandez2009fixed, hahn2011bias and jackknife dj2015 corrections.

The objects of interest in binary choice models are usually the APEs. Let us define the population APE as

equation[equation omitted — 80 chars of source]

where $\mu_{it}(\beta_0,\alpha_{i0}) = F'(x_{it}'\beta_0 + \alpha_{i0}) \beta_{0}$ and $F'(\cdot)$ is the first derivative of the probit/logit link function. The ML plug-in estimator of $\mu_0$ is readily available as

equation[equation omitted — 111 chars of source]

and its asymptotic expansion is such that $\hat{\mu} = \mu_0 + O_p(1/T)$, where the $O_p()$ term represents the bias arising from IPP hahn2011bias. Unlike the ML estimator of $\beta_0$, any plug-in APE estimator does not converge at the rate $1/\sqrt{NT}$, but more slowly, as stated by Theorem 5.1 by dj2015. Define $\mu_i = T^{-1}\sum_t \mu_{it}(\beta_0,\alpha_{i0})$ and $\sigma^2_\mu = \underset{N \to \infty}{\mathrm{lim}} N^{-1} \sum_{i=1}^N (\mu_i - \mu_0)^2$. Then they show that as $N,T \to$ $\infty$ with $N/T \to \rho$, $0 <\rho < \infty$, we have:

equation*[equation* omitted — 124 chars of source]

The above expression clarifies that the plug-in ML APE estimator converges at the rate $1/\sqrt{N}$ and the IPP bias, now captured by the term $O_p(1/\sqrt{T})$, is now asymptotically negligible, as it vanishes as $T \to \infty$. However, this bias may still be present in finite samples, especially when the observed time series is short or the IPP is particularly severe (e.g., Nickell's bias in dynamic models). Therefore, the use of analytical or jackknife bias corrections is advised for APEs fernandez2009fixed, dj2015.

Our main concern in this context is the FE estimate $\hat{\alpha}_i$ in finite samples.

commentDefine $y_i=(y_{i1},\ldots, y_{iT})$ as the vector for the response variable configuration of individual $i$.

Whenever $\sum_{t=1}^T y_{it}=0$ or $T$, meaning that there is no variability in the dependent variable, the ML estimate of $\alpha_{i0}$ does not exist finite, which is an instance of CS. \begin {example} As an example, consider a static FE logit model without covariates: it is easy to see that the individual likelihood $\ell_i= \alpha_i\sum_{t}^{T}y_{it} - T\log\left[1+\exp(\alpha_i)\right]$ is maximized at $ \hat{\alpha}_i=\log(\frac{p_i^*}{1 - p_i^*})$, where $ p_i^*=\sum_{t}^{T} y_{it}/T$. Therefore, the ML estimate of the individual intercept is not finite when $p_i^*$ is either $0$ or $1$. \end {example} Statistical software usually removes subjects in CS from the dataset. Although this reduction has no impact on the estimates of structural parameters $\beta_0$ in absence of cross-sectional dependence, the quantities computed using the predicted probabilities exhibit a bias that depends on the intensity of the CS problem.

commentConsider the predicted probabilities: \[ p_{it} =F\left(x_{it}'\beta_0 + \alpha_{i0}\right) \] where $F(\cdot)$ is the probit/logit link function.

Non-finite estimates of $\alpha_{i0}$ lead to an estimated probability $F(x_{it}'\hat{\beta} + \hat{\alpha}_{i})$ exactly equal to zero or one. Consider expression (ref) in presence of CS: \[ \widehat{\mu}^\ast = \frac{1}{N^\ast T} \sum_{i \not\in D}\sum_{t}\mu_{it}(\hat{\beta}, \hat{\alpha_i}), \] where $D=\{ i: \sum_{t}^Ty_{it} = 0 \,\ \text{or} \sum_{t}^Ty_{it} = T \}$ is the subset of individuals for whom no transitions are observed in the outcome variable and $N^\ast < N$ the number of individuals who are not in CS. Data reduction causes the overestimation of $\mu_0$, because discarded units tend to have a large index functions in absolute value and, in turn, small individual population partial effects, usually close to zero. After the removal of problematic units, the distribution of the estimated PEs (in absolute value) becomes left-truncated as the smaller values are excluded. Because $\widehat{\mu}^\ast$ is computed using only the PEs of individuals who do not belong to $D$, the APE, conditional on this restricted sample, is systematically greater than $\mu_0$. The resulting quantification of the effects of interest is then imprecise.

exampleConsider a static logit model including a single binary explanatory variable $x_{it}$ and $T=2$. The PE for a generic individual $i$ is \[ PE_i(\alpha_{i0}, \beta_0)= F(\alpha_{i0}+ \beta_0)-F(\alpha_{i0}). \] Consider the case in which $x_i = \{x_{i1},x_{i2}\} = \{0,1\}$. Conditional on $x_i$, the probability of not observing a change in outcomes is \begin{gather*} P(i \in D \mid\alpha_{i0}) = P(y_{i1} = 0, y_{i2} = 0 \mid\alpha_{i0}) + P(y_{i1} = 1, y_{i2} = 1 \mid\alpha_{i0}) = \\ [1-F(\alpha_{i0})][1-F(\alpha_{i0}+\beta_0)] + F(\alpha_{i0})F(\alpha_{i0}+\beta_0). \end{gather*} \begin{figure}[h!b] \caption{ Probability of CS and PE } \begin{tikzpicture} \begin{axis}[ width=11cm, height=6cm, domain=-6:5.5, samples=600, xlabel={$\alpha_{i0}$}, ylabel=, ymin=0, ymax=1, axis lines=left, axis line style={thick}, tick style={thick}, grid=both, major grid style={gray!20}, clip=false, legend style={ at={(0.5,1.05)}, anchor=south, legend columns=1, draw=none, fill=none, font=, inner ysep=10pt } ] \pgfmathdeclarefunction{F}{1}{\pgfmathparse{1/(1+exp(-#1))}} \pgfmathdeclarefunction{psel}{1}{ \pgfmathparse{F(#1)*(1-F(#1+1)) + (1-F(#1))*F(#1+1)} } \pgfmathdeclarefunction{Delta}{1}{ \pgfmathparse{F(#1+1)-F(#1)} } \addplot [very thick, red!70!black, dashed] ({x},{1 - psel(x)}); \addlegendentry{$P(i \in D \mid\alpha_{i0}, \; \beta_0 = 1)$ } \addplot [very thick, blue!70!black] ({x},{Delta(x)}); \addlegendentry{$PE_i(\alpha_{i0}) = F(\alpha_{i0}+1)-F(\alpha_{i0})$} \node[anchor=south west, font=, color=blue!70!black] at (axis cs:-0.7,0.25) ; \node[anchor=south east, font=, color=red!70!black] at (axis cs:-5,0.95) ; \node[anchor=south east, font=, color=red!70!black] at (axis cs:3,0.95) ; \end{axis} \end{tikzpicture} \end{figure} In Figure (ref), we illustrate $PE_i(\alpha_{i0})$ and the probability of being dropped $P(i \in D \mid \alpha_{i0})$ for $\beta_0=1$. As we can see, individuals who are more likely to be dropped due to complete separation correspond precisely to those with extreme fixed effects and negligible partial effects.

Heavy data separation has an additional effect on the estimation of the APEs. Since $\hat{\mu}$ converges to the real APE at the rate $N^{-1/2}$, finite sample performance of the estimator crucially relies on the availability of a large number of individuals. However, unless $N$ is large, when events are extremely rare or outcomes are very persistent, $N^\ast \ll N$ might be too small for the asymptotic approximation to deliver good coverage. In this vein, $\hat{\mu}^\ast$ can also lead to misleading inference.

Finally, ML estimation in presence of CS leads to trivial forecasts for discarded subjects, potentially giving rise to misclassification instances, regardless of the threshold $\tau \in [0,1]$ used to build the test-set confusion matrix. In certain contexts, this is not without consequences: one example is that of rare events (low-probability, high impact). In fact, the predicted probability for a subject will be non-zero only if another event has been experienced by the same unit in the past, thus preventing a meaningful forecast of a first-ever occurrence.

commentConsider, for instance, a scenario in which the rare event is observing $y_{it}=1$, meaning that the ML estimates $\hat{\alpha_i} = - \infty$ for all units for which $\sum_{t}^Ty_{it} = 0$. In formulae: \[ \hat{p}_{it}=\mathbb{I}\left[F(x_{it}'\hat{\beta} - \infty ) > \tau\right] = 0 \,\, \forall \,\, t \,\, \mathrm{and} \,\, \tau, \] where $\tau \in [0,1]$ is a cut-off point. Predictions for units in CS are tied to their past history and lead to biased outcomes whenever the dependent variable changes for the first time in the year of forecast. At the same time, discarding problematic units is often not an option, as analysts require forecast for every unit in the dataset, given also that in presence of rare events the percentage of dropped units is substantial.

Grouped fixed-effects estimation with complete separation

In this Section, we illustrate how the use of the GFE approach blm2022 mitigates the issues arising with CS by limiting the number of subjects dropped due to the lack of outcome variation.

The GFE approach is based on the idea that individual UH $\alpha_{i0}$ can be approximated by a smaller set of group-specific parameters. Grouped structures of heterogeneity, which are assumed to be discrete in the population, are becoming increasingly popular in the FE literature hahnmoon2010,bm2015,lumsdaine2023estimation,mugnier2025simple. In contrast, blm2022's approach is in the same spirit of contributions that employ clustered structures to approximate general forms - both continuous and discrete - of UH beyhum2024inference, freemanweidner2023 and it is, to the best of our knowledge, the only viable for nonlinear models. The GFE estimation procedure consists of two steps:

enumerate• {\bf Classification step} The individual heterogeneity $\alpha_{i0}$ is {\em discretized} by {\em kmeans} clustering, which uses the vector of the $J$ individual averages $\bar{x}_{i} = T^{-1}\sum_{t} x_{it}$. The algorithm partitions individuals into $K$ groups, with $K \ll N$, such that \begin{equation} (\bar{x}_{\hat{k}=1}, \ldots, \bar{x}_{\hat{k}=K},\hat{k}_1, \ldots, \hat{k}_N) = \mathrm{argmin} \sum_{i=1}^N || \bar{x}_i - \bar{x}_{k_i} ||^2, \end{equation} where $\bar{x}_{k}$ is the mean of $\bar{x}_i$ in group $k$. • {\bf Estimation step} Consists of the ML estimation of the model \[ y_{it} = \mathbbm{1}(x_{it}' \beta_0 + \alpha_{\hat{k}_i} + u_{it} > 0), \] where $\alpha_{\hat{k}_i} = \alpha_k \mathbbm{1}(i \in k)$, $k = 1, \ldots,K$. These are the cluster-specific FE, related to the group-membership dummies, and are estimated jointly with the structural parameters yielding $(\tilde{\beta}', \tilde{\alpha}_1, \ldots, \tilde{\alpha}_K)'$.

Regularity conditions for the validity of {\em kmeans} clustering and details on the asymptotic properties of the GFE estimator are given in blm2022, suppblm2022.

It is worth highlighting that the moments used for the {\em kmeans} clustering have to be informative about the UH: for $T \to \infty$, blm2022 clarify that they have to be functions of the UH such that one can separate two individuals with a different level of UH by comparing their vector of moments.\footnote{In particular, Assumption 2 in blm2022 requires moments to be injective. Denote the individual UH as $\xi_{i0}$, of unspecified form, such that $\alpha_{i0} = \alpha(\xi_{i0})$, where $\alpha(\cdot)$ is a Lipschitz-continuous function. Then there exist moments $h_i = (1/T) \textstyle{\sum_t} h(y_{it},x_{it}')$ and a Lipschitz-continuous function $\phi(\cdot)$ such that $\mathrm{plim}_{T \to \infty} h_i = \phi(\xi_{i0})$; moreover, there exist a Lipschitz-continuous function $\psi(\cdot)$ such that $\xi_{i0} = \psi^{-1}(\phi(\xi_{i0}))$.} In principle, not only the moments of the regressors, but also $\bar{y}_i = T^{-1} \sum_{t=1}^T y_{it}$ can be used. However, when outcomes are highly persistent or very rare, $\bar{y}_i$ may exhibit small to no variability between subjects when the time series is short, which hampers its ability to inform us about different types of UH. For this reason, in our subsequent simulation study and empirical applications, we avoid using the individual averages of the dependent variable in clustering.

The smaller number of FE to estimate decreases the likelihood of dealing with a set of individuals clustered in a group where no variability in the outcome variable is observed. In practice, if individuals with no transitions are clustered together with non-problematic ones, the ML estimate of the related group-specific intercept exists finite. As a consequence, the number of observations discarded due to CS is lower.

Example 1 (continued) Consider again the example about the static panel logit model without covariates. Applying the GFE procedure leads to three possible solutions for $\alpha_{\hat{k}}$: \[ \tilde{\alpha}_{k} = \log\left(\frac{p_{\hat{k}}}{1 - p_{\hat{k}}}\right) =

cases- \infty & if \quad p_{\hat{k}}= 0, \\ finite & if \quad p_{\hat{k}} \in (0,1),\\ \infty & if \quad p_{\hat{k}}= 1,

\] where $p_{\hat{k}} = \frac{\sum_{i \in k} \sum_{t} y_{it}}{T\sum_{i}\mathbbm{1}(i \in k)}$ is the average of the outcomes in group $k$. In Table (ref), we consider a panel composed by $N=2$ individuals and $T=2$ time periods. The second individual does not show any state transition, therefore it is not possible to have a finite estimate for $\hat{\alpha}_{2}$. If, instead, the two individuals are clustered in the same group ($k=1$), within-cluster variability allows us to obtain a finite ML estimate of their shared intercept $\tilde{\alpha}_{1}$.}\footnote{In the example $p_{1}=3/4$ and $\tilde{\alpha}_{1}=\log{3}\approx 1.1$.}

table[table omitted — 785 chars of source]

By limiting the number of FE to be estimated and relying on an increased sample size, the GFE approach introduces a regularization that mitigates the consequences of CS. Unlike popular contributions, such as phillips2016 and wang2024panel, which focus on regularization as a way to credibly identify latent grouped patterns, GFE regularization pertains to the estimation of a more parsimonious model.

The amount of information loss depends on the within-cluster variability implied by the classification step, which makes the choice of the number of groups $K$ crucial in this context. At the same time, in the framework put forward by blm2022 clustering is an approximation device for an unspecified form of UH, thus the granularity of the discretization is closely tied to the quality of the approximation. blm2022 propose a rule for choosing $K$ that reflects such trade-off between accuracy and parsimony. Specifically, they set \[ K = \underset{K \geq 1}{\mathrm{min}} \{K: \widehat{Q}(K) \leq \gamma \widehat{V}_{\bar{x}} \}, \] where $\widehat{Q}(K)$ here indicates the {\em kmeans} objective function in (ref) and $\widehat{V}_{\bar{x}}$ is an estimate of the variability between the moment vectors and the individual UH. In practice, they propose choosing the smallest number of groups such that the variability of $\bar{x}_i$ with respect to the centroids is less than or equal to the variability of $\bar{x}_i$ with respect to the individual UH. This choice is also governed by the specification of the user-specified hyper-parameter $\gamma$, which is bound in $(0,1]$ and such that smaller values require a larger number of groups in order to make the between-centroids variability smaller or equal than $\widehat{V}_{\bar{x}}$. However, no tuning procedure for $\gamma$ is suggested. We here provide further guidance on the choice of $K$ on the basis of the asymptotic expansion for the APEs, that are the main objects of interest in the present context.

commentThe number of groups $K$ is not known in practice and must be determined. To address this, blm2022 proposed a rule\footnote{Technical details on the rule are available in Appendix XXX.} that selects the smallest $K$ such that the objective function of the k-means problem (i.e., the variability of $\bar{x}_i$ with respect to the centroids) is less than or equal to the variability of $\bar{x}_i$ with respect to individual UH, multiplied by a user-specified hyperparameter $\gamma \in (0,1]$. The rule is conceived to balance the trade-off between the number of groups, that is, the goodness of the approximation of the UH, and the severity of the IPP bias, as is clear by the role of $K$ in Equation (ref). In empirical settings, the choice of $\gamma$ is not trivial and blm2022 does not provide a recommendation on its selection.

Consider the GFE plug-in APE estimator $\tilde{\mu}$. Then assuming that suitable regularity conditions hold,\footnote{These conditions are contained in Assumptions 1-3 in blm2022 and Assumption S1 in suppblm2022.} the asymptotic expansion of the APE estimator suppblm2022 implies that as $N,K,T\to \infty$, $K/NT \to 0$, with $N/T \to \rho$,

equation[equation omitted — 252 chars of source]

The above expression shows that $\tilde{\mu}$ has three sources of bias, represented by the $O_p(\cdot)$ terms: the $O_p(\sqrt{N}/K^2 )$ term arises from the approximation error due to the discretization of the UH using the kmeans procedure; the $O_p( 1/\sqrt{T})$ term originates as a classification-step IPP bias, due to the use of $N$ averages, $\bar{x}_{i}$, for $NT$ observations; finally, the $O_p(K/T\sqrt{N})$ term represents the second step IPP bias due to the estimation of $K$ cluster-specific intercepts using $NT$ observations.

That all bias terms in expression (ref) will be asymptotically negligible is guaranteed by $K$ growing at certain rates in relation to $N$ and $T$. In particular, blm2022 suggest that setting $K$ proportional to or greater than $\mathrm{min}(\sqrt{T}, N)$ guarantees that the clustering approximation bias will be of order $O_p(1/T)$ for the GFE estimator, thus $O_p(1/\sqrt{T})$ in the above expansion. In this respect, it is worth noting that setting $K$ in this way generates a constant bias in the asymptotic distribution of the GFE estimator blm2022, giving rise to a negligible term in the distribution of the APE GFE estimator. In addition, we argue that the second-step IPP may also be asymptotically negligible as long as $K$ is chosen to be smaller than $T\sqrt{N}$. This refinement of the rule indirectly suggests how to tune $\gamma$, which should be chosen within a range of values yielding $\sqrt{T} \ll K < T\sqrt{N}$.

As the GFE approach limits the number of units dropped from the data, the asymptotic distribution offers a better approximation of the sampling one of $\tilde{\mu}$ in finite samples, with respect to the FE approach, thus providing a more accurate coverage. In addition, $\tilde{\mu}$ is less likely to overestimate $\mu_0$, as more units, including those with small PE, are retained in the dataset. Finally, the GFE approach manages to provide nontrivial predictions for units without transition in the outcome variable, as long as they are clustered in groups where outcome variability is observed. This allows practitioners to get finite predictions for every unit.

commentWe propose a simple procedure to select $\gamma$ which is rooted in the asymptotic behavior of the GFE plug-in APE estimator and suggest setting the parameter to get a number of groups for which: i) $K \gg \sqrt{T}$ to have negligible approximation and first step IPP bias and ii) $K < T\sqrt{N}$ so that second step IPP bias is not severe. The refinement of the rule leads to a range of values of $\gamma$ in which the biases are balanced and further choice within the interval is left to the analyst. In practice, we suggest performing the GFE estimation for different values of $\gamma$ and finding the aptest value of $K$ in the interval.

Simulation study

Static logit model

We study the finite sample performance of the GFE approach by estimating a static logit model in presence of CS. We generate data from the model

equation*[equation* omitted — 102 chars of source]

where $\alpha_{i} \sim N (\nu_{\alpha}, 1)$. The two regressors are generated as $x_{it,j} =\,\ N(0,1) + \alpha_i$, for $j= 1,2$, and $\beta_1=\beta_2=1$. The error term $u_{it}$ follows an i.i.d standard logistic distribution. We study panels of $N=(100,200)$ individuals observed for $T=(8,16)$ time occasions. We control for the degree of CS by setting $\nu_{\alpha} = 1,2$, with a proportion of subjects without individual outcome variation ranging from 40% to 80%. For each scenario, we run 1000 Monte Carlo simulations.

The number of groups $K$ chosen for the GFE approach is implied by a set of equally spaced values of the hyperparameter $\gamma = (0.1, 0.4, 0.7, 1)$. Larger values of $\gamma$ correspond to fewer groups, with $\gamma=1$ yielding the smallest $K$ and therefore the strongest reduction in the number of units incurring in CS. Each value of $\gamma$ implies that $K$ is within $\sqrt{T}$ and $T\sqrt{N}$ in each scenario considered.

We compare the plug-in GFE APE estimator for $x_1$ with an infeasible APE estimator that computes average effects only for subjects with outcome variation over time. We also compare the performance of the GFE approach with four alternative APE estimators: the FE plug-in ML estimator, the analytical and iterated jackknife bias corrected APE estimators by HN2004, and the APE estimator that plugs in ML Firth-regularized estimates.

Tables (ref)-(ref) report the mean and median ratios between the estimated and the real population APE, the APE standard deviation ("S.D."), and the empirical size of a two-sided $t$-test\footnote{We use analytical standard errors obtained via Delta Method.} centered in the population APE at significance levels 0.05 and 0.1 ("p .05" and "p .10"). We also report the percentage of observations removed due to CS and the average number of groups (K) implied by the chosen values of $\gamma$.

First of all, it is worth noting that the infeasible estimator systematically presents a ratio much greater than one, clearly showing that removing observations in CS unavoidably leads to an overestimation of the APE. Coherently, this bias decreases in the percentage of subjects without variability in the response configuration (denoted by the % of CS for the ML estimator), which gets smaller as $T$ increases and $\nu_\alpha$ is set to $1$.

The plug-in ML estimator of the APE does not apparently exhibit an upward bias, as it is likely to be offset by the IPP one, which can still shift the sampling distribution when $T$ is small dj2015. Nevertheless, for this estimator, coverage issues arise when the percentage of units in CS is elevated. The upward bias in the APE estimator shows up as soon as the IPP bias is reduced by either an analytical or jackknife correction, thus also affecting coverage accuracy. Finally, the APE estimator that plugs in ML Firth-regularized estimates shows an unsatisfactorily finite sample performance.\footnote{We find that, in the scenarios considered, the FE regularized estimates are systematically smaller than the true individual intercepts, which clearly leads to larger estimated individual partial effects and to an upward bias in the APE estimator.}

commentNo alternative estimator has comparable performance: the FE estimator has an overall limited bias, but this evidence is due to the IPP bias which has a downward sign; Moreover, the coverage is not precise. As discussed in section (ref), the FE and the debiased estimators tend to overestimate the APE and the quantification of it is rather inaccurate, especially when $T$ is short and the CS is strong. In addition to this, the empirical size is rather distant from nominal levels, making the inference not reliable.

The regularization entailed by the GFE approach effectively reduces the instances of complete separation for all the values of $\gamma$, and thus the number of groups considered. Regarding its finite sample performance, overall the mean and median ratios display smaller biases with respect to the alternative estimators considered, and the larger number of observations retained help to improve the finite-sample coverage.

The performance of the GFE estimator sensitively varies with the number of groups considered in the classification step. In fact, the bias of the ratio increases with the value of the hyper-parameter: this is a result of the number of groups yielded by $\gamma$ not being large enough to provide an adequate approximation of the underlying UH distribution, even though the average $K$ across simulations complies with the guidelines to choose the number of groups, i.e, $K > \sqrt{T}$. This is expected in our design, as the UH is normally distributed and its support is, for instance, approximated only by roughly $6$ to $8$ points when $\gamma = 1$. For this reason, it is advisable to choose a value of $\gamma$ which implies $K \gg \sqrt{T}$.

An increase in bias should also be expected for very small values of the hyper-parameter, as a larger number of groups operating the discretization could give rise to an IPP bias in finite samples. However, this issue does not arise in the scenarios considered as, with $\gamma$ as small as $0.1$, the implied number of GFEs to estimate does not seem to be large enough for such bias to show up prominently.

Finally, it should be noted that the finite-sample coverage of the GFE estimator does not improve with larger sample sizes. This is likely due to the fact that, on average, the number of groups increases only slightly or remains stable when $N$ doubles, so that the approximation bias does not decay, while the confidence interval shrinks instead.

commentThe proposed approach significantly reduces the number of CS instances, the mitigation increasing with the $\gamma$ parameter. For instance, in Table (ref) the percentage of observations dropped by the GFE estimator is half the ML one for the highest value of $\gamma$. In turn, the increased number of observations reverberates in the finite sample performance of the plug-in GFE APE, which exhibits low bias. At the same time, the inference is more precise, as is clear from values of the empirical size reported in the Tables, which are close to the nominal size and almost always in the confidence interval.
table[table omitted — 3,350 chars of source]
table[table omitted — 3,355 chars of source]

Dynamic logit model

We also study the finite sample performance of the GFE approach by estimating a dynamic logit model in presence of CS. For $i =1,\dots, N$ and $t=1,\dots,T$, we generate the outcome variable as

equation*[equation* omitted — 123 chars of source]

where $\theta_1=\theta_2=1$ and $\beta=0.5$. The two regressors and the time-invariant FE are generated as $ x_{it,j} =\,\ N(0,1) + \alpha_i$ for $j=1,2$, with $\alpha_i \sim N(\nu_{\alpha},1)$, respectively.

We study panels of $N=(100, 200)$ individuals observed for $T=(8,16)$ time occasions. We control for the degree of CS with two values of $\nu_{\alpha}=(0,-1)$, which results in a percentage of 24% to 50% of units without outcome variation. We run 1000 Monte Carlo simulations for each scenario. We report simulation statistics for the APE estimator of $y_{t-1}$ in Tables (ref)-(ref). The estimators analyzed and the values of $\gamma$ are the same selected for static design, with two exceptions: (i) we do not include the APE estimator that plugs in ML Firth-regularized estimates, since its employment in dynamic settings lacks a theoretical background, and (ii) we use bias-correction methods suited for dynamic models, namely the analytical one of fernandez2009fixed and the half-panel jackknife estimator dj2015.

As in the static case, the infeasible estimator systematically overestimates the population APE with the bias decreasing in $T$ and in the values considered for $\nu_{\alpha}$. The FE plug-in ML estimator and both the analytical and jackknife bias corrected APE estimators exhibit poor performance. When the dimension $T$ is short, the IPP is severe and, in turn, the overestimation of the APE is offset by a strong downward bias. Bias corrections manage to improve the mean and median ratios as $T$ increases, although the coverage remains overall inaccurate.

The GFE approach preserves its regularizing properties in the dynamic setting and manages to reduce the instances of CS for all values of $\gamma$. However, the finite sample properties of the GFE plug-in estimator suggest that regularization does not fully offset the stronger IPP bias that arises in dynamic settings, causing the APE to be systematically underestimated when $\gamma = 0.1$. Accordingly, this bias decreases with larger $T$. However, the bias already decreases and the empirical coverage attains its nominal values with intermediate values of $\gamma$, such as $0.4$ and $0.7$, especially when $T=16$.

Appendix (ref) contains additional simulation evidence related to a data generating process that violates the assumption of the stationarity of regressors blm2022. Table (ref) shows, however, that the results on the finite sample properties of the proposed estimator are robust to the inclusion of a trending regressor.

table[table omitted — 3,030 chars of source]
table[table omitted — 3,034 chars of source]

Empirical applications

Female labor force participation

We revisit the empirical application on inter-temporal labor supply decisions of women, also illustrated in dj2015. Data are related to the employment status of $N=1461$ married women aged between 18 and 60 in 1985, whose husbands were always employed in the period from 1981 to 1988, observed for $T=8$ years (PSID waves 15-22). We estimate a dynamic logit model and include control variables such as the number of kids of different ages, the logarithm of the yearly income of the husband, the age, and the age squared.

The employment status exhibits strong inter-temporal correlation: 143 women are unemployed for the whole period, while 719 women are always employed. Therefore 862 units, that is around $60\%$ of the sample, do not exhibit any outcome variation and are dropped due to CS when a FE model is estimated. We should expect the GFE to keep increasingly more units as we increase the value of the hyper-parameter $\gamma$.

In Table (ref) we compare the results obtained by the GFE approach with four alternative APE estimators: the plug-in pooled estimator, the FE plug-in ML estimator, the analytical bias corrected APE estimators by fernandez2009fixed and the half-panel jackknife APE estimator by dj2015. For what concerns GFE, we report estimates for $\gamma=0.4,0.6,1$ and the GFE APE estimates when $K$ is fixed and equal to 5.

The APEs for every variable across each estimation method are in line with the corresponding economic intuition, but the magnitude of the the APE for $y_{i,t-1}$ is rather different across estimators. The pooled estimator indicates a strong positive effect, which likely reflects an upward omitted variable bias from ignoring UH. In contrast, the FE ML estimator yields a much lower APE, which may indicate a downward bias due to the IPP. Consequently, the analytical and jackknife bias-corrected estimators mitigate this issue, giving an estimated APE for $y_{i,t-1}$ twice the one obtained by the FE ML estimator.

In order to select the proper value of the GFE hyperparameter, we follow the proposed rule that implies $K <\sqrt{N}T \approx 305$ and $K\gg \sqrt{T}\approx 3$. Out of the values of $\gamma$ giving rise to the estimates in Table (ref), only $\gamma=0.6,1$ are compliant with the rule, while $\gamma=0.4$ violates it. Also notice that, while greater than $\sqrt{T}$, $K=5$ is too close to the lower bound, as the results are identical to the pooled model, thus indicating that the number of approximating points is too small to guarantee a good description of the UH. When $\gamma=0.4$, the approximation of the UH is likely to be sufficient but the estimated number of grouped FE is too large to control the IPP bias. In this vein, the choice of $\gamma=0.6$, where the number of parameters is 1/4 with respect to FE estimation and the approximation of the UH is delivered by 233 support points, is suggested in this case.

As summarized by Figure (ref), which reports the plug-in GFE APE estimator of $y_{i,t-1}$ for 20 values of $\gamma$, the GFE approach always gives a quantification of the effect that is greater than ML alternatives, in line with the findings in the simulation study. Moreover, for increasing values of the hyperparameter, the plug-in GFE APE moves towards the plug-in pooled estimator, although in this case the number of groups with $\gamma=1$ is still sizable and equal to 117.

Regarding the GFE, the percentage of dropped observations is decreasing in $\gamma$ and the proposed approach stops dropping units for values of the hyperparameter higher than $0.8$. Figure (ref) shows the decreasing trend of discarded units for increasing values of $\gamma$.

table[table omitted — 2,286 chars of source]
figure[figure omitted — 739 chars of source]
figure[figure omitted — 518 chars of source]

An early warning system for banking crises

An early warning system for banking crises is a binary choice model where the outcome variable takes value 1 if a banking crisis occurs in country $i$ at time $t$ and 0 in non-critical periods. The probability of a crisis is modeled as a function of lagged macroeconomic and financial indicators that are supposed to warn about the likelihood of a crisis in advance. The dataset in exam, issued by laeven2018systemic, consists of a balanced panel of $N=33$ countries observed over the years 1986 - 2015 ($T=30$). laeven2018systemic give the definition of a banking crisis for a large set of countries and identify 69 crisis episodes over 990 data points, so we have a panel dataset where the dependent variable is an extremely rare event. In addition to the one period lagged dependent variable $y_{i,t-1}$, macroeconomic variables used in the analysis, available as International Financial Statistics (International Monetary Fund) or World Development Indicators (World Bank) are: real GDP growth, the log of per capita GDP, inflation, real interest rate, the ratio of M2 (broad money) to foreign exchange reserves, the growth rate of real domestic credit and the growth rate of foreign assets. All explanatory variables are lagged by one period and further description of the dataset can be found in pigini2021penalized and caggiano2016comparing.

CS instances are a major issue in forecasting crises: in fact, the FE logit model cannot be used to predict the occurrence of a crisis for countries that never experienced one, as the estimates of the FEs would not be finite. In order to illustrate how the proposed approach can circumvent this problem, Figure (ref) depicts the empirical density of in-sample predicted probability of crisis for both ML and the GFE estimator ($\gamma=0.5$). The ML estimator drops 13 countries out of 33 due to CS and, as a result, we observe a large probability mass in 0. In contrast, the GFE approach drops only one country so that the empirical density of the forecast probability turns out to be right-shifted compared to the ML one, thus allowing non-trivial predictions of crises events for countries without outcome variation.

figure[figure omitted — 382 chars of source]

We also perform a one-step-ahead forecast exercise. Using an expanding training set stopping at years 2006 to 2010, we estimate a dynamic logit model and compute the out-of-sample predicted probabilities using the next year. The last forecast year is 2011, as every year after that does not present any crisis in the dataset. The cut-offs used to compile confusion matrices are chosen by optimizing the in-sample sum of specificity and sensitivity. We compare the forecasting performance of the GFE approach to that of the FE ML estimator and analytical bias-corrected ML estimator fernandez2009fixed. We experiment with 4 values of $\gamma=(0.005,0.1,0.5,1)$.

Figure (ref) reports out-of-sample F1 score for ML and GFE with $\gamma=0.5$ for all forecast years: the latter strictly outperforms the former, achieving perfect classification in two out of five scenarios (2009 and 2011). The better F1 score for GFE is strictly due to the higher rate of false negatives detected. For the sake of clarity, it is interesting to note that the number of groups found by the GFE procedure with $\gamma=0.5$ in the first step varies in time over the training sets - ranging from 7 to 9.

figure[figure omitted — 304 chars of source]

The complete set of results of this forecast exercise is reported in Table (ref) in Appendix (ref).

comment, where we report customary statistics such as the number of true positives, true negatives, false positives, false negatives and the F1 score together with the percentage of countries in CS and the number of groups for GFE. Both standard ML and GFE estimators manage to forecast a good amount of crises - 5 out of 7 in the 5 years - although the GFE exhibits a better performance in correctly predicting noncritical events. Further considerations on these results should take into account the cost of an undiagnosed crisis with respect to that of a false alarm, which is something we prefer not to take a stand on, since it is heavily dependent on policy-makers' priors.

Overall, the forecast performance of at least one of the GFE estimators considered in the exercise is always better than the one of the ML or BC estimators. In order to provide guidance on the choice of the hyperparameter, we suggest the choice of $\gamma=0.5$ as it is the value that gives the best performance in compliance with the proposed rule, since $\gamma=1$, which would be slightly better in terms of F1 score, violates the lower bound.

Conclusions

This paper motivates the use of the recently developed GFE approach to perform regularized estimation of binary choice FE models in presence of severe CS. In such settings, FE models exhibit several deficiencies, including biased estimates of APE, inaccurate coverage, and the inability to generate meaningful predictions for units affected by CS.

We provide a simulation study concerning both static and dynamic specifications of logit models. Our results show that, by estimating a smaller number of FE, the proposed approach reduces the instances of CS and yields unbiased APE estimates with improved coverage properties, relative to the available alternatives. Moreover, by keeping all units grouped in clusters with response variability, the GFE approach enables predictions for a much larger number of subjects in the sample.

We also provide two illustrative examples, namely an analysis of determinants of labor force participation and a logit-based early warning system for rare bank crises. The first one shows that, by tuning the hyperparameter so as to provide a trade-off between a good approximation of the UH and a limited number of FEs to estimate, the GFE quantification of the APE for the lagged dependent variable diverges from the ML and bias corrected ones, while mitigating the potential omitted variable bias possibly exhibited by the pooled model. The second example focuses on forecasting and illustrates how the GFE approach manages to offer predictions for units that never experience a financial crisis in the training set.

Acknowledgments

We are grateful to Pavel \v{C}i\v{z}ek, Fulvio Corsi, Riccardo Lucchetti, Silvia Sarpietro, Laura Serlenga, Amrei Stammann, Rainer Winkelmann, to the participants of : the 12th Workshop of Econometrics and Empirical Economics, 13th IAAE annual conference, 30th IPDC, seminar at the University of Pisa, for their helpful comments and suggestions. Claudia Pigini and Alessandro Pionati would like to acknowledge the financial support by the European Union - Next Generation EU - Prin 2022, Project Code: 2022TZEXKF 02; Project CUP: I53D23002800006; Project Title: “Hidden Markov Models for Early Warning Systems” and financial coverage D.D. MUR 47/2025, CUP I33C25000280001; Project Title: “The use of the Grouped Fixed Effects estimator in panel data analysis addressing unobserved heterogeneity”.

footnotesize