A Simple and Computationally Trivial Estimator for Grouped Fixed Effects Models

A Simple and Computationally Trivial Estimator for Grouped Fixed Effects Models

\thispagestyle{empty}

\thispagestyle{empty}

\pagenumbering{arabic}

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Introduction} Suppose a sample of panel data $\{(y_{it},x_{it}):1\leq i\leq N, 1\leq t\leq T\}$ is observed and consider a linear regression model with grouped fixed effects:

where $i$ denotes cross-sectional units, $t$ denotes time periods, $y_{it} \in \mathbb R$ is a dependent variable, and $x_{it}\in\mathbb R^K$ is a vector of explanatory covariates uncorrelated with the zero-mean random variable $v_{it}\in\mathbb R$ but possibly arbitrarily correlated with the unobserved group membership variable $g_i\in\{1,\ldots,G\}$ and the group-time effect $\alpha_{g_it}\in\mathbb R$.

This paper focuses on the estimation of and inference on the unknown slope parameter $\beta\in\mathbb R^K$ and the group-time effects $(\alpha_{1t},\ldots,\alpha_{Gt})'\in\mathbb R^{G}$, as well as the consistent estimation of the group memberships $g_i\in\left\{1,\ldots,G\right\}$ and the number of groups $G$, within an asymptotic framework in which $N/T^\nu\to0$ for some constant $\nu>0$, as $N$ and $T$ diverge while $K$ and $G$ remain fixed.

Since they are special cases of interactive fixed effects with factor loadings confined to a finite set,\footnote{Note that $\alpha_{g_it}=\lambda_i'f_t$ for any $\lambda_i'\equiv\left(c_1\mathbf{1}\{g_i=1\},\ldots,c_G\mathbf{1}\{g_i=G\}\right)$, $f_t'\equiv\left(\alpha_{1t}/c_1,\ldots,\alpha_{Gt}/c_G\right)$, and $c\equiv(c_1,\ldots,c_G)'\in(\mathbb R\backslash\{0\})^{G}$. Thus $\lambda_1/c,\ldots, \lambda_N/c$ lie in the finite set of vertices of the unit simplex of $\mathbb R^G$. Reciprocally, if $\tilde f_t\in\mathbb R^r$ for some $r\in\mathbb N\backslash\{0\}$ and $\tilde \lambda_i\in\Lambda\subset\mathbb R^r$ with $\left\vert\Lambda\right\vert=\tilde G$, then there exist $(\tilde g_1,\ldots,\tilde g_N)'\in\{1,\ldots, \tilde G\}^N$ and $(\tilde\alpha_{11},\ldots,\tilde\alpha_{\tilde GT})'\in\mathbb R^{\tilde GT}$ such that $\tilde \lambda_i'\tilde f_t=\tilde\alpha_{\tilde g_it}$. } grouped fixed effects (GFE hereafter) provide a parsimonious yet flexible device to accommodate cross-sectional correlations and a few unrestricted trends of unobserved heterogeneity. Since their introduction in economics HahnMoon2010, BM2015, GFE have gained considerable interest in both methodological and applied work SSP2016, cheng2019clustering, BLM2019,GU2019, ChetverikovManresa2021, BonhommeLamadonManresa2022, mugnier_JMP,Janys2024.

Treating GFE as interactive fixed effects, however, leads to two main issues. First, parametric-rate inference for the slope parameter is generally not available when $T$ grows very slowly with $N$ Bai2009, moon2019nuclear, BeyhumGautier2023. Second, although Higgins2022 shows that parametric-rate inference remains possible under some circumstances, his proposed method, like most other interactive fixed effects methods, solves a high-dimensional non-convex least-squares problem typically subject to local minima. Both issues can lead to poor inferences, especially in microeconometric datasets where $T$ may be much smaller than $N$.

Addressing the combinatorial nature of GFE, on the other hand, introduces two main difficulties. First, GFE estimators, defined as pseudo joint maximum likelihood estimators optimizing over all possible partitions of cross-sectional units in $G$ groups, encounter a challenging non-convex and combinatorial optimization problem. If P$\neq$NP, exact solutions in polynomial time are unattainable for most real-world datasets of interest.\footnote{GFE estimators are instances of minimum sum-of-squares clustering, which is NP-hard unless both $G$ and $T$ are fixed InabaKatohImai1994, Aloise2009NPhardnessOE,MAHAJAN201213.} Second, existing GFE estimators typically require the number of groups, or an upper bound for it, to be known to the analyst. While BIC-type criteria are known to achieve consistent model selection when $N$ and $T$ grow at the same rate, to the best of my knowledge, no such result is available in the asymptotic regimes considered here.

This paper provides a three-step estimation procedure free of these limitations. The first step leverages the low-rank factor structure of the linear GFE model to offer a computationally simple preliminary estimator of the slope coefficient. Though some emphasis is put on smooth and convex regularized nuclear-norm estimation, the interactive fixed effects literature provides several candidates further discussed. The second step examines residual correlations within triads of units to construct distances between units that, for each pair of units, are shown to be asymptotically zero if and only if both units belong to the same group, whenever group-time effects are well separated. An agglomerative clustering based on distance thresholding, in turn, is a natural choice. The third step computes an ordinary least squares regression (OLS) that controls for interactions of time and estimated group dummies.

While the novel procedure combines disparate ideas from the existing panel data literature, the combination itself, along with its appealing large-sample statistical properties, appears to be new. Specifically, consistent model selection -- i.e., consistent estimation of the number of groups -- is achieved jointly with the estimation of other parameters by merging similar units into groups at a rate governed by time-dependence conditions. An advantage over traditional BIC-type selection methods is that it does not require estimating an arbitrary range of possibly misspecified models, determined by an upper bound $G_{\max}\geq G$, and applying some selection rule. The latter is often computationally cumbersome for GFE estimators. Since $G$ is unknown, however, regularization is necessary, and there is no free lunch. A computationally simple, data-driven procedure for determining the minimal tolerated distance to merge units is proven to be theoretically valid and demonstrates reasonable performance in finite samples.

Unlike standard GFE estimators, the proposed method is both computationally “trivial” and exact, as the combinatorial optimization problem of GFE -- often solved using heuristics that provide local solutions -- is replaced by a sequential agglomerative clustering procedure with polynomial-time complexity. Fast implementations of agglomerative clustering procedures are already available in standard software used by economists (e.g., R or MATLAB). In addition, the preliminary estimation step can be performed via convex optimization, for which computationally efficient algorithms also exist.\footnote{A full implementation in MATLAB is available at \href{https://github.com/martinmugnier/TPWD-Estimators}{https://github.com/martinmugnier/TPWD-Estimators}.}

The main theoretical contribution of this paper is to provide standard regularity conditions under which the first two steps yield consistent estimation of well-separated groups. As in similar methods, this leads to an estimator of the slope parameter and the group-fixed effects asymptotically equivalent to the infeasible regression controlling for the true group memberships. Compared to existing approaches, the proposed method is computationally simple, endogenously determines the number of groups using a theoretically valid data-driven selection rule, and relies on minimal assumptions about the covariates -- requiring only sufficient conditional variation through a restricted eigenvalue condition. This generalization rules out cross-sectional correlation in idiosyncratic errors, as the new estimator uses triad-specific comparisons rather than cross-sectional averages to recover the group structure.

While the main text focuses on a simple linear model with a homogeneous slope to convey the core proof ideas, the versatility of the new approach is illustrated through several extensions with a growing number of groups, heterogeneous slope coefficients (in time or across units), and multiplicative network models (e.g., gravity trade equations), all discussed in Appendix Section (ref). This paper does not address inference on group memberships. Dzemski2024 develop pointwise valid inference methods in a model with group-specific slopes, but ruling out time-varying group-specific effects, provided a preliminary estimator of the slopes is available. If this approach could be adapted to accommodate a homogeneous slope and time-varying effects, the estimator proposed in the present paper could plausibly be used.

The finite-sample performance of the new estimator is compared to that of state-of-the-art alternatives across various Monte Carlo experiments calibrated to the empirical application: potentially misspecified GFE, spectral, and post-spectral estimators with different pre-specified group numbers as proposed in BM2015 and ChetverikovManresa2021, nuclear-norm regularized (NNR) and nuclear-norm (NN) estimators from moon2019nuclear, and the interactive fixed effects (IFE) estimator studied in Bai2009, initialized with either the NNR estimator or random draws. Without covariates, the new method outperforms all estimators except the well-specified GFE, closely matching its performance for moderate values of \(T\). It achieves remarkably homogeneous clustering in terms of Precision, Recall, and Rand Index, which are defined later. With a scalar covariate, the new estimator significantly improves bias, root mean square error, and coverage compared to all alternatives, except the well-specified GFE, GFE with a well-specified BIC criterion (\(G_{\max} \geq G\)), and the well-specified IFE in terms of bias. For \(T \geq 20\), its performance approaches that of the infeasible pooled OLS regression, with confidence intervals based on a consistent asymptotic variance estimator aligning closely with the nominal 95% level.

The usefulness of the method is demonstrated by revisiting BM2015’s analysis of the statistical relationship between income and democracy in a large panel of countries from 1970 to 2000. This association may be confounded by critical junctures in history that led to similar unobserved development paths for some countries and different ones for others. While the authors find statistically significant approximations of the GFE estimates for the effect of lagged log-income per capita on a democracy index ranging from 0.061 to 0.089, depending on the pre-specified number of groups, and suggest fewer than 10 groups, the proposed method identifies 4 groups, an effect of 0.07, and a larger cumulative income effect of 0.258 (vs. from 0.104 to 0.151). The preliminary regularized nuclear-norm estimator delivers point estimates of 0.016 and 0.078, respectively.

\paragraph{Related Literature.}

This paper contributes to the extensive literature on estimating panel data models with interactive fixed effects Pesaran2006, Bai2009, MoonWeidner2015, MoonWeidner2017, BonhommeLamadonManresa2022, armstrong2022robust, BeyhumGautier2023. While convex nuclear-norm regularized estimators moon2019nuclear, chernozhukov2019inference are computationally simpler than non-convex least squares estimators Bai2009, they may converge at slower-than-parametric rates. The first contribution of this paper is to leverage off-the-shelves bounds on the rate of convergence for such estimators, as derived in this literature, to construct computationally simple estimators of the special case of grouped fixed effects models that converge at the parametric rate. ChetverikovManresa2021 independently employed similar ideas, though they impose some factor structure on the covariates and a known upper bound on the number of groups in order to apply spectral clustering techniques. Differently, the method proposed in this paper relies on a plausibly weaker restricted eigenvalue condition and achieves model selection simultaneously. Using the nuclear-norm regularized estimator in the first step, it requires two regularization hyperparameters, which serve as data-driven substitutes for a known upper bound $G_{\max}\geq G$ and the BIC (or AIC) model selection criteria often used to select the number of factors in interactive fixed effects models Bai2003. These hyperparameters, however, depend on fundamental sampling properties of the data, such as cross-sectional and time dependence, rather than prior information about the maximal dimension of the model.\footnote{Alternatively, if a tuning-free nuclear-norm estimator is used in the first step, such as moon2019nuclear's nuclear norm estimator, the method requires only one hyperparameter for the clustering step. While the available assumptions ensuring a sufficiently fast rate of convergence for the nuclear-norm estimator are relatively high-level compared to those for the regularized nuclear-norm estimator, both methods exhibit nearly identical finite-sample properties in the Monte Carlo simulations reported in Section (ref).}

This paper also contributes to the rapidly growing literature on estimating grouped fixed effects models. BM2015's GFE estimator, an extension of $k$-means clustering to handle covariates, solves an NP-hard optimization problem. Algorithms that provide fast solutions may fail to converge to the true estimate values defined as a global minimum, and the same drawback applies to extensions and other non-convex estimators SSP2016, Ando_Bai_2022, mugnier_JMP,Lumsdaine2023. In contrast, the inferential theory developed in this paper is valid for a computationally simple estimator that replaces a known upper bound on the number of groups with the willingness to merge units based on estimated pairwise distances. Since inference is conducted on a true population parameter, these results contrast with those of Pollard_1981, Pollard_1982, who provide asymptotic theory for the solution to the population $k$-means sum-of-squares problem in the cross-sectional case, which corresponds to a pseudo-true value. Similarly, Lewis_et_al2022 introduced a fuzzy clustering procedure that performs well in simulations, but whose large sample properties, however, have not been derived in the large \(N,T\) framework and thus apply to a pseudo-true value. Its implementation requires the number of groups to be specified by the analyst.

The closest approach is that of ChetverikovManresa2021's spectral and post-spectral estimators. While computationally straightforward, these estimators require prior knowledge about the number of groups or a consistent estimator (not provided by the authors) and the theoretical validity of the spectral (resp. post-spectral) estimator crucially rests upon a factor (resp. grouped) structure for the covariates. Such assumptions bring the model closer to random or correlated random effects in the spirit of Pesaran2006 and could be restrictive in practice. Differently, the method proposed in this paper does not impose a specific model for the covariates or a known upper bound on the number of groups. Recently, Mehrabani2023 adapts “sum-of-norm” convex clustering Hocking2011, Tan2015 to a linear panel data model with latent group structure but time-constant group effects. Extending this approach to accommodate time-varying effects seems challenging. Albeit close in spirit, the proposed procedure differs from the binary segmentation algorithm developed in KeLiZhang2016 and WANG2021272. Another approach involves applying spectral clustering to some dissimilarity matrix Ng_2002, vonluxburg2007tutorial, ChetverikovManresa2021, brownlee_Lugosi2022, LuGuVolgushev_2022. One limitation is to introduce additional complexity and tuning parameters, as \(L\) eigenvectors of the dissimilarity matrix need to be computed and clustered (and \(L\) must be chosen), typically by approximating an NP-hard $k$-means solution, which this paper aims to avoid in order to ensure valid inference.

Finally, this paper establishes a connection between the mature statistical and operation research literature on clustering problems and the recent grouped fixed-effects literature. Although agglomerative or hierarchical clustering methods are well established in the former, their adaptation to the latter is relatively new Vogt2016,Chen2019, MammenWilke2022. The advantage lies both in computational efficiency and valid inference, breaking the high-dimensionality of $k$-means by considering agglomerative approaches and leveraging the discrete structure of the econometric model. Hence, this paper can be seen as the first application and analysis of an agglomerative clustering method in the econometric panel data model (ref).\footnote{Since the first arXiv version of this paper, FreemanWeidner2023 also propose a hierarchical clustering algorithm but do not explicitly verify that it meets the approximation conditions for their large sample theory to be valid. Differently from $k$-means BonhommeLamadonManresa2022,GrafLuschgy2002, it remains unclear when these conditions are met.} The pairwise distance used in the clustering step has already been employed in the mathematical statistics and econometric literature to study topological properties of the graphon Lovasz2012,ZhangLevinaZhu2017, Zeleneev2020, Auerbach2022. More generally, dyad, triad, or tetrad comparisons have proven useful in a variety of other econometric contexts HONORE1994241,Graham2017, Charbonneau2017, jochmans_2017.

The rest of the paper is organized as follows. Section (ref) outlines the three-step estimation procedure. Section (ref) discusses the large-sample properties, including uniform consistency for the grouping structure and asymptotic normality at parametric rates. Section (ref) presents the results of a small-scale Monte Carlo simulation. Section (ref) reports the results of the empirical application. Section (ref) concludes. All proofs are provided in the Appendix. Additional material can be found in the Supplementary Material mugnier2022SM, with section numbers denoted as S.1, etc.

\paragraph{Notation.} Let $\mathbf{1}\{ \cdot \}$ denote the indicator function. For any set $S\subset \mathbb R^p$ (for any $p\geq 1$), let $S^*\equiv S\backslash\{0\}$ and $\vert S\vert $ denote the cardinality of $S$. For any set $I\subset\mathbb N$, all $k\in\mathbb N^*$ such that $k\leq\vert I\vert$, let $\mathscr P_k(I)$ denote the set of subsets of $I$ with cardinality $k$. The operators $\overset{p}{\to}$, $\overset{d}{\to}$, and $\mathrm{plim}$ denote convergence in probability, convergence in distribution, and probability limit, respectively. All vectors are column vectors. $\Vert\cdot\Vert$ denotes the Euclidean norm. For an $m\times n$ real matrix $A$ of rank $r$, let $A'$ denote its transpose, $\Vert A\Vert_{\textrm{F}}\equiv[{\rm Tr}(A'A)]^{1/2}$ its Frobenius norm, $\Vert A\Vert_\infty\equiv\sigma_1(A)$ its spectral norm, $\Vert A\Vert_{\max}\equiv\max_{i=1,\ldots,m;j=1,\ldots,n}\left\verta_{ij}\right\vert$ its max norm, and $\Vert A \Vert_1\equiv\sum_{i=1}^r\sigma_i(A)$ its nuclear norm, where $\sigma_1(A)\geq\cdots\geq\sigma_r(A)>0$ denote the positive singular values of $A$.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Three-step estimation}

In this section, I introduce the three-step triad pairwise-differencing (TPWD) estimator $\widehat\theta$ for the parameter $\theta\equiv(G, g_1,\ldots, g_N,\alpha_{11},\ldots,\alpha_{GT}, \beta')'\in\Theta$ of Model (ref), where $$\Theta\equiv\bigcup_{g\in\mathbb N^*}\Theta_{g} \quad \text{with} \quad \Theta_{g}\equiv\left\{g\right\}\times\left\{1,\ldots,g\right\}^N\times\mathcal A^{gT}\times\mathcal B,$$ for some subsets $\mathcal B\subset\mathbb R^K$ and $\mathcal A\subset\mathbb R$. The dependence of $\Theta$ on $N,T$ is omitted. Sections (ref)--(ref) describe the three key components of $\widehat \theta$: a preliminary consistent estimator of $\beta$, a measure of pairwise distances between units, and an agglomerative clustering algorithm. Section (ref) formally defines $\widehat \theta$.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Preliminary consistent estimation of the slope coefficient}

The first component of the three-step TPWD estimator is a preliminary consistent estimator $\widehat \beta^1$ of $\beta$:

Intuitively, since covariates can be arbitrarily correlated with the unobserved grouped fixed effects, the clustering problem is simpler in the approximate “pure” grouped fixed effects model: $$y_{it}-x_{it}' \widehat\beta^1= \alpha_{g_it}+v_{it}+o_p(1), \quad i=1,\ldots,N, \; t=1,\ldots,T,$$ as $\min(N,T)\to\infty$.

The interactive fixed effects literature offers several computationally simple estimators that verify (ref) under various identifying assumptions moon2019nuclear, chernozhukov2019inference, beyhum2019squareroot, BeyhumGautier2023. The large-sample properties of \( \widehat{\theta} \), established in Section (ref), do not depend on the convergence rate of the preliminary estimator, which is typically slower than \( \sqrt{NT} \), as long as a lower bound on this rate is known and used to calibrate the clustering tuning parameter. Such bounds are available for all the references mentioned above, provided certain regularity conditions hold.

As an illustrative example, consider below the nuclear-norm regularized estimator. Under a known upper bound on the number of groups and a factor structure on the covariates, another example is the correlated random effects estimators proposed by ChetverikovManresa2021. Simulation results reported in Section (ref), however, suggest that while ChetverikovManresa2021's spectral norm estimator may outperform the nuclear-norm regularized estimator under a well-specified covariate model, the iterated version of the triad pairwise-differencing estimator based on the nuclear-norm regularized estimator may outperform ChetverikovManresa2021's spectral and post-spectral estimators for small values of $T$.\footnote{For large values of \(T\), which is not the typical setting for which GFE models were developed, both estimators exhibit equivalent performance (although spectral and post-spectral estimators require the number of groups as input); in this case, however, many other estimators compete.}

\paragraph{Nuclear-norm regularization.}

Let $Y\equiv(y_{it})_{i=1,\ldots,N;t=1,\ldots,T}\in\mathbb R^{N\times T}$ and $X_k\equiv(x_{it,k})_{i=1,\ldots,N;t=1,\ldots,T} \in \mathbb R^{N\times T}$ for all $k\in\left\{1,\ldots,K\right\}$. For all $v=(v_1,\ldots,v_K)'\in\mathbb R^K$, let $v\cdot X\equiv\sum_{k=1}^KX_kv_k$. For any value of the regularization parameter $\psi_{NT}\in(0,\infty)$, let $Q_{\psi_{NT}}$ denote the nuclear-norm regularized concentrated objective function:

for all $\beta\in\mathbb R^K$. A nuclear-norm regularized (NNR) estimator is defined as a solution to the following minimization problem:

Under regularity conditions, $\widehat\beta^1(\psi_{NT})$ is unique with probability approaching one moon2019nuclear.

Instead of setting $\psi_{NT}=0$ and optimizing over all $G^{\rm guess}$ ($G^{\rm guess}<<N$) unobserved grouped trend patterns $$\Gamma\in\left\{\Lambda F': \Lambda\in \left\{0,1\right\}^{N\times G^{\rm guess}},F\in\mathbb R^{T\times G^{\rm guess}},\sum_{j=1}^{G^{\rm guess}}\Lambda_{ij}=1\quad i=1,\ldots,N\right\}\ldots,N},$$ which leads to the NP-hard problem solved by BM2015's GFE estimator, the NNR estimator penalizes the nuclear norm (the sum of singular values) of an unrestricted matrix $\Gamma\in\mathbb R^{N\times T}$ of individual trends. This convexifies the rank-constrained problem solved by Bai2009's interactive fixed effects (IFE) estimator -- which sets $\psi_{NT}=0$ and optimizes over $\Gamma\in\{\Lambda F': \Lambda\in\mathbb R^{N\times G^{\rm guess}}, F\in\mathbb R^{T\times G^{\rm guess}} \}$ -- thus solving the “local minima” problem.\footnote{The nuclear norm $\Vert \Gamma\Vert_1$ is the convex envelope of rank$(\Gamma)$ over the set of matrices with spectral norm at most one.}

The tuning parameter $\psi_{NT}$ performs model regularization without the need to set both $G^{\rm guess}\geq G$ and a model selection rule. As proved in Section (ref), under regularity assumptions and a rate condition on $\psi_{NT}$, the interactive fixed effects structure of grouped fixed effects is sufficient for $\widehat\beta^1(\psi_{NT})$ to achieve a convergence rate arbitrarily close to $\sqrt{\min(N,T)}$. In the Monte Carlo experiments, the theoretically justified choice $\psi_{NT}\equiv\log(\log(T))/\sqrt{16\min(N,T)}$ is employed.

\paragraph{Computation.} Minimization problem (ref) is convex and can be efficiently solved using modern optimization techniques. Algorithm 1 below follows a straightforward iterative strategy that alternates between the proximal gradient algorithm proposed by Mazumder2010 and a projection step. For any real matrix $A$ of rank $r$ and any scalar $\tau\in[0,\infty)$, let $\mathbf S_\tau(A)\equiv UD_\tau V'$, where $D_\tau\equiv{\rm diag}[(\sigma_1(A)-\tau)_+,\ldots,(\sigma_r(A)-\tau)_+]$, $UDV'$ is the compact singular value decomposition of $A$, $D\equiv{\rm diag}(\sigma_1(A),\ldots, \sigma_r(A))$, and $t_+\equiv \max(t,0)$. For all $\beta\in\mathbb R^K$ and $\Gamma\in\mathbb R^{N\times T}$, let \[ L_{\psi_{NT}}(\beta,\Gamma)\equiv\frac{1}{2NT}\Vert{Y-\beta\cdot X - \Gamma}\Vert^2_{\textrm{F}}+\frac{\psi_{NT}}{\sqrt{NT}}\Vert \Gamma\Vert_1. \]

Algorithm 1 alternates between soft-thresholding the singular values of the residual matrix $Y-\beta^{\rm old}\cdot X$ to obtain $\Gamma^{\rm new}$, and performing a pooled OLS regression of the residual matrix $Y-\Gamma^{\rm new}$ on $X$ to obtain $\beta^{\rm new}$, iterating until convergence. Commonly employed statistical software includes fast routines for both tasks (e.g., R or MATLAB function svd). By convexity and due to the uniqueness of the minimum of the objective function $L_{\psi_{NT}}$, if Algorithm 1 converges, it will converge to $\widehat\beta^1(\psi_{NT})$. In addition, the optimization error can be made arbitrarily small by selecting a sufficiently small value for $\varepsilon$, which justifies abstracting from optimization errors hereafter.

Since Algorithm 1 may be slow in practice, an alternative approach, based on Lemma 1 from moon2019nuclear, uses a built-in optimization tool to minimize the closed-form expression of the convex concentrated objective function \[ Q_{\psi_{NT}}(\beta)=\sum_{r=1}^{\min(N,T)}q_{\psi_{NT}}\left(\sigma_r\left(\frac{Y-\beta\cdot X}{\sqrt{NT}}\right)\right), \] where, for any $s\in[0,\infty)$, \[ q_{\psi_{NT}}(s)=\left\{

\right. \] For example, an implementation using MATLAB's fminsearch is significantly faster and provides estimates that are nearly identical to those obtained from Algorithm 1 with $\varepsilon=10^{-9}$.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Pairwise distance between cross-sectional units}

The second component of the TPWD estimator is a measure of distance between any two units $i$ and $j$, informative of whether $i$ and $j$ belong to the same group. This can be constructed by exploiting the linear structure of the model and the preliminary consistent estimator $\widehat \beta^1$. Let $\widehat v_{it}\equiv y_{it}-x_{it}'\widehat\beta^1$ denote the first-step residual. The empirical distance between $i$ and $j$, denoted $\widehat d_{\infty,1}^2(i,j)$, is given by

Let $\widehat D\equiv (\widehat d_{\infty,1}^2(i,j))_{(i,j)\in\{1,\ldots,N\}^2}$ denote the symmetric dissimilarity matrix that collects all pairwise distances.

The distance $\widehat d_{\infty,1}^2$ is borrowed from the statistical literature on graphon estimation Lovasz2012, zhang2017estimating, Zeleneev2020, Auerbach2022. While alternative distances could be of interest, e.g.,

the distance $\widehat d_{\infty,1}^2$ performs well in simulations and is convenient to prove large sample properties. For example, following the present proof techniques, using $\widehat d_{\infty,2}^2$ would rule out asymptotically negligible groups in order to ensure consistency of the subsequent clustering step, while $\widehat d_2^2$ could be contaminated by heteroskedastic error variances which would prevent consistency. Similarly, the distance could be based on a self-normalized sum to account for heteroskedastic variances across pairs, which might provide some finite-sample improvements. This is not required for deriving asymptotic results under heteroskedastic errors and will not be pursued further here.

Why is the empirical distance $\widehat d^2_{\infty,1}(i,j)$ informative about whether $g_i=g_j$? The intuition is as follows.\footnote{See also p. 14 in Zeleneev2020 in a network setting.} Since $\widehat\beta^1\overset{p}{\to}\beta$ as $\min(N,T)\to\infty$, it holds with probability approaching one that $\widehat v_{it}\approx \alpha_{g_it}+v_{it}$. Under weak time dependence, tails, and cross-sectional independence restrictions on the error terms, it then holds “uniformly” over $i$, $j$, and $k\in \left\{1,\ldots,N\right\}\backslash\left\{i,j\right\}$ that

If $g_i=g_j$, then $\alpha_{g_it}-\alpha_{g_jt}=0$ and

Reciprocally, if

then necessarily $g_i=g_j$. To see this, suppose that $g_i\neq g_j$. Then, under the weak condition that each group has at least two units asymptotically, there exist $k^*,l^*\in\left\{1,\ldots,N\right\}\backslash\left\{i,j\right\}$ such that $g_{k^*}=g_i$ and $g_{l^*}=g_j$. Equation (ref) implies

Differencing (ref)--(ref) and using that $\alpha_{g_{k^*t}}=\alpha_{g_it}$ and $\alpha_{g_{l^*t}}=\alpha_{g_jt}$ yields

a contradiction if groups are “well separated”, i.e., if for all $(g,\widetilde g)\in\left\{1,\ldots,G\right\}^2$ such that $g\neq \widetilde g$, there exists a constant $c_{g,\widetilde g}>0$ such that $$\frac{1}{T}\sum_{t=1}^T(\alpha_{gt}-\alpha_{\widetilde gt})^2 \geq c_{g,\widetilde g}>0.$$ Section (ref) formalizes this equivalence result by establishing uniform asymptotic control over the remainders in the stochastic approximations.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Agglomerative clustering based on thresholding distances}

Given a preliminary consistent estimator $\widehat\beta^1$ and a dissimilarity matrix $\widehat D$ based on $\widehat\beta^1$, the last component of the TPWD estimator is an agglomerative clustering algorithm that builds clusters from the dissimilarity matrix.

The main methodological contribution of this paper is to frame the clustering task as the equivalent high-dimensional multiple-testing problem associated with the null hypothesis $H_{0,i,j}:g_i=g_j$ for all $N(N-1)/2$ distinct ordered pairs of units.\footnote{Similar ideas have recently been developed by Dzemski2024 for different purposes.} If the researcher knew the adjacency matrix of the graph spanned by group memberships, $W\equiv (\mathbf{1}\{ g_i=g_j \})_{(i,j)\in\{1,\ldots,N\}^2}$, which is block diagonal up to some permutation of its rows and columns, then they would know the groups and vice-versa. Although \( W \) is not directly observable in practice, it can be consistently estimated based on the results of pairwise tests.

Specifically, Section (ref) and Lemma (ref) provide sufficient conditions under which any sequence of matrices $\widehat W(c_{NT})\equiv \mathbf{1}\{ \widehat D\leq c_{NT} \}$ with $c_{NT}\to0^+$ converges in max norm to $W$: as $N$ and $T$ tend jointly to infinity, $$\big\Vert\widehat W(c_{NT})-W\big\Vert_{\max}=\max_{(i,j)\in\{1,\ldots,N\}^2}\left\vert\widehat W_{ij}(c_{NT})-W_{ij}\right\vert=o_p(1).$$ The result follows from establishing that both the Type I and Type II errors of the test $\phi_{NT}(i,j)=1-\widehat W_{ij}(c_{NT})$ of $H_{0,i,j}$ tend to zero uniformly across $(i,j)$ (or, equivalently, that both its power and confidence level tend to 1). It motivates using thresholding tests based on the entries of $\widehat D$ and basic agglomerative merging rules that involve a number of operations independent of $G$ and bounded by a polynomial in $N$. In fact, it is even more natural to apply agglomerative clustering algorithms based on merging units whose weigthed pairwise distances fall below some threshold. Such an agglomerative step addresses the practical problem of aggregation arising from the fact that, in finite samples, $\widehat W(c_{NT})$ is generally not block diagonal even after an arbitrary number of permutations of its rows and columns.

This idea lies at the heart of Hierarchical Agglomerative Clustering (HAC) algorithms hastie2009elements, which rely on various “linkage” functions to measure distances between disjoint clusters $A$ and $B$, where $A,B\subset\{1,\ldots,N\}$, based on the dissimilarity matrix $\widehat D$, and various rules to stop the sequential agglomerative merging process, i.e., to cut the induced dendrogram from $N$ singleton clusters, each containing one unit, to a unique cluster of $N$ units, or vice versa for divisive algorithms. Table (ref) summarizes popular linkage functions.

The main theoretical contribution of this paper is to be the first to highlight that such HAC-type algorithms can be successfully applied to GFE models and to provide a distance matrix and cut-off rule which, together with sufficient regularity conditions, yield clustering consistency as $\min(N,T)\to\infty$.

Algorithm 2 below starts with $N$ singleton groups and iteratively merges the closest groups based on their linkage value, continuing until only one group remains or the smallest linkage value exceeds a thresholding parameter $c_{NT}\in[0,\infty)$. This process guarantees a final partition into $\widehat G\in\{1,\ldots, N\}$ non-empty groups. In case of ties, the pair $(j^\star,l^\star)$ with the smallest $j^\star$ and smallest $l^\star$ is selected for merging.\footnote{Asymptotically, or if observed variables are continuously distributed, ties do not occur.}

Algorithm 2 can be efficiently implemented using built-in functions available in popular software (e.g., hclust and cutree in R, or cluster and linkage in MATLAB).

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{A three-step triad pairwise-differencing estimator} The TPWD estimator $\widehat\theta\in\Theta$ of $\theta\in\Theta$ is obtained as follows. Fix $(\psi_{NT},c_{NT})\in[0,\infty)^2$, choose a linkage function (e.g, from Table (ref)), and perform the following steps.

If $\mathcal B=\mathbb R^K$ and $\mathcal A=\mathbb R$, the projection step is a pooled OLS regression of $y_{it}$ on $x_{it}$ and the interactions of estimated groups and time dummies.

{ Remark 1 (Regularization Path):} Consider a finite sample of fixed dimensions $N$ and $T$. If all random variables except group memberships are continuous, then as $c_{NT}\to0$, $\widehat G\to N$ and each group contains a single unit in the limit. Conversely, as $c_{NT}\to +\infty$, $\widehat G\to1$ and a single group contains all units in the limit. Given the low computational cost of the method, a regularization path can be obtained by varying $c_{NT}$ between these extremes. Notably, the clustering step can be vectorized, significantly reducing the computational burden compared to iterative loop-based approaches.

{ Remark 2 (Choice of Tuning Parameters):} Section (ref) offers theoretical guidance on selecting $\psi_{NT}$ and $c_{NT}$. Section (ref) proposes a simple data-driven selection rule that performs well across various Monte Carlo simulations. To further enhance the finite-sample performance, the first step can be re-run with $\widehat v_{it}=y_{it}-x_{it}'\widehat\beta$ (instead of using $\widehat v_{it}= y_{it}-x_{it}'\widehat\beta^1$) to obtain new group assignments $\widehat g_1,\ldots,\widehat g_N$. The second and third steps are then repeated, and the process iterated until convergence. The asymptotic results presented in the following section hold for all subsequent iterations, and Monte Carlo simulations further suggest that the iterative procedure achieves notable improvements in precision. Simulations also suggest that removing the tuning parameter $\psi_{NT}$ by using the NN estimator as the preliminary estimator, instead of the NNR estimator, \[ \widehat\beta^1 \in\operatorname*{arg\,min}_{\beta\in\mathbb R^K}\Vert Y-\beta\cdot X\Vert_1, \] leads to nearly identical results.

{ Remark 3 (Computation):} The full estimation procedure requires $O(N^3T)$ operations, which is a substantial improvement relative to the NP-hard $k$-means problem underlying BM2015's GFE estimator. Whether the computational cost of a consistent clustering algorithm could be further reduced, for example to $O(N^2T)$, remains to the best of my knowledge an open question. While a current limitation of the method is rather the memory required to store triad differences of residuals (looping over $N(N-1)/2$ rather than considering a full vectorization when $N$ is large), it seems hardly possible to cluster units without some measure of distance between them. When unobserved heterogeneity is assumed to be time-constant, the computation cost can be reduced from $O(N^3T)$ to $O(N^2T)$, and the preliminary estimator can be replaced with any standard fixed effects differencing estimator arellano_bond_1991,wooldridge2010econometric; see Section S.2 of the Supplementary Material for further discussion.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Large-sample properties}

In this section, I establish the consistency and asymptotic normality at the parametric rate of the TPWD estimator introduced in Section (ref), under conditions similar to those in the existing literature. Consider the data generating process:

where $g_i^0\in\left\{1,\ldots,G^0\right\}$ denotes group membership, and the $0$ superscripts refer to true parameter values. The asymptotic framework is such that $N$ and $T$ diverge jointly to infinity, denoted by $\min(N,T)\to\infty$. The number of groups $G^0$ is fixed relative to $(N,T)$ but unknown. The discussion on the case of an increasing sequence $G^0=G^0_{NT}$ is relegated to Section S.1 of the Supplementary Material.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Consistency of estimated group memberships} Consider the following assumptions.

Assumption (ref) requires $\widehat \beta^1$ to be consistent for $\beta^0$ at a rate bounded below by $r_{NT}^{-1}$. This rate can be slow. Examples of computationally simple estimators satisfying this condition under low-level conditions are provided in Sections (ref) and Appendix Section (ref). Assumption (ref) ensures that both Type I and Type II classification errors are minimized during the clustering step. It is satisfied for the average, complete, and single linkage functions displayed in Table (ref). While the specific choice of a linkage function is asymptotically innocuous, it may matter in finite samples (see, e.g., Table (ref)). Assumption (ref) allows $T$ to grow considerably more slowly than $N$ (if $\nu\gg1$). It requires the sequence of clustering tuning parameters to vanish, but not too quickly, at a rate of convergence bounded below by $r_{NT}^{-1}$ and strictly slower than $T^{1/2}$. Probability limits are used because this tuning parameter is allowed to be data-driven and, as such, can be random. Assumptions (ref)(ref)--(ref) and (ref)(ref) collect standard moment, tail, and decaying temporal dependence conditions. These assumptions do not impose homoskedasticity but only require uniform bounds on the unconditional variances. Assumption (ref)(ref) requires groups to be well-separated. Although it is strictly weaker than the “strong factors” assumption commonly found in the interactive fixed effects literature, it may break down in very large panels, e.g., in the presence of converging trends.\footnote{Consider the interactive fixed effect version of the model described in Footnote (ref), where $\lambda_i$ belongs to the set of vertices of $\mathbb R^{G^0}$ and $f_t$ collects the group-specific effects at time $t$. Common strong factor assumptions impose $\frac1N\sum_{i=1}^N\lambda_i\lambda_i'\overset{p}{\to}\Sigma_\Lambda>0$ and $\frac1T\sum_{t=1}^Tf_tf_t'\overset{p}{\to}\Sigma_F>0$. Letting $e_g$ denote the $G^0$-vector with 1 in the $g$th coordinate and 0 everywhere else, and $e_{g,\tilde g}:=e_g-e_{\tilde g}$, the continuous mapping theorem implies

Hence, groups are well separated under strong factor assumptions. However, if the time trend of a given group is a homothety of another, the strong factor assumption fails but the two groups will be well-separated as long as the homothety is not the identity. Similarly, \[ \frac{1}{N}\sum_{i=1}^N\mathbf 1\{g_i^0=g\}=e_g'\left(\frac1N\sum_{i=1}^N\lambda_i\lambda_i'\right)e_g\overset{p}{\to}e_g'\Sigma_\Lambda e_g>0. \] Hence, Assumption (ref)(ref) is weaker than the strong factor assumption.} Assumption (ref)(ref) rules out cross-sectional correlation in the error term. Assumption (ref)(ref) allows for asymptotically negligible groups, but requires that each group contains at least two members with probability approaching one. Assumption (ref)(ref) is BM2015's Assumption 2(e). It holds if covariates have bounded support or satisfy dependence and tail conditions similar to those of $v_{it}$. All results below are understood up to group relabelling.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Asymptotic distribution} The next assumption is useful to establish the asymptotic distribution of $\widehat\beta$ and $\widehat{\alpha}_{gt}$.

Assumption (ref) ensures that the infeasible least squares estimator has a standard asymptotic distribution. Assumption (ref)(b) is satisfied if the $x_{it}$ are strictly exogenous or predetermined and observations are independent across units. As a special case, lagged outcomes may thus be included in $x_{it}$. The assumption does not allow for spatial lags such as $y_{i-1t}$.

Consistent plug-in estimates of the asymptotic variances can easily be constructed BM2015.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Choice of tuning parameters} Assumptions (ref) and (ref) provide theoretical guidance for selecting the tuning parameters $\psi_{NT}$ and $c_{NT}$. This guidance, however, is based on asymptotic approximations. In applications with finite sample sizes, an infinite number of choices consistent with the theory will lead to different results. This issue is analogous to bandwidth selection in nonparametric density estimation, or to selecting the number of factors or groups in interactive fixed effects panel data models using AIC and BIC criteria as done in NgBai_2002 and BM2015. The latter involves two tuning parameters: an upper bound for the number of factors or groups and a penalization function characterized by its asymptotic behaviour.\footnote{The eigenvalue ratio test of AhnHorenstein2013 requires only a known upper bound, but it validity is proved under rectangular array asymptotics.}

A data-driven rule is proposed below, whose advantage is to eliminate the risk of specifying an incorrect upper bound for the number of groups. This type of misspecification may lead to poor inference, both in finite samples and asymptotically, as demonstrated in Section (ref).

The selection rule for the thresholding parameter $c_{NT}$ assumes the availability of a consistent estimator $\widehat{\beta}^1$ for $\beta^0$, with a rate of convergence of at least $r_{NT}^{-1}\equiv\sqrt{\min(N,T)}/\log(\log(T))$. By Proposition (ref) in Appendix Section (ref), the nuclear-norm regularized estimator $\widehat\beta^1(\psi_{NT})$, with $\psi_{NT}=\log(\log(T))/\sqrt{16\min(N,T)}$, satisfies this requirement under weak conditions.\footnote{The values obtained for $\widehat\beta^1(\psi_{NT})$ across all Monte Carlo simulations are very close to the nuclear norm estimates, whose computation does not require any tuning parameter.}

The selection rule for the thresholding parameter $c_{NT} $ also requires an estimate of the noise dispersion. Under cross-sectional homoskedasticity in the sense that \[ \frac{1}{T}\sum_{t=1}^Tv_{it}^2\overset{p}{\to}\sigma^2_v, \quad i=1,\ldots, N, \] as $\min(N,T)\to\infty$, an estimator of $\sigma^2_v$ is \[ \widehat \sigma_v^2\equiv \min_{(i,j)\in\{1,\ldots,N\}^2, i\neq j}\frac{1}{2T}\sum_{t=1}^T(\widehat v_{it}-\widehat v_{jt})^2. \] It is not difficult to show that $\widehat\sigma_v^2\overset{p}{\to} \sigma_v^2$ under the maintained assumptions Zeleneev2020. It turns out, however, that this estimator exhibits a substantial downward finite-sample bias, which converges to zero very slowly in Monte Carlo simulations: prohibitively large values of $T$ are needed to obtain reasonable performance. An asymptotically valid method to diminish the bias is to consider the max-min estimator \[ \check \sigma_v^2\equiv \max_{i\in\{1,\ldots,N\}}\min_{j \in\{1,\ldots,N\}, j\neq i}\frac{1}{2T}\sum_{t=1}^T(\widehat v_{it}-\widehat v_{jt})^2. \] Again, it is not difficult to show that $\check\sigma_v^2\overset{p}{\to} \sigma_v^2$ under cross-sectional homoskedasticity and the maintained assumptions. This estimator yields important improvements in Monte Carlo simulations with small to moderate sample sizes. The data-driven selection rule for the thresholding parameter $c_{NT} $ is then \[c_{NT}\equiv\frac{1.35\times\check\sigma_v\times\log(T)}{\max(K,1)\sqrt{\min(N,T)}}. \] The theoretical validity of this choice follows from $\check \sigma_v=O_p(1)$ under the maintained assumptions, so that $r_{NT}$ and $c_{NT}$ verify Assumptions (ref) and (ref).\footnote{An alternative approach could be based on limiting overfitting for different value of $c_{NT}$ as is done in BonhommeLamadonManresa2022. This requires computing a regularization path and, therefore, is computationally more costly.}

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Monte Carlo simulations} This section reports the results of some Monte Carlo simulations calibrated to the empirical application. The finite sample properties of the TPWD estimator are assessed across two distinct data generating processes (DGP): a pure GFE model without covariates and a full GFE model with a scalar covariate. In each DGP, the variance of the error is set to the one from BM2015's calibrated simulations. Grouped fixed-effects and group memberships are set so that the minimum distance between groups matches those of BM2015's calibrated simulations for $G=4$. The signal-to-noise ratios of the group membership variables match approximately across all DGPs. The simulation design is not identical in order to let the sample size grow.\footnote{BM2015's simulations are based on the dataset from the empirical application, where $N=90$ and $T=7$.} The sample sizes and number of groups vary across $(N,T)\in\{90, 180\} \times \{7, 10, 20, 40\}$ and $G\in\{3, 4\}$, respectively.

To mirror the empirical findings of BM2015, where the observed outcome is a scalar, continuous measure of a country's democratic regime status at time $t$, the grouped fixed effects are set as follows: for all $t\in\{1,\ldots,T\}$, \[ \alpha_{1t}\equiv 1, \quad \alpha_{2t}\equiv \frac{t-1}{T-1}, \quad \alpha_{3t}\equiv0, \quad \alpha_{4t}\equiv \frac{\mathbf{1}\{ t\geq\lfloor T/2\rfloor \}(t-\lfloor T/2 \rfloor)}{T-\lfloor T/2\rfloor}. \] The first group could be labeled as “high democracy,” the second as “early transition,” the third as “low democracy,” and the fourth as “late transition.” In the baseline setting, groups are balanced: \[ g_i = 1 + \sum_{g=1}^{G-1}\mathbf{1}\{ i > g\lfloor N/G\rfloor \}, \quad i=1,\ldots,N. \]

The performance of the TPWD estimator is compared with benchmark alternatives: the grouped fixed-effects (GFE) estimator from BM2015, the spectral and post-spectral estimators from ChetverikovManresa2021, the nuclear-norm (NN) and nuclear-norm regularized (NNR) estimators from moon2019nuclear, the least squares interactive fixed effects (IFE) estimator from Bai2009, and the infeasible pooled OLS regression that uses the true group memberships (Oracle).

Sections S.3.1--3.5 contain additional results for models with unbalanced groups, more groups, unit-specific effects, higher signal-to-noise ratio, or time-invariant unobserved heterogeneity. In all experiments, results are averaged across $500$ Monte Carlo samples.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Pure grouped fixed effects model} Consider a pure GFE model without covariates: $$y_{it}=\alpha_{g_it}+v_{it}, \quad i=1,\ldots,N,\; t=1,\ldots,T,$$ where $v_{it}\overset{iid}{\sim}\mathcal N(0,(1/3)^2)$ across units and time periods. Let $\widehat\pi_g\equiv N^{-1}\sum_{i=1}^N\mathbf{1}\{ g_i=g \}$. Since groups are balanced, $\widehat\pi_g=1/G$ and the average signal-to-noise ratio of group membership is \[ \frac{\sum_{g=1}^G{\rm Var}(\mathbf{1}\{ g_i=g \})}{\sigma_v^2}=\frac{\frac{1}{G}\sum_{g=1}^G\widehat\pi_g(1-\widehat\pi_g)}{(1/3)^2} =2\mathbf{1}\{ G=3 \}+\frac{27}{16}\mathbf{1}\{ G=4 \}. \] This is to compare with $1.89\mathbf{1}\{ G=3 \}+1.69\mathbf{1}\{ G=4 \}$, the average signal-to-noise ratio in BM2015's calibrated simulations.

The TPWD estimator with $\widehat\beta^1=0$ and $c_{NT}=1.35\check\sigma_v\log(T)/\sqrt{\min(N,T)}$ is compared to versions without covariates of ChetverikovManresa2021's post-spectral estimator with $g\in\{2,3,4,10\}$ user-specified number of groups (Post-Spectral$^{\bar G=g}$) and BM2015's grouped fixed-effects estimator with $g\in\{2,3,4,10\}$ user-specified number of groups, i.e., their “Algorithm 1” with 1,000 randomly generated initialization points (GFE$^{\bar G=g}$).

The performance of all estimators is assessed in terms of the average estimated number of groups (if applicable), $500^{-1}\sum_{b=1}^{500}\widehat G^{(b)}$, the average root mean square error (RMSE) of the estimated grouped fixed effects,

and the average clustering accuracy, measured using the average Precision (P) rate, Recall (R) rate, and Rand Index (RI): \[ P(\widehat g)\equiv\frac{1}{500}\sum_{b=1}^{500}p(\widehat g^{(b)}), \quad R(\widehat g)\equiv\frac{1}{500}\sum_{b=1}^{500}r(\widehat g^{(b)}), \quad RI(\widehat g)\equiv\frac{1}{500}\sum_{b=1}^{500}ri(\widehat g^{(b)}), \] where, by denoting the number of false positives (FP), true positives (TP), false negatives (FN), and true negatives (TN) as

the following three measures of clustering accuracy are invariant to cluster relabelling:

The precision rate measures the proportion of correctly matched unit pairs among all matched pairs. The recall rate measures the proportion of correctly matched unit pairs among all pairs of units belonging to the same population cluster. The Rand index measures the proportion of correct decisions (to match or not to match) among all matching decisions taken by the clustering algorithm (i.e., among all possible pairs of units). It summarizes both Type I and Type II errors in the prediction of group memberships.\footnote{Given the null hypothesis $H_{0,i,j}:g_i=g_j$, $\forall (i,j)$, introduced earlier, the definitions of $FP,TP,FN$, and $TN$ should be flipped, e.g., $FP(\widehat g)\equiv\sum_{i<j}\mathbf{1}\{ \widehat g_i\neq \widehat g_j \} \mathbf{1}\{ g_i= g_j \}$. Conceptually, however, a merging HAC algorithm initialized at $N$ singleton groups “starts” with the null $\widetilde H_{0,i,j}:g_i\neq g_j$, $\forall (i,j)$, while a divisive algorithm initialized at a unique group of $N$ units “starts” with $H_{0,i,j}$, $\forall (i,j)$. The multiple-testing problem was stated that way following the standard practice of stating the null hypothesis as an equality. }

Table (ref) reports the average estimated number of groups (if applicable), the RMSE of the estimated grouped fixed effects, and the execution time in seconds for each estimator. The performance of TPWD in terms of RMSE uniformly dominates that of Post-Spectral$^{\bar G=g}$ across all values of $(G,N,T)$ and user-specified $g\in\{2,3,4,10\}$.\footnote{Note that Post-Spectral$^{\bar G=10}$ cannot be computed for $T=7$ because it requires computing the first ten largest eigenvalues of a $7\times7$ matrix.} For small $T\in\{7,10\}$, the RMSE of the well-specified (green-shaded) Post-Spectral$^{\bar G=G}$ is about twice as large as that of TPWD. Except for $T=7$ or $(G,N,T)\in\{$(3,180,10), (4,90,10),(4,180,10)$\}$, the performance of TPWD in terms of RMSE dominates that of each misspecified GFE estimator. The performance of TPWD in terms of RMSE is close but dominated by that of the well-specified (green-shaded) GFE$^{\bar G=G}$ across almost all values of $(G,N,T)$. For $G=3$ (resp. $G=4$), it is within a $0.071$ (resp. $0.020$) distance.

For $(G,N,T)=(3,90,7)$, the TPWD estimator is computed in $1/20$ seconds, has an RMSE of $0.150$, and estimates $6.654$ groups on average. The well-specified post-spectral estimator, Post-Spectral$^{\bar G=3}$, has an RMSE of $0.305$ and estimates $2.838$ groups on average, the well-specified GFE estimator, GFE$^{\bar G=3}$, has an RMSE of $0.088$, the misspecified GFE estimators, GFE$^{\bar G=2}$, GFE$^{\bar G=4}$, and GFE$^{\bar G=10}$, have RMSEs of $0.270$, $0.129$, and $0.216$, respectively, and the Oracle estimator has an RMSE of $0.060$. As $T$ grows, the number of groups estimated by TPWD converges to the ground truth, and its RMSE steadily decreases to reach the same level of performance as the Oracle estimator, with precision 10$^{-2}$ provided $T\geq 20$ (for $G=3$) or $T\geq 40$ (for $G=4$). While GFE estimators exhibit a similar pattern, post-spectral estimators reach this level of precision only for $G=3$ and $T\geq 40$. Unreported simulations suggest that GFE estimators are quite sensitive to the choice of initializers: if those are randomly drawn around zero (e.g., from a normal distribution), convergence breaks down in the simulations. This is likely because the clustering problem solved by GFE estimators becomes more difficult due to its NP-hardness as $G$ and $T$ grow. In particular, when $G=3$, the well-specified GFE estimator never attains the Oracle performance and does not seem to converge with $T$ (numerical errors dominate statistical errors). When $G=4$, its convergence is slow so that it does not attain the Oracle performance and similar numerical challenges can be expected for larger values of $T$.

For $(G,N,T)=(3,90,40)$, the TPWD estimator is computed in $1/3$ seconds, has an RMSE of $0.061$, and estimates $3.012$ groups on average. The well-specified post-spectral estimator, Post-Spectral$^{\bar G=3}$, has an RMSE of $0.062$ and estimates $3.000$ groups on average, the well-specified GFE estimator, GFE$^{\bar G=3}$, has an RMSE of $0.061$, and the misspecified GFE estimators, GFE$^{\bar G=2}$, GFE$^{\bar G=4}$, and GFE$^{\bar G=10}$, have RMSEs of $0.243$, $0.084$, and $0.143$, respectively. For $(G,N,T)=(4,90,40)$, the TPWD estimator is computed in $1/3$ seconds, has an RMSE of $0.077$, and estimates $3.986$ groups on average. In comparison, Post-Spectral$^{\bar G=4}$ has an RMSE of $0.119$ and estimates $3.472$ groups on average, and GFE$^{\bar G=4}$, GFE$^{\bar G=2}$, GFE$^{\bar G=3}$, GFE$^{\bar G=10}$, have RMSEs of $0.072$, $0.227$, $0.117$, and $0.140$, respectively. Similar patterns are observed for $N=180$, where TPWD requires at most 4 seconds to compute.\footnote{Experiments were run in parallel on the computing cluster of ENSAE Paris, but the estimator itself was not parallelized and takes a similar amount of time to compute on any professional laptop.}

In some settings, the well-specified post-spectral estimator gets the number of groups closer (e.g., if $G=3$), but this is by construction since the number of groups is an input of the algorithm. The TPWD estimator may slightly overestimate or underestimate the number of groups, depending on the choice of the thresholding parameter in finite samples, but this has little impact on the RMSE. Intuitively, similar units are grouped together based on a finer grid, rather than by imposing a fixed value of \( G \).

The results of Table (ref) are best understood by examining the performance of group membership estimates. Table (ref) displays the average Precision rate, Recall rate, and Rand index of all reported estimators of group memberships. For all user-specified $g\in\{2,3,4,10\}$, the performance of Post-Spectral$^{\bar G=g}$ in terms of Precision rate and Rand index is dominated by that of TPWD across all values of $(G,N,T)$. The Rand index of TPWD is most often slightly dominated by but within a $0.141$ distance of that of the well-specified GFE estimator across all values of $(G,N,T)$. This is not the case for post-spectral estimators. For $G=4$, post-spectral estimators have a better Recall rate than TPWD, which is partly explained by the fact that they take the number of groups as input and contrain the estimated number of groups to be less or equal to it. This is also the case for GFE$^{\bar G=g}$ with $g\leq 4$. Similarly, misspecified GFE estimators with $g> G$ almost always achieve the best Precision rate, as they allow for finer clustering at the expense of lower recall. TPWD provides more accurate clustering in terms of the Rand Index by achieving a better balance between precision and recall, while simultaneously determining the number of groups endogenously as a by-product. Although it may occasionally identify a slightly larger number of groups (leading to lower recall), the units within these groups tend to share the same true group membership (resulting in higher precision). As shown in Table (ref), this trade-off is sufficient to yield a lower RMSE.

\@startsection{subsection}{2}{0mm}{-1.2\baselineskip}{1\baselineskip}{\normalfont}{Full grouped fixed effects model} Consider now a full GFE model with a scalar covariate: $$y_{it}=x_{it}\beta+\alpha_{g_it}+v_{it}, \quad i=1,\ldots,N, \; t=1,\ldots,T,$$ where $\beta=1$, $x_{it}=0.5\alpha_{g_it}+u_{it}$, $u_{it}\overset{iid}{\sim}\mathcal N(0,1/(2\sqrt{3})^2)$, $v_{it}\overset{iid}{\sim}\mathcal N(0,(1/3)^2)$, and $u_{it}$ and $v_{it}$ are mutually independent across units and time periods. In the dataset from the empirical application and BM2015's calibrated simulations, the lagged income covariate has variance 1.152, resulting in a signal-to-noise ratio of 10.368. Here, the variance of the covariate is less than 0.1, which for $(N,T)=(90,7)$ leads to an average signal-to-noise ratio for the slope coefficient of $$9\left(\frac{0.25}{T}\sum_{t=1}^T\frac1G\sum_{g=1}^G\widehat\pi_g(1-\widehat\pi_g)\alpha_{gt}^2+\frac{1}{12}\right)=0.9769\mathds1\{G=3\}+0.9218\mathds1\{G=4\},$$ i.e., more than ten times smaller. The level of correlation between the covariate and the fixed effect is similar in both papers.

The performance obtained after one (TPWD) or four (Iterated TPWD) iterations of the TPWD estimator which uses $\widehat\beta^1(\psi_{NT})$ -- the NNR estimator described in Section (ref) with $\psi_{NT}=\log(\log(T))/\sqrt{16\min(N,T)}$ -- and the data-driven thresholding rule $c_{NT}=1.35\check\sigma_v\log(T)/\sqrt{\min(N,T)}$ is compared with that of the NNR and NN estimators proposed in moon2019nuclear, the well-specified spectral and post-spectral estimators proposed in ChetverikovManresa2021, the GFE estimator (with or without BIC selection of $G$) proposed in BM2015, the well-specified IFE estimator proposed in Bai2009, initialized at NNR or random draws, and the infeasible pooled OLS regression that uses the true group memberships (Oracle). I refer to the original papers, table notes, and the code made available online for details regarding the implementation of each estimator. For each estimator, the bias and root mean square error (RMSE) of the estimated slope coefficient are reported:

When applicable, the coverage rate of a $95\%$-level confidence interval for $\beta$ based on a consistent estimator of the large-$N$, large-$T$ asymptotic variance clustered at the unit level, the RMSE of grouped fixed effect estimates, the estimated number of groups, the Precision and Recall rates, and the Rand index of the estimated clustering are reported.

Table (ref) shows the results in terms of bias, RMSE, estimated number of groups, and coverage rates. These findings confirm the theoretical guarantees stated in Proposition (ref) and Corollary (ref). As the sample size increases, both the TPWD and iterated TPWD estimators exhibit consistency, and the coverage rates of the confidence intervals -- constructed using a consistent estimator of the variance of their asymptotic normal distribution -- converge to their nominal levels. Remarkably, for small values of $T$, four iterations significantly reduce bias and improve coverage. As noted in moon2019nuclear and ChetverikovManresa2021, the convergence of the NNR and NN estimators is indeed slow. Yet, contrary to what was conjectured in ChetverikovManresa2021, this does not prevent the iterated TPWD estimator to reach near-Oracle performance even for small values of $T$ for which post-spectral estimators perform quite poorly. Remarkably, the iterated TPWD estimator uniformly dominates the NNR, NN, and well-specified post-spectral estimators across all values of $(G,N,T)$, in terms of all metrics. While its bias is sometimes slightly greater -- and often slightly smaller -- than that of the spectral estimator (which relies on the factor structure assumption for the covariates), the difference remains within a margin of 0.03. However, its RMSE is uniformly smaller across all values of \( (G, N, T) \) except for $G=4$ and $T\in\{20,40\}$. In the following, the focus is exclusively on the iterated TPWD, which will be referred to simply as TPWD.

For $(G,N,T)=(3,90,7)$, TPWD has a bias of $0.028$ and its coverage rate already peaks at $80.8\%$. The bias of the well-specified post-spectral estimator is $0.447$ and its coverage rate is $3.4\%$. The well-specified GFE estimator has a bias of $0.013$ and coverage rate of $90.8\%$. The infeasible oracle estimator has a bias less than $10^{-3}$ and displays a coverage rate of $93.4\%$. The NNR and NN estimator have biases of $0.732$ and $0.736$ respectively. As $T$ increases, the performance of each estimator steadily improves, and estimates of the number of groups converge to the true number of groups. For $(G,N,T)=(3,90,20)$, TPWD has a bias of $0.001$ and coverage rate of $93.2\%$. The bias of the well-specified post-spectral estimator is $0.041$ and its coverage rate is $77.8\%$. Both the well-specified GFE and Oracle estimator have a bias less than $10^{-3}$ and coverage rate of $93\%$. The NNR and NN estimator have biases of $0.417$ and $0.418$ respectively. As $N$ increases, the discrepancy in the performance of the TPWD and post-spectral estimators remains stable. For $(G,N,T)=(3,180,10)$, TPWD has $88.8\%$ coverage to contrast with $13\%$ and $92.8\%$ for the post-spectral and GFE estimators respectively. For $(G,N,T)=(3,180,20)$, TPWD has $94.6\%$ coverage to contrast with $86\%$ and $94.8\%$, respectively.

While the coverage rate converges to the prescribed level for all estimators as $N$ and $T$ reach their maximum value if $G=3$, this is not the case for the post-spectral estimator if $G=4$. For $(G,N,T)=(4,90,7)$, TPWD has $66.2\%$ coverage to contrast with $18\%$ and $80.8\%$ for the post-spectral and GFE estimators respectively. For $(G,N,T)=(4,90,40)$, TPWD has $91.4\%$ coverage to contrast with $27\%$ and $93.4\%$, respectively. For $(G,N,T)=(4,180,7)$, TPWD has $75.8\%$ coverage to contrast with $11.8\%$ and $74.8\%$, respectively. For $(G,N,T)=(4,180,40)$, TPWD has $86.4\%$ coverage to contrast with $10\%$ and $94.6\%$, respectively.

In line with Corollary 1 in Bai2009, in this setting with i.i.d. errors and homoskedasticity, the performance of the well-specified IFE estimators is comparable to that of the GFE estimator in terms of bias, even for small values of $T$, but the former displays systematic under-coverage in comparison with the latter. Similarly, the BIC selection approach of BM2015 performs quite well, as conjectured by the authors. These three approaches, however, are not immune to misspecifying $G$ or the upper bound $G_{\max}$.

Table (ref) reports clustering accuracy metrics. The three measures are little affected by the introduction of a single covariate and the slow rate of convergence of the preliminary consistent estimator, even for small values of $T$ (compare with Table (ref)).

In summary, given that I am not aware of any theoretical result that would guarantee consistency of the BIC selection approach for some of the asymptotic regimes considered in this paper (with possibly $T<<N$), that naive cross-validation for $G$ is not recommended for unsupervised clustering algorithms tasks with sample dependent parameters, and that heuristics such as the “elbow method” or “gap curve” are generally not theoretically justified hastie2009elements, I would recommend to use the TPWD estimator for inference, especially when $T$ is small, $N$ is large, or $G$ is expected to be large. In particular, TPWD scales remarkably well. In an unreported Monte Carlo simulation with $N=2000$, $T=7$, and $G=4$, which may be typical in microeconometric applications, the TPWD point estimate of $\beta=1$ is $1.015$ and it takes less than 2 minutes to compute.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Empirical application: income and (waves of) democracy} Understanding the statistical relationship between income and democracy has been a longstanding issue in political science and economics Lipset1959,Barro1999. Using panel data for $N=90$ countries observed at $T=7$ points in time over the period 1970--2000, Acemoglu2008 found that the statistically significant positive association between income on democracy vanishes when country fixed effects are included in the regression. They argued that these results are consistent with countries having embarked on divergent paths of economic and political development at certain points in history, or critical junctures. Some of the examples they mention are the end of feudalism, the industrialization age, or the process of colonization. In this perspective, the fixed effects are meant to capture these highly persistent historical events. BM2015 proposed to test this assumption by computing an approximation of the GFE estimator using alternative minimization and reporting the results for several numbers of groups. They argued that the true number of groups would be less than $10$, reporting statistically significant income effects -- elasticities of a measure of democracy to lagged income per capita -- between $0.061$ and $0.089$ .

This section provides a reassessment of their results, consistently estimating the number of groups by applying the TPWD estimator to their preferred specification: a regression model of democracy (measured by the Freedom House indicator) on lagged democracy and lagged log-GDP per capita with unrestricted group-specific time patterns of heterogeneity $\alpha_{g_it}$:

The data is obtained from the balanced subsample of Acemoglu2008.\footnote{Available at: \href{https://www.aeaweb.org/articles?id=10.1257/aer.98.3.808}{https://www.aeaweb.org/articles?id=10.1257/aer.98.3.808}.} The preliminary estimator is a nuclear-norm regularized (NNR) estimator with tuning parameter set to the theoretically valid rule $\psi_{NT}=\log(\log(T))/\sqrt{16\min(N,T)}$. The data-driven thresholding rule of the TPWD estimator is as described in Section (ref).

Table (ref) displays the NNR estimates, the TPWD$^{k{\rm it}}$ estimates at iteration $k\in\left\{1,2,3,4\right\}$ (until convergence), and the GFE$^{\bar G=g}$ estimates with user-specified number of groups $g\in\left\{2,3,10\right\}$. After three iterations, the TPWD estimator converges to four estimated groups and delivers a significant income effect of $0.070$. This point estimate is relatively close to the GFE$^{\bar G=2}$ and GFE$^{\bar G=10}$ estimates of $0.061$ and $0.075$, respectively. The estimated cumulative income effect ($\beta_2/(1- \beta_1)$) is 0.258, also significant, and more than twice the GFE$^{\bar G=10}$ estimate of $0.104$. The preliminary NNR estimator delivers point estimates of $0.016$ and $0.078$, respectively.

\@startsection{section}{2}{0mm}{-1.5\baselineskip}{1\baselineskip}{\normalfont}{Conclusion} Grouped fixed effects models are plagued with an underlying combinatorial classification problem, rendering estimation and inference difficult. This paper proposes a novel strategy for the constructive identification of all model parameters, including the number of groups. The method simultaneously solves the model selection, classification, and estimation problems. The corresponding three-step estimator has polynomial computational cost and is straightforward to implement, requiring only smooth convex optimization and elementary arithmetic operations. It builds on an initial off-the-shelf consistent estimator of the slope coefficients and applies thresholding to suitable pairwise-differencing transformations of the residualized regression equations. Under mild conditions, the proposed estimator is shown to be uniformly consistent for the latent grouping structure and asymptotically normal as both dimensions of the panel grow jointly. Importantly, the number of groups is consistently estimated without prior knowledge of its support, and the time dimension may grow at a much slower rate than the cross-sectional dimension.

Beyond its stronger large-sample properties -- established under relatively weaker assumptions than in the existing literature -- Monte Carlo simulations demonstrate its finite-sample competitiveness, if not superiority, compared to spectral clustering and grouped fixed effects estimators.

Several open questions remain for future research: Could this approach be used to construct a formal test of the grouping assumption? Can similar differencing strategies be extended to nonlinear structural models or potential outcome frameworks? Might agglomerative methods help address the problem of weak factors in latent group structures? And is it possible to develop finite-sample or uniform inference procedures?