# A Simple and Computationally Trivial Estimator for Grouped Fixed Effects
Models111An earlier version of this paper was circulated under the title “Make
the Difference! Computationally Trivial Estimators for Grouped Fixed Effects
Models”. I thank Stéphane Bonhomme, Xavier D’Haultfœuille, Elena Manresa,
Pauline Rossi, Ao Wang, Martin Weidner, Andrei Zeleneev, and seminar
participants at CREST, LMU Munich, UChicago, Oxford EET 2023, and Bristol ESG
2022 for helpful comments. This research is supported by the French National
Research Agency grants ANR-17-CE26-0015-041, ANR-18-EURE-0005,
ANR-11-LABX-0047, and the European Research Council grant
ERC-2018-CoG-819086-PANEDA.

Martin Mugnier Department of Economics, University of Oxford,
<EMAIL_ADDRESS>
###### Abstract

This paper introduces a new fixed effects estimator for linear panel data
models with clustered time patterns of unobserved heterogeneity. The method
avoids non-convex and combinatorial optimization by combining a preliminary
consistent estimator of the slope coefficient, an agglomerative pairwise-
differencing clustering of cross-sectional units, and a pooled ordinary least
squares regression. Asymptotic guarantees are established in a framework where
$T$ can grow at any power of $N$, as both $N$ and $T$ approach infinity.
Unlike most existing approaches, the proposed estimator is computationally
straightforward and does not require a known upper bound on the number of
groups. As existing approaches, this method leads to a consistent estimation
of well-separated groups and an estimator of common parameters asymptotically
equivalent to the infeasible regression controlling for the true groups. An
application revisits the statistical association between income and democracy.

Keywords: panel data, time-varying unobserved heterogeneity, grouped fixed
effects, agglomerative clustering
JEL Codes: C14, C23, C38.

##  1 Introduction

Suppose a sample of longitudinal or 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:

$\left\\{\begin{array}[]{ll}y_{it}&=x_{it}^{\prime}\beta+\alpha_{g_{i}t}+v_{it}\\\
\mathbb{E}[v_{it}]&=\mathbb{E}[\alpha_{g_{i}t}v_{it}]=\mathbb{E}[x_{itk}v_{it}]=0\end{array}\right.,\quad
i=1,\ldots,N,\;t=1,\ldots,T,\;k=1,\ldots,K,$ (1.1)

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 group-
time-specific effect $\alpha_{g_{i}t}\in\mathbb{R}$.

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

Special cases of factor models with factor loadings confined to a finite
set,222Note that $\alpha_{g_{i}t}=\lambda_{i}^{\prime}F_{t}$ for any
$\lambda_{i}^{\prime}\equiv\left(c_{1}\mathbf{1}\\{g_{i}=1\\},\ldots,c_{G}\mathbf{1}\\{g_{i}=G\\}\right)$,
$F_{t}^{\prime}\equiv\left(\alpha_{1t}/c_{1},\ldots,\alpha_{Gt}/c_{G}\right)$,
and
$c\equiv(c_{1},\ldots,c_{G})^{\prime}\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},\tilde{\lambda}_{i}\in\mathbb{R}^{r}$ for some
$r\in\mathbb{N}\backslash\\{0\\}$ with
$\left|\left\\{\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{N}\right\\}\right|=\tilde{G}$,
then there exist
$(\tilde{g}_{1},\ldots,\tilde{g}_{N})^{\prime}\in\\{1,\ldots,\tilde{G}\\}^{N}$
and
$(\tilde{\alpha}_{11},\ldots,\tilde{\alpha}_{\tilde{G}T})^{\prime}\in\mathbb{R}^{\tilde{G}T}$
such that
$\tilde{\lambda}_{i}^{\prime}\tilde{F}_{t}=\tilde{\alpha}_{\tilde{g}_{i}t}$.
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 (Hahn and Moon,
2010; Bonhomme and Manresa, 2015), GFE have gained considerable interest in
both methodological and applied work (e.g., Su, Shi, and Phillips, 2016;
Cheng, Schorfheide, and Shao, 2019; Bonhomme, Lamadon, and Manresa, 2019; Gu
and Volgushev, 2019; Chetverikov and Manresa, 2021; Mugnier, 2022b; Janys and
Siflinger, 2023).

Treating GFE as interactive fixed effects, however, leads to two main issues.
First, parametric rate inference for the common parameters is generally not
available when $T$ grows very slowly with $N$ (e.g., Bai, 2009; Moon and
Weidner, 2019; Beyhum and Gautier, 2023). Second, although Higgins (2022)
shows that parametric rate inference remains possible under some
circumstances, his proposed method solves a high-dimensional non-convex least-
squares problem typically subject to local minima. Similarly, addressing the
combinatorial nature of GFE 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. Second, GFE estimators typically require $G$ or
$G_{\max}\geq G$ to be known to the analyst.

This paper provides a three-step estimation procedure free of these
limitations. The first step leverages the low-rank factor structure of GFE to
compute a computationally simple estimator of the slope coefficient. I focus
on smooth and convex regularized nuclear norm estimation, but the interactive
fixed effects literature provides many candidates further discussed. The
second step examines the correlation in pairwise differences of residuals
within triads of units to build distances between units that are shown to be
asymptotically zero if and only if two units belong to the same groups,
whenever group-time effects are well separated. An agglomerative clustering
based on thresholding distances, in turn, is natural. The third step computes
a pooled ordinary least squares regression (OLS) of the outcome on the
covariates and interaction of time and estimated group dummies.

While the estimator combines existing disparate ideas from the panel data
literature, this particular combination and its appealing large sample
statistical properties are truly novel. Specifically, I show that model
selection (i.e., estimation of $G$) can be achieved simultaneously without the
need to choose a pre-specified range of misspecified models (an upper bound
$G_{\max}$) but simply merging close units at a rate governed by time
dependence conditions. Since $G$ is unknown, regularization is needed and
there is no free lunch. Bonhomme and Manresa (2015) discuss the use of AIC/BIC
information criteria applied to a computationally difficult estimator but do
not provide asymptotic guarantees within the asymptotic sequences considered.
In contrast, I propose a computationally trivial data-driven procedure which
is shown to be theoretically valid and performs relatively well in finite
samples. Under standard conditions, I show that the first two steps achieve
consistent estimation of well-separated groups. As similar methods, this leads
to an estimator of common parameters asymptotically equivalent to the
infeasible regression controlling for the true groups. Compared to the
existing literature, the new approach is computationally trivial, generates
the number of groups endogenously with a theoretically valid data-driven
selection rule and is based on minimal assumptions on the covariates
(sufficient conditional variation through a restricted eigenvalue condition).
This generalization of Bonhomme and Manresa (2015) comes at the price of
ruling out cross-sectional correlations in the error terms because the new
estimator relies on triad-specific comparisons, instead of cross-sectional
averages, to estimate the grouping structure.

While I focus on a simple linear model with a homogeneous slope to best convey
the core proof ideas, the versatility of the new approach is illustrated by
several extensions with a growing number of groups or covariates,
heterogeneous slope coefficients (in time or units), or multiplicative network
models (e.g., gravity trade equations) discussed in Appendix B. This paper
does not address the question of conducting inference on groups. Dzemski and
Okui (2018) provide pointwise valid inference methods given a preliminary
estimator $\widehat{\beta}$ is available. To this end, one can use the
estimator proposed in the present paper.

I investigate the finite-sample performance of the new estimator across
different Monte Carlo experiments calibrated to the empirical application. I
compare the new estimator to (possibly misspecified) GFE estimators with
several prespecified numbers of groups and the nuclear norm regularized (NN-
reg) estimator. In a pure setting without covariates, the new method
asymptotically outperforms all GFE estimators and achieves homogeneous
clustering performance in terms of Precision, Recall, and Rand Index defined
later. In a full setting with a scalar covariate, the new estimator
significantly improves upon NN-reg in terms of bias and root mean square
error. For $T\geq 20$, it achieves similar performance as the infeasible
pooled OLS regression, and confidence intervals based on an estimator of the
asymptotic variance attain their pre-specified $95\%$ level.

Finally, I showcase the usefulness of the method by revisiting Bonhomme and
Manresa (2015)’s study of the statistical association between income and
democracy for a large panel of countries between 1970-2000. This association
could be confounded by critical junctures in history that led to similar
unobserved development paths for some countries and different for others.
While the authors find statistically significant (approximated) GFE estimates
of the effect of lagged income per capita on a measure of democracy ranging
between 0.061 and 0.089 and suggest that the number of groups is less than 10,
I estimate $3$ groups, an effect of 0.105 and a much higher cumulative income
effect of 0.243 (v.s. 0.104 to 0.151). The preliminary regularized nuclear
norm estimator delivers point estimates of 0.016 and 0.078 respectively.

#### Related Literature.

This paper contributes to a vast literature on estimating interactive fixed
effects panel models (e.g., Pesaran, 2006; Bai, 2009; Moon and Weidner, 2015,
2017; Bonhomme, Lamadon, and Manresa, 2022; Armstrong, Weidner, and Zeleneev,
2023; Beyhum and Gautier, 2023). Convex nuclear norm regularized estimators
(e.g., Moon and Weidner, 2019; Chernozhukov, Hansen, Liao, and Zhu, 2019),
though computationally simpler than non-convex least squares (e.g. Bai, 2009),
may converge slower than parametric rates. I use bounds on the rate of
convergence of such estimators derived in this literature to show how
computationally simple estimators of the special case of grouped fixed effects
models converge at a parametric rate. Chetverikov and Manresa (2021)
independently used similar ideas, although they impose a factor structure on
the covariates and a known upper bound on the number of groups to apply
spectral clustering techniques. Differently, I only rely on a restricted
eigenvalue condition and achieve model selection simultaneously. The
regularization parameter is a data-driven substitute for a known upper bound
$G_{\max}\geq G$ and an AIC/BIC model selection criterion often used to select
the number of factors in interactive fixed effects models (e.g., Bai, 2003).

This paper also contributes to a rapidly growing literature on estimating
grouped fixed effects models. Bonhomme and Manresa (2015)’s GFE estimator, an
extension of $k$-means clustering to handle covariates, solves an NP-hard
optimization problem. Algorithms that provide fast solutions may not converge
to its true value, and the same holds true for extensions and other non-convex
estimators (e.g., Su, Shi, and Phillips, 2016; Ando and Bai, 2022; Mugnier,
2022b; Lumsdaine, Okui, and Wang, 2023). In contrast, the inferential theory
developed in this paper is valid for a computationally simple estimator which
substitutes a known upper bound on the number of groups with willingness to
merge units into groups based on estimated pairwise distances. Similarly,
because inference is on a true population parameter, my results contrast with
Pollard (1981, 1982) which provides an asymptotic theory for the solution to
the population $k$-means sum of squares problem in the cross-sectional case,
i.e., only for a pseudo-true value. While the closest approach of Chetverikov
and Manresa (2021)’s spectral and post-spectral estimators is computationally
straightforward, its implementation requires prior knowledge about the number
of groups (or a consistent estimator which is 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 Pesaran (2006) and could be restrictive in practice.
Differently, I do not impose a model for the covariates (only sufficient
variation) neither an upper bound for the number of groups. Similarly, Lewis,
Melcangi, Pilossoph, and Toner-Rodgers (2022) recently introduced a fuzzy-
clustering procedure that is shown to work well on simulations but whose large
sample properties are not investigated and which requires the number of groups
to be known to the researcher. Mehrabani (2023) adapts “sum-of-norm” convex
clustering (e.g., Hocking, Joulin, Bach, and Vert, 2011; Tan and Witten, 2015)
to a linear panel data model with latent group structure but time-constant
group effects only. Extending his approach to accommodate time-varying effects
seems difficult. Albeit close in spirit, the proposed procedure is different
from the binary segmentation algorithm developed in Wang and Su (2021), or the
pairwise comparisons method proposed in Krasnokutskaya, Song, and Tang (2022).
A different approach would apply spectral clustering to some dissimilarity
matrix (see, e.g., Ng, Jordan, and Weiss, 2002; von Luxburg, 2007; Chetverikov
and Manresa, 2021; Brownlees, Guðmundsson, and Lugosi, 2022; Yu, Gu, and
Volgushev, 2022). A disadvantage is that it adds one layer of
complexity/tuning parameters: $L$ eigenvectors of the dissimilarity matrix
have to be computed and clustered (and $L$ chosen), which is usually done by
approximating a NP-hard $k$-means solution that this paper aims to avoid 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. Though agglomerative (hierarchical)
clustering methods are well known in the former, adapting them to the latter
is relatively new (a few exceptions in different models being Vogt and Linton,
2016; Chen, 2019; Mammen, Wilke, and Zapp, 2022). The advantage lies in
considerable gains of time and valid inference. I suspect the idea of breaking
the high-dimensionality of $k$-means by considering agglomerative approaches
in discrete econometric models to be very versatile and fruitful. Hence, this
paper can be seen as the first application and analysis of an agglomerative
clustering method to the econometric panel data model (1.1).333Since the first
arXiv version of this paper, Freeman and Weidner (2023) 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 (e.g. Bonhomme, Lamadon, and Manresa, 2022; Graf and
Luschgy, 2002), it seems not obvious when these conditions are met. The
measure of 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 (e.g., Zhang, Levina, and Zhu, 2017;
Lovász, 2012; Zeleneev, 2020; Auerbach, 2022). More generally, dyad, triad, or
tetrad comparisons have proven useful in a variety of other econometric
contexts (see, e.g., Honoré and Powell, 1994; Graham, 2017; Charbonneau, 2017;
Jochmans, 2017).

The rest of the paper is organized as follows. Section 2 introduces the three-
step estimation procedure. Section 3 presents large sample properties
including uniform consistency for the grouping structure and asymptotic
normality at parametric rates. Section 4 reports the results of a small-scale
Monte Carlo exercise calibrated to the empirical application of Section 5.
Section 6 concludes. All proofs are in the Appendix. Additional material is in
the Supplemental Material (Mugnier, 2022a), with section numbers S.1, etc.

#### Notation.

$\mathbf{1}\\{\cdot\\}$ denotes the indicator function. For any set
$S\subset\mathbb{R}^{p}$ (for any $p\geq 1$), let $S^{*}\equiv
S\backslash\\{0\\}$ and $|S|$ denote the cardinal of $S$. For all
$\mathcal{I}\subset\mathbb{N}^{*}$, all $k\in\mathbb{N}^{*}$, let
$\mathscr{P}_{k}(\mathcal{I})$ denote the set of subsets of $\mathcal{I}$ with
cardinal $k$. The operators $\overset{p}{\to}$, $\overset{d}{\to}$, and
$\mathrm{plim}$ denote respectively, convergence in probability, convergence
in distribution, and probability limit. All vectors are column vectors.
$\big{\|}\cdot\big{\|}$ denotes the Euclidean norm. For an $m\times n$ real
matrix $A$ of rank $r$, I write the transpose $A^{\prime}$, the Frobenius norm
as $\big{\|}A\big{\|}_{\textrm{F}}\equiv\sqrt{{\rm Tr}(A^{\prime}A)}$, the
spectral norm as $\big{\|}A\big{\|}_{\infty}\equiv\sigma_{1}(A)$, its max norm
as
$\big{\|}A\big{\|}_{\max}\equiv\max_{i=1,\ldots,m;j=1,\ldots,n}\left|a_{ij}\right|$,
and its nuclear norm as
$\big{\|}A\big{\|}_{1}\equiv\sum_{i=1}^{r}\sigma_{i}(A)$, where
$\sigma_{1}(A)\geq\cdots\geq\sigma_{r}(A)>0$ are the positive singular values
of $A$.

##  2 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^{\prime})^{\prime}\in\Theta$
of model (1.1), where

$\Theta\equiv\bigcup_{g\in\mathbb{N}^{*}}\Theta_{g},\Theta_{g}\equiv\left\\{g\right\\}\times\left\\{1,\ldots,g\right\\}^{N}\times\mathcal{A}^{gT}\times\mathcal{B}\text{
for some }\mathcal{B}\subset\mathbb{R}^{K},\mathcal{A}\subset\mathbb{R}.$

Sections 2.1-2.3 describe the three key ingredients of $\widehat{\theta}$: a
preliminary consistent estimator of $\beta$, a measure of pairwise distances
between units, and an agglomerative clustering algorithm. Section 2.4 formally
defines $\widehat{\theta}$.

###  2.1 Preliminary consistent estimation of the slope coefficient

The first ingredient is a preliminary consistent estimator
$\widehat{\beta}^{1}$ of $\beta$:

$\big{\|}\widehat{\beta}^{1}-\beta\big{\|}=o_{p}(1)\quad\text{ as
}\min(N,T)\to\infty.$

Intuitively, as covariates can be arbitrarily correlated with grouped fixed
effects, the clustering problem is simpler for the approximate “pure” grouped
fixed effects model:

$y_{it}-x_{it}^{\prime}\widehat{\beta}^{1}=\alpha_{g_{i}t}+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 provides various computationally
simple estimators under different identifying assumptions (e.g., Moon and
Weidner, 2019; Chernozhukov, Hansen, Liao, and Zhu, 2019; Beyhum and Gautier,
2019, 2023). Large sample properties developed in Section 3 do not depend on
the rate of convergence of the computationally simple preliminary estimator
which is generally less than $\sqrt{NT}$. Only an upper bound on this rate
must be known to implement the method, which is the case for all the above-
mentioned approaches under regularity conditions.

As an illustrative example, consider the nuclear norm regularized estimator.
Let $Y\equiv(y_{it})_{i=1,\ldots,N;t=1,\ldots,T}\in\mathbb{R}^{N\times T}$,
$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\\}$, and $v\cdot
X\equiv\sum_{k=1}^{K}X_{k}v_{k}$ for all $v\in\mathbb{R}^{K}$. For all
regularization parameter $\psi\in(0,\infty)$, let $Q_{\psi}$ denote the
nuclear norm regularized concentrated objective function such that, for all
$\beta\in\mathbb{R}^{K}$,

$Q_{\psi}(\beta)\equiv\min_{\Gamma\in\mathbb{R}^{N\times
T}}\left\\{\frac{1}{2NT}\big{\|}Y-\beta\cdot
X-\Gamma\big{\|}^{2}_{\textrm{F}}+\frac{\psi}{\sqrt{NT}}\big{\|}\Gamma\big{\|}_{1}\right\\}.$

Define

$\widehat{\beta}^{1}(\psi)\in\operatorname*{arg\,min}_{\beta\in\mathbb{R}^{K}}Q_{\psi}(\beta).$
(2.1)

Minimization problem (2.1) is convex and easily solvable with modern
optimization techniques. Under regularity conditions,
$\widehat{\beta}^{1}(\psi)$ is unique and called the nuclear norm regularized
estimator (see, e.g., Moon and Weidner, 2019). Instead of controlling for
unobserved discrete trends $\Gamma=[\alpha_{g_{i}t}]_{i,t}$, which leads to an
NP-hard problem, the nuclear norm regularized estimator penalizes the nuclear
norm (the sum of singular values) of an unrestricted matrix $\Gamma$ of
individual trends. Under some assumptions (especially on $\psi$), the
interactive fixed effects structure of grouped fixed effects is sufficient for
$\widehat{\beta}^{1}(\psi)$ to reach a convergence rate of at least
$\sqrt{\min(N,T)}$ as proved in Section A. In the Monte Carlo experiments, the
theoretically valid choice $\psi\equiv\log(\log(T))/\sqrt{16\min(N,T)}$ is
used.

If the researcher is willing to assume a known upper bound on the number of
groups and a factor structure on the covariates, other examples include the
correlated random effects estimators proposed in Chetverikov and Manresa
(2021) or Armstrong, Weidner, and Zeleneev (2022)’s approach (which only
requires an upper bound on $G$).

###  2.2 Pairwise distance between cross-sectional units

Equipped with a preliminary consistent estimator $\widehat{\beta}^{1}$, the
second ingredient is a measure of distance between any two units $i$ and $j$
informative of whether $i$ and $j$ belong to the same group or not. This can
be achieved by using the linear (or multiplicative) structure of the model.
Define the first-step residual $\widehat{v}_{it}\equiv
y_{it}-x_{it}^{\prime}\widehat{\beta}^{1}$. An empirical distance
$\widehat{d}_{\infty}^{2}(i,j)$ between $i$ and $j$ is

$\widehat{d}_{\infty}^{2}(i,j)\equiv\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|.$

The distance $\widehat{d}_{\infty}^{2}$ is borrowed from the statistical
literature on graphon estimation (see, e.g., Lovász, 2012; Zeleneev, 2020;
Auerbach, 2022). Other distances could be employed, e.g.,

$\widehat{d}_{2}^{2}(i,j)\equiv\frac{1}{N}\sum_{k=1}^{N}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|,$

and I focus on one that works well in practice and turns out to be convenient
to prove large sample properties (for instance, $\widehat{d}_{2}^{2}$ would
rule out asymptotically negligible groups to obtain consistent clustering
estimation while $\widehat{d}_{\infty}^{2}$ allows for such groups). I collect
all pairwise distances in the symmetric dissimilarity matrix
$\widehat{D}\equiv(\widehat{d}_{\infty}^{2}(i,j))_{(i,j)\in\\{1,\ldots,N\\}^{2}}$.

Why is the empirical distance $\widehat{d}^{2}_{\infty}(i,j)$ informative
about whether $g_{i}=g_{j}$? The intuition is as follows.444See also p.14 in
Zeleneev (2020) in a network setting. Because
$\widehat{\beta}^{1}\overset{p}{\to}\beta$, it holds with probability
approaching one that $\widehat{v}_{it}\approx\alpha_{g_{i}t}+v_{it}$. Under
weak time dependence, tails, and cross-sectional independence restrictions on
the error terms of the form $\mathbb{E}[(v_{it}-v_{jt})v_{kt}]=0$, it then
holds “uniformly” over $i$,$j$, and
$k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}$ that

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\approx\frac{1}{T}\sum_{t=1}^{T}(\alpha_{g_{i}t}-\alpha_{g_{j}t})\alpha_{g_{k}t}.$

If $g_{i}=g_{j}$, then $\alpha_{g_{i}t}-\alpha_{g_{j}t}=0$ and

$\displaystyle\widehat{d}_{\infty}^{2}(i,j)$
$\displaystyle=\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|$
$\displaystyle\approx\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}\underbrace{(\alpha_{g_{i}t}-\alpha_{g_{j}t})}_{=0}\alpha_{g_{k}t}\right|$
$\displaystyle=0.$

Reciprocally, if

$\widehat{d}_{\infty}^{2}(i,j)=\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|\approx
0,$ (2.2)

then necessarily $g_{i}=g_{j}$. In fact, suppose that $g_{i}\neq g_{j}$. Then,
provided each group has at least two units asymptotically (which is weak)
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 (2.2) therefore implies

$\displaystyle\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{k^{*}t}\approx\frac{1}{T}\sum_{t=1}^{T}(\alpha_{g_{i}t}-\alpha_{g_{j}t})\alpha_{g_{k^{*}}t}$
$\displaystyle\approx 0,$ (2.3)
$\displaystyle\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{l^{*}t}\approx\frac{1}{T}\sum_{t=1}^{T}(\alpha_{g_{i}t}-\alpha_{g_{j}t})\alpha_{g_{l^{*}}t}$
$\displaystyle\approx 0.$ (2.4)

Differencing (2.3)-(2.4) and using that $\alpha_{g_{k^{*}t}}=\alpha_{g_{i}t}$
and $\alpha_{g_{l^{*}t}}=\alpha_{g_{j}t}$ yields

$\frac{1}{T}\sum_{t=1}^{T}(\alpha_{g_{i}t}-\alpha_{g_{j}t})^{2}\approx 0,$

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 $c_{g,\widetilde{g}}>0$ such that

$\frac{1}{T}\sum_{t=1}^{T}(\alpha_{gt}-\alpha_{\widetilde{g}t})^{2}\geq
c_{g,\widetilde{g}}>0.$

Section 3 formalizes this equivalence result by proving uniform asymptotic
control of the remainders in the stochastic approximations.

###  2.3 Agglomerative clustering based on thresholding distances

Given a preliminary consistent estimator $\widehat{\beta}^{1}$ and a distance
$\widehat{d}_{\infty}^{2}$, the last ingredient is an agglomerative clustering
algorithm that builds clusters from the dissimilarity matrix
$\widehat{D}\equiv(\widehat{d}_{\infty}^{2}(i,j))_{(i,j)\in\\{1,\ldots,N\\}^{2}}$.
Section 3 (see also Lemma C.1) provides sufficient conditions under which any
matrix $\widehat{W}(c_{NT})\equiv\mathbf{1}\\{\widehat{D}\leq c_{NT}\\}$ with
$c_{NT}\to 0$ converges in max norm to
$W\equiv(\mathbf{1}\\{g_{i}=g_{j}\\})_{(i,j)\in\\{1,\ldots,N\\}^{2}}$:

$\big{\|}\widehat{W}(c_{NT})-W\big{\|}_{\max}=\max_{(i,j)\in\\{1,\ldots,N\\}^{2}}\left|\widehat{W}_{ij}(c_{NT})-W_{ij}\right|=o_{p}(1)\quad\text{
as }\min(N,T)\to\infty.$

In light of this result solving the problem of estimating a sparse (block
diagonal) graphon, it is natural to propose agglomerative clustering
algorithms based on merging units whose (weighted) pairwise distances fall
below some threshold. I focus on the following intuitive and simple algorithm.

Agglomerative Clustering Algorithm:
Set the value of willingness-to-pair parameter $c_{NT}\in(0,\infty)$. Then,

  1. 1.

If there are no unclustered pairs $(i,j)$ such that
$\widehat{d}_{\infty}^{2}(i,j)\leq c_{NT}$, assign each unclustered units to
singleton groups.

  2. 2.

Else:

     * 2.1.

Create a new group consisting of the pair(s) of units with the smallest
pairwise distance. Let $k\geq 2$ denote the number of such reference units.

     * 2.2.

Add to this group all units whose weighted pairwise distance to those $k$
reference units falls below $c_{NT}$.

     * 2.3.

Go to Step 1.

This procedure always produces $\widehat{G}\in\\{1,\ldots,N\\}$ non-empty
groups.

###  2.4 A three-step triad pairwise-differencing estimator

The TPWD estimator $\widehat{\theta}\in\Theta$ for $\theta\in\Theta$ is
obtained as follows.

  1. 1.

Preliminary Slope Estimation Step: Compute $\widehat{\beta}^{1}$ as in Section
2.1.

  2. 2.

Agglomerative Clustering Step:

    1. 2.a.

Set $c_{NT}\in(0,\infty)$.

    2. 2.b.

Compute $\widehat{D}$ as in Section 2.2.

    3. 2.c.

Compute $\left\\{\widehat{g}_{1},\ldots,\widehat{g}_{N}\right\\}$ and
$\widehat{G}=\left|\left\\{\widehat{g}_{1},\ldots,\widehat{g}_{N}\right\\}\right|$
as in Section 2.3.

  3. 3.

Projection Step:
Compute:

$\left(\widehat{\beta}^{\prime},\widehat{\alpha}_{11},\ldots,\widehat{\alpha}_{\widehat{G}T}\right)\in\operatorname*{arg\,max}_{\left(\beta^{\prime},\alpha_{11},\ldots,\alpha_{\widehat{G}T}\right)\in\mathcal{B}\times\mathcal{A}^{\widehat{G}T}}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(y_{it}-x_{it}^{\prime}\beta-\alpha_{\widehat{g}_{i}t}\right)^{2}.$

Provided $\mathcal{B}=\mathbb{R}^{K}$ and $\mathcal{A}=\mathbb{R}$, the
projection step to obtain
$\left(\widehat{\beta}^{\prime},\widehat{\alpha}_{11},\ldots,\widehat{\alpha}_{\widehat{G}T}\right)^{\prime}$
is a pooled OLS regression of $y_{it}$ on $x_{it}$ and the interactions of
estimated groups and time dummies.

Remark 1: Consider a finite sample of fixed dimensions $N$ and $T$. If all
random variables except group memberships are continuous, then to the extreme
where $c_{NT}\to 0$, $\widehat{G}\to N$ and each group contains a single unit.
To the extreme where $c_{NT}\to+\infty$, $\widehat{G}\to 1$ and a single group
contains all units. Given the low CPU time of the method, an entire
regularization path can be reported by the researcher by making $c_{NT}$ vary
between these two regimes. In particular, the clustering step can be made all
vectorized, greatly reducing computational burden compared to running
loops.555MATLAB code is provided on the author’s website:
https://martinmugnier.github.io/research.

Remark 2: In Section 3, I provide theory for the choice of $c_{NT}$. In
Section 3.3, I provide a simple data-driven selection rule that works well in
finite samples. It is possible to further improve the finite sample
performance by re-running the first step with
$\widehat{v}_{it}=y_{it}-x_{it}^{\prime}\widehat{\beta}$ in place of
$\widehat{v}_{it}=y_{it}-x_{it}^{\prime}\widehat{\beta}^{1}$ to obtain new
$\widehat{g}_{1},\ldots,\widehat{g}_{N}$, and then re-running the second and
third steps and iterating again until some convergence criterion is achieved.
The asymptotic results hold for any of the subsequent iterates and Monte Carlo
performance suggest important gains of precision.

Remark 3: The computation routine requires $O(N^{3}T)$ operations, and it is
quite fast compared to the NP-hard $k$-means problem that forms the basis of
Bonhomme and Manresa (2015)’s GFE estimator. Whether the computational cost of
an exact algorithm can be lowered to $O(N^{2}T)$ remains to the best of my
knowledge an open problem. A current limitation of the method is rather the
memory required to store triad differences of residuals, but it seems hardly
possible to cluster units without any measure of distance between them. When
unobserved heterogeneity is assumed to be time-constant, the $O(N^{3}T)$
computation cost can be reduced to $O(N^{2}T)$ and the preliminary estimator
can be replaced with any standard differencing fixed effects estimator such as
Arellano and Bond (1991). See also Wooldridge (2010). This particular case, as
well as some extensions of model (1.1), is discussed in the Supplemental
Material S.2.

##  3 Large sample properties

In this section, I derive the large sample properties of the TPWD estimator
introduced in Section 2 under conditions similar to the existing literature.
Consider the data generating process:

$y_{it}=x_{it}^{\prime}\beta^{0}+\alpha_{g_{i}^{0}t}^{0}+v_{it},\quad
i=1,\ldots,N,\;t=1,\ldots,T,$ (3.1)

where $g_{i}^{0}\in\left\\{1,\ldots,G^{0}\right\\}$ denotes group membership,
and where the $0$ superscripts refer to true parameter values. I consider
asymptotic sequences $(N,T)$ where $N$ and $T$ diverge jointly to infinity
which I denote $\min(N,T)\to\infty$. The number of groups $G^{0}$ is fixed
relative to $N,T$ but unknown, and I defer the discussion on the case of an
increasing sequence $G^{0}=G^{0}_{NT}$ to the Supplemental Material S.1.

###  3.1 Consistency of estimated group memberships

Consider the following assumptions.

###### Assumption 1

$\big{\|}\widehat{\beta}^{1}-\beta^{0}\big{\|}=O_{p}(r_{NT})$ for some
deterministic sequence $r_{NT}\to 0$ as $N$ and $T$ tend to infinity.

###### Assumption 2

There exist constants $(C,\nu,\kappa)^{\prime}\in(0,+\infty)^{2}\times(0,1/2)$
such that

  1. (a)

$NT^{-\nu}\to 0$ as $\min(N,T)\to\infty$.

  2. (b)

$c_{NT}\overset{p}{\to}0$, $c_{NT}r_{NT}^{-1}\overset{p}{\to}\infty$, and
$\Pr(c_{NT}T^{\kappa}\geq C)\to 1$ as $\min(N,T)\to\infty$.

###### Assumption 3

There exist constants $a,b,c,d_{1},d_{2},M>0$ and a sequence $\tau(t)\leq
e^{-at^{d_{1}}}$ such that:

  * (a)

$\mathcal{A}$ is a compact subset of $\mathbb{R}$.

  * (b)

For all
$(i,t)\in\left\\{1,\ldots,N\right\\}\times\left\\{1,\ldots,T\right\\}$:
$\Pr(\left|v_{it}\right|>m)\leq e^{1-(m/b)^{d_{2}}}$ for all $m>0$.

  * (c)

For all $(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}$ such that
$g\neq\widetilde{g}$:

$\mathrm{plim}_{\min(N,T)\to\infty}T^{-1}\sum_{t=1}^{T}(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})^{2}=c_{g,\widetilde{g}}\geq
c.$

  * (d)

For all
$(i,j,k,g,\widetilde{g})\in\mathscr{P}_{3}(\left\\{1,\ldots,N\right\\})\times\left\\{1,\ldots,G^{0}\right\\}^{2}$
such that $g\neq\widetilde{g}$, $\left\\{v_{it}\right\\}_{t}$,
$\left\\{(v_{it}-v_{jt})v_{kt}\right\\}_{t}$,
$\left\\{\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right\\}_{t}$, and
$\left\\{(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})v_{it}\right\\}_{t}$ are
strongly mixing processes with mixing coefficients $\tau(t)$. Moreover,
$\mathbb{E}[v_{it}v_{jt}]=0$.

  * (e)

$\lim_{\min(N,T)\to\infty}\Pr(\min_{g\in\left\\{1,\ldots,G^{0}\right\\}}\;\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}\geq
2)=1$.

  * (f)

As $N$ and $T$ tend to infinity:

$\sup_{i\in\left\\{1,\ldots,N\right\\}}\;\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}\geq
M\right)=o(T^{-\delta})\text{ for all }\delta>0.$

Assumption 1 requires $\widehat{\beta}^{1}$ to be consistent for $\beta^{0}$
at a rate bounded by $r_{NT}$. This rate can be slow. Examples of
computationally simple estimators satisfying this condition are given in
Sections 2.1 and A. Assumption 2 allows $T$ to grow considerably more slowly
than $N$ (when $\nu\gg 1$). It requires that the tuning parameter decreases to
zero, but not too fast, at a rate bounded below by $T^{-1/2}$ and strictly
slower than $r_{NT}$. I use probability limit because I allow $c_{NT}$ to be
data-driven and hence random. Assumptions 3(a)-(b) and 3(d) collect standard
moment, tail, and dependence conditions. They do not impose homoscedasticity
but only require uniform bounds on the unconditional variances. Assumption
3(c) requires groups to be well-separated. Assumprion 3(d) rules out cross-
sectional correlation in the error term. Assumption 3(e) allows for
asymptotically negligible groups but requires that each group has at least two
members with probability approaching one. Assumption 3(f) is a slight
reinforcement of Bonhomme and Manresa (2015)’s Assumption 2(e). It holds if
covariates have bounded support or if they satisfy dependence and tail
conditions similar to $v_{it}$. All results below are understood up to group
relabeling.

######  Proposition 3.1 (Sup-Norm Classification Consistency)

Let Assumptions 1-3 hold. Then, as $N$ and $T$ tend to infinity,

$\max_{i\in\left\\{1,\ldots,N\right\\}}\left|\widehat{g}_{i}-g_{i}^{0}\right|\overset{p}{\to}0,$
(3.2)

and

$\widehat{G}-G^{0}\overset{p}{\to}0.$ (3.3)

###  3.2 Asymptotic distribution

The next assumption is useful to establish the asymptotic distribution of
$\widehat{\beta}$ and $\widehat{\alpha}_{gt}$.

###### Assumption 4

  * (a)

For all $g\in\left\\{1,\ldots,G^{0}\right\\}$:
$\mathrm{plim}_{\min(N,T)\to\infty}\frac{1}{N}\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}=\pi_{g}>0$.

  * (b)

For all
$(i,j,t)\in\left\\{1,\ldots,N\right\\}^{2}\times\left\\{1,\ldots,T\right\\}$:
$\mathbb{E}[x_{jt}v_{it}]=0$.

  * (c)

There exist positive definite matrices $\Sigma_{\beta}$ and $\Omega_{\beta}$
such that

$\displaystyle\Sigma_{\beta}$
$\displaystyle=\underset{\min(N,T)\to\infty}{\mathrm{plim}}\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}(x_{it}-\overline{x}_{g_{i}^{0}t})(x_{it}-\overline{x}_{g_{i}^{0}t})^{\prime},$
$\displaystyle\Omega_{\beta}$
$\displaystyle=\underset{\min(N,T)\to\infty}{\mathrm{plim}}\frac{1}{NT}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T}\sum_{s=1}^{T}\mathbb{E}\left[v_{it}v_{js}(x_{it}-\overline{x}_{g_{i}^{0}t})(x_{js}-\overline{x}_{g_{j}^{0}s})^{\prime}\right],$

where
$\overline{x}_{gt}\equiv\left(\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}\right)^{-1}\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}x_{it}$.

  * (d)

As $N$ and $T$ tend to infinity:
$\frac{1}{\sqrt{NT}}\sum_{i=1}^{N}\sum_{t=1}^{T}(x_{it}-\overline{x}_{g_{i}^{0}t})v_{it}\overset{d}{\to}\mathcal{N}(0,\Omega_{\beta})$.

  * (e)

For all
$(g,t)\in\left\\{1,\ldots,G^{0}\right\\}\times\left\\{1,\ldots,T\right\\}$:

$\lim_{\min(N,T)\to\infty}\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathbb{E}\left[\mathbf{1}\\{g_{i}^{0}=g\\}\mathbf{1}\\{g_{j}^{0}=g\\}v_{it}v_{jt}\right]=\omega_{gt}>0.$

  * (f)

For all
$(g,t)\in\left\\{1,\ldots,G^{0}\right\\}\times\left\\{1,\ldots,T\right\\}$: as
$N$ and $T$ tend to infinity,

$\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}v_{it}\overset{d}{\to}\mathcal{N}(0,\omega_{gt}).$

Assumption 4 ensures that the infeasible least squares estimator has a
standard asymptotic distribution. Assumption 4(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}$.

######  Corollary 3.2 (Asymptotic Distribution)

Let Assumptions 1-4 hold. Then, as $N$ and $T$ tend to infinity,

$\sqrt{NT}(\widehat{\beta}-\beta^{0})\overset{d}{\to}\mathcal{N}\left(0,\Sigma_{\beta}^{-1}\Omega_{\beta}\Sigma_{\beta}^{-1}\right),$
(3.4)

and, for all $t$:

$\sqrt{N}(\widehat{\alpha}_{gt}-\alpha_{gt}^{0})\overset{d}{\to}\mathcal{N}\left(0,\frac{\omega_{gt}}{\pi_{g}^{2}}\right),\;g=1,\ldots,G^{0},$
(3.5)

where $\Sigma_{\beta},\Omega_{\beta},\omega_{gt}$, and $\pi_{g}$ are defined
in Assumption 4.

Consistent plug-in estimates of the asymptotic variances can easily be
constructed (see, e.g., Supplemental Material in Bonhomme and Manresa, 2015).

###  3.3 Choice of the tuning parameter

Assumption 2 provides theoretical guidance for choosing the tuning parameter
$c_{NT}$. Because this guidance is based on asymptotic arguments, it has the
convenience of not requiring any unknown constant but the drawback of
depending on rates. As is the case with bandwidth selection in nonparametric
density estimation, infinitely many different values consistent with the
theory can be chosen in applications where only a finite sample is observed.

Below, I provide a theoretically valid data-driven selection rule that is
shown to work well on the Monte Carlo experiments reported in Section 4. I
assume that a consistent estimator $\widehat{\beta}^{1}$ that converges at a
rate of at least $\sqrt{\min(N,T)}$ is available. For instance, by Proposition
A.1, the nuclear norm regularized estimator $\widehat{\beta}^{1}(\psi)$ with
$\psi=\log(\log(T))/\sqrt{16\min(N,T)}$ verifies this assumption under weak
conditions. The data-driven rule is as follows:

  1. 1.

Compute $\widehat{v}_{it}\equiv y_{it}-x_{it}^{\prime}\widehat{\beta}^{1}$.

  2. 2.

Compute
$\widehat{\sigma}^{1}\equiv\sqrt{\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(\widehat{v}_{it}-\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\widehat{v}_{it}\right)^{2}}$.

  3. 3.

Set $c_{NT}\equiv\widehat{\sigma}^{1}\log(T)/\sqrt{T}$.

Theoretical validity follows from $\widehat{\sigma}^{1}=O_{p}(1)$ under the
maintained assumptions, so that $c_{NT}$ verifies Assumption 2.

##  4 Monte Carlo simulations

In this section, I investigate the finite sample properties of the TPWD across
two Monte Carlo experiments calibrated to the empirical application: a pure
GFE model without covariates, and a full GFE model with a scalar covariate. I
consider sample sizes
$(N,T)\in\left\\{90,180\right\\}\times\left\\{7,10,20,40\right\\}$ and numbers
of groups $G\in\left\\{3,4\right\\}$. Mimicking the empirical results of
Bonhomme and Manresa (2015), I let $\alpha_{1t}\equiv 1$,
$\alpha_{2t}\equiv\frac{t-1}{T-1}$, $\alpha_{3t}\equiv 0$, and
$\alpha_{4t}\equiv\mathbf{1}\\{t\geq\lfloor T/2\rfloor\\}\frac{t-\lfloor
T/2\rfloor}{T-\lfloor T/2\rfloor}$. In their application, the first group is
referred to as “high democracy”, the second to “transition”, the third to “low
democracy”, and the fourth one to “late transition”. I consider balanced
groups:

$g_{i}=1+\sum_{g=1}^{G-1}\mathbf{1}\\{i>g\lfloor N/G\rfloor\\},\quad
i=1,\ldots,N,$

and refer to Section S.3.1 for similar results with unbalanced groups.
Similarly, I refer to Sections S.3.2-3.6 for additional results for models
with more groups, lagged outcomes, unit-specific effects, higher signal-to-
noise ratio, or time-invariant unobserved heterogeneity. For all experiments,
results are averaged across $500$ Monte Carlo samples.

###  4.1 Pure GFE model

First, consider a pure GFE model without covariates:

$y_{it}=\alpha_{g_{i}t}+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. The signal-to-noise ratio is one. I evaluate the performance of
the TPWD estimator with $\widehat{\beta}^{1}=0$ and
$c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$ by considering the following
three criteria. First, I report the (average) root mean square errors (RMSE)
of the estimated grouped fixed effects, which is measured by

$\displaystyle{\rm RMSE}(\widehat{\alpha})$
$\displaystyle\equiv\frac{1}{500}\sum_{b=1}^{500}\sqrt{\frac{1}{NT}\sum_{i=1}^{N}\sum_{i=t}^{T}(\widehat{\alpha}_{\widehat{g}_{i}t}-\alpha_{g_{i}t})^{2}}.$

Second, I report the average estimated number of group $\widehat{G}$. Third, I
assess clustering accuracy using three metrics from the statistical learning
literature on binary classification. Define

$\displaystyle FP$
$\displaystyle\equiv\sum_{i<j}\mathbf{1}\\{\widehat{g}_{i}=\widehat{g}_{j}\\}\mathbf{1}\\{g_{i}\neq
g_{j}\\},$ $\displaystyle TP$
$\displaystyle\equiv\sum_{i<j}\mathbf{1}\\{\widehat{g}_{i}=\widehat{g}_{j}\\}\mathbf{1}\\{g_{i}=g_{j}\\},$
$\displaystyle FN$
$\displaystyle\equiv\sum_{i<j}\mathbf{1}\\{\widehat{g}_{i}\neq\widehat{g}_{j}\\}\mathbf{1}\\{g_{i}=g_{j}\\},$
$\displaystyle TN$
$\displaystyle\equiv\sum_{i<j}\mathbf{1}\\{\widehat{g}_{i}\neq\widehat{g}_{j}\\}\mathbf{1}\\{g_{i}\neq
g_{j}\\}.$

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

$\displaystyle P$ $\displaystyle\equiv\frac{TP}{TP+FP},$ $\displaystyle R$
$\displaystyle\equiv\frac{TP}{TP+FN},$ $\displaystyle RI$
$\displaystyle\equiv\frac{TP+TN}{TP+FP+FN+TN}.$

The precision rate (P) measures the proportion of correctly matched pairs of
units among all matched pairs of units. The recall rate (R) measures the
proportion of correctly matched pairs of units among all pairs of units that
belong to the same population cluster. The Rand Index (RI) measures the
proportion of correct matching/non-matching decisions among all decisions made
by the clustering algorithm.

Table 1 shows the results for the first two criteria and compares them to a
Lloyd’s approximation of the GFE estimator with user-specified number of
groups $g\in\left\\{2,3,10\right\\}$ and $10,000$ initializers. TPWD has the
lowest RMSE in all settings for which $T\geq 20$. For settings with $T<20$,
TPWD does always better than mispecified GFE estimators and has a RMSE twice
as large as the well-specified GFE. The number of groups is estimated
consistently in less than 0.4s when $N=90$, and less than 3s when $N=180$.
TPWD estimates are exact, which is not the case of the GFE estimator
approximated by a Lloyd’s type algorithm (see MATLAB documentation on
$k$-means).

Table 1: Estimation of the grouped fixed effects (pure GFE model)

[h]  TPWD GFEG=2 GFEG=3 GFEG=10 $G$ $N$ $T$ $\widehat{G}$ RMSE CPU time RMSE
RMSE RMSE 3 90 7 4.486 0.154 0.066 0.272 0.096 0.206 10 3.784 0.123 0.090
0.264 0.084 0.180 20 3.086 0.074 0.164 0.250 0.087 0.150 40 3.002 0.061 0.331
0.244 0.088 0.134 3 180 7 5.450 0.145 0.521 0.271 0.077 0.189 10 4.270 0.106
0.683 0.261 0.059 0.160 20 3.246 0.059 1.249 0.246 0.061 0.127 40 3.008 0.043
2.427 0.241 0.073 0.107 4 90 7 4.744 0.176 0.065 0.265 0.128 0.215 10 3.946
0.176 0.085 0.250 0.131 0.192 20 3.576 0.165 0.173 0.236 0.124 0.157 40 3.652
0.156 0.352 0.233 0.126 0.135 4 180 7 5.502 0.145 0.559 0.205 0.124 0.208 10
4.434 0.139 0.742 0.193 0.119 0.184 20 3.952 0.133 1.409 0.191 0.105 0.141 40
3.978 0.122 2.764 0.192 0.113 0.112 • Notes: This table reports the estimated
number of group ($\widehat{G}$), the root mean square error (RMSE) and the
execution time in seconds (CPU time) for the triad pairwise distance estimator
(TPWD) with $\widehat{\beta}^{1}=0$ and
$c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$. It reports the RMSE for a
Lloyd’s approximation to the grouped fixed effects estimator with user-
specified number of groups $g\in\left\\{2,3,10\right\\}$ and 10,000
initializers (GFEG=g). Results averaged across 500 Monte Carlo samples.

Table 2 shows the result for the last three criteria. The accuracy of TPWD
classification in terms of Precision, Recall and Rand Index quickly improves
with $T$, homogeneously across the three measures. For $G=3$, the TPWD
estimator performs always better in the limit than the GFE estimator,
irrespective of the specified number of groups. This may suggest that Lloyd’s
algorithm approximating the GFE estimator does not return the global optimum.
For $G=4$ the performance of TPWD in terms of RI is close to the well-
specified GFE estimator. When $g=10$, the recall of the GFE estimator is quite
low (less than 50%); when $g=2$ the precision of the GFE estimator is quite
low (always less than $60\%$). In contrast, the TPWD estimator achieves a
data-driven balance between Precision and Recall.

Table 2: Classification accuracy (pure GFE model)

[h]  TPWD GFEG=2 GFEG=3 GFEG=10 $G$ $N$ $T$ P R RI P R RI P R RI P R RI 3 90 7
0.892 0.834 0.913 0.568 0.911 0.746 0.957 0.967 0.973 0.964 0.364 0.788 10
0.930 0.910 0.948 0.574 0.932 0.753 0.958 0.980 0.974 0.988 0.393 0.801 20
0.987 0.986 0.991 0.586 0.977 0.768 0.940 0.979 0.963 0.999 0.420 0.811 40
1.000 1.000 1.000 0.590 0.994 0.773 0.935 0.980 0.961 1.000 0.433 0.815 3 180
7 0.906 0.812 0.910 0.571 0.905 0.745 0.974 0.975 0.983 0.968 0.355 0.784 10
0.946 0.912 0.954 0.579 0.933 0.754 0.979 0.988 0.987 0.990 0.373 0.792 20
0.988 0.985 0.991 0.591 0.982 0.771 0.961 0.986 0.976 1.000 0.408 0.805 40
1.000 1.000 1.000 0.595 0.998 0.776 0.935 0.979 0.960 1.000 0.436 0.814 4 90 7
0.568 0.792 0.803 0.384 0.931 0.621 0.616 0.946 0.841 0.733 0.366 0.814 10
0.563 0.823 0.800 0.387 0.953 0.623 0.621 0.959 0.843 0.802 0.411 0.833 20
0.578 0.856 0.810 0.390 0.981 0.623 0.627 0.985 0.848 0.905 0.484 0.8630 40
0.604 0.853 0.827 0.398 0.992 0.632 0.625 0.991 0.847 0.968 0.538 0.884 4 180
7 0.716 0.827 0.779 0.564 0.964 0.656 0.724 0.907 0.806 0.775 0.206 0.622 10
0.714 0.908 0.798 0.569 0.985 0.663 0.726 0.929 0.814 0.832 0.227 0.637 20
0.706 0.950 0.802 0.565 0.991 0.658 0.731 0.983 0.831 0.929 0.268 0.667 40
0.744 0.960 0.835 0.565 0.999 0.658 0.724 0.981 0.824 0.988 0.315 0.695 •
Notes: This table reports the precision rate (P), recall rate (R), and Rand
Index (RI) for the triad pairwise distance estimator (TPWD) with
$\widehat{\beta}^{1}=0$ and $c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$, as
well as for a Lloyd’s approximation to the grouped fixed effects estimator
with user-specified number of groups $g\in\left\\{2,3,10\right\\}$ and 10,000
initializers (GFEG=g). Results are averaged across 500 Monte Carlo samples.

###  4.2 Full GFE model

In this section, I turn to a full GFE model with a scalar covariate:

$y_{it}=x_{it}\beta+\alpha_{g_{i}t}+v_{it},\quad i=1,\ldots,N,\;t=1,\ldots,T,$

where $\beta=1$, $x_{it}=0.5\alpha_{g_{i}t}+u_{it}$,
$u_{it}\overset{iid}{\sim}\mathcal{N}(0,(1/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. The signal-to-noise
ratio is one.

I evaluate the performance of the TPWD estimator using
$\widehat{\beta}^{1}(\psi)$ the nuclear-norm regularization described in
Section 2.1 with $\psi=\log(\log(T))/\sqrt{16\min(N,T)}$ as a preliminary
estimator and the data-driven rule
$c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$. I report the TPWD estimate
obtained after four iterations, where convergence is achieved. Beyond the root
mean square error (RMSE) of the estimated grouped fixed effects, I report that
of the estimated regression coefficient, which is measured by

$\displaystyle{\rm RMSE}(\widehat{\beta})$
$\displaystyle=\sqrt{\frac{1}{500}\sum_{b=1}^{500}\big{\|}\widehat{\beta}^{(b)}-\beta\big{\|}^{2}},$

and the coverage rate of a 95$\%$-level confidence interval for $\beta$ based
on an estimation of the large-$T$ asymptotic variance of the TPWD estimator
clustered at the unit level. I compare the TPWD estimator to the nuclear-norm
regularized estimator (NN-reg) with regularization parameter
$\psi=\log(\log(T))\min(N,T)^{-1/2}/4$ and the infeasible pooled OLS
regression with known group memberships (Infeasible OLS).

Table 3 shows the results in terms of bias, RMSE, and coverage of a
95$\%$-level confidence interval (when available). For $G=3$, the bias of the
TPWD is estimator is only $10\%$ to $0.1\%$ of that of NN-reg, which converges
very slowly but do not prevent TPWD from achieving high precision gains.
Similarly, the TPWD coverage rate reaches its specified level for $T=40$ and
is already quite close for $T=20$. The clustering task becomes more difficult
with $G=4$, and the price to pay for model selection is that bias exceeds
variance in finite sample so that $95\%$-level confidence intervals have poor
finite sample coverage (less than 50%). Under additional assumptions on the
covariate structure and knowledge of the number of groups, one solution to
obtain improved CI is to improve the rate of convergence of the preliminary
estimator using the (post-)spectral estimator (e.g. Chetverikov and Manresa,
2021). If these assumptions are overly restrictive, and given the low bias,
another possibility is to use the bootstrap.

Table 3: Estimation of the slope coefficient and grouped fixed effects (full
GFE model)

[h]  NN-reg TPWD Infeasible OLS $G$ $N$ $T$ Bias $\widehat{\beta}$ RMSE
$\widehat{\beta}$ Bias $\widehat{\beta}$ RMSE $\widehat{\beta}$ $.95$ RMSE
$\widehat{\alpha}_{gt}$ $\widehat{G}$ Bias $\widehat{\beta}$ RMSE
$\widehat{\beta}$ $.95$ RMSE $\widehat{\alpha}_{gt}$ 3 90 7 0.550 0.561 0.084
0.123 0.560 0.178 4.494 -0.001 0.002 0.953 0.061 10 0.448 0.451 0.060 0.095
0.638 0.144 3.718 0.000 0.001 0.971 0.061 20 0.306 0.313 0.013 0.045 0.878
0.083 3.102 -0.001 0.001 0.972 0.061 40 0.234 0.234 -0.001 0.000 0.944 0.061
3.010 -0.001 0.000 0.972 0.061 3 180 7 0.551 0.557 0.068 0.089 0.466 0.160
5.360 0.002 0.001 0.965 0.044 10 0.443 0.449 0.045 0.071 0.638 0.123 4.244
0.001 0.001 0.972 0.043 20 0.299 0.300 0.008 0.032 0.894 0.060 3.206 0.000
0.000 0.976 0.043 40 0.219 0.226 0.000 0.000 0.952 0.043 3.006 0.000 0.000
0.976 0.043 4 90 7 0.552 0.563 0.112 0.138 0.396 0.197 4.620 -0.002 0.002
0.950 0.071 10 0.450 0.452 0.121 0.138 0.258 0.195 3.772 0.000 0.001 0.966
0.070 20 0.308 0.315 0.109 0.123 0.144 0.178 3.454 -0.001 0.001 0.964 0.070 40
0.235 0.237 0.097 0.105 0.050 0.168 3.552 -0.001 0.000 0.974 0.071 4 180 7
0.498 0.504 0.059 0.071 0.506 0.155 5.514 0.002 0.001 0.962 0.051 10 0.397
0.404 0.067 0.078 0.352 0.148 4.388 0.001 0.001 0.968 0.049 20 0.267 0.268
0.068 0.071 0.126 0.141 3.930 0.000 0.000 0.976 0.050 40 0.195 0.202 0.059
0.063 0.088 0.130 3.962 0.000 0.000 0.977 0.050 • Notes: This table reports
results for the nuclear norm regularized (NN-reg) estimator with
regularization parameter $\psi=\log(\log(T))\min(N,T)^{-1/2}/4$; triad
pairwise differencing estimator (TPWD) obtained after four iterations starting
at NN-reg and setting $c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$; and
infeasible pooled OLS regression with known group memberships (Infeasible
OLS). .95 denotes coverage rate of a 95%-level confidence interval for $\beta$
based on large $N,T$ approximations and an estimator of the asymptotic
variance clustered at the unit level. Results are averaged across 500 Monte
Carlo samples.

Table 4 reports the results regarding clustering accuracy. For $G=3$, the
three measures are little affected by the slow rate of convergence of the
preliminary consistent estimator, even for small values of $T$. For $G=4$,
Precision decreases significantly (between 0.53 and 0.72) while Recall remains
at a high level (between 0.77 and 0.95). This mitigates the drop in the Rand
Index.

Table 4: Classification accuracy of the TPWD estimator (full GFE model)

[h]  $G$ $N$ $T$ Precision Recall Rand Index 3 90 7 0.836 0.798 0.881 10 0.886
0.884 0.923 20 0.974 0.975 0.983 40 1.000 1.000 1.000 3 180 7 0.872 0.796
0.894 10 0.916 0.894 0.938 20 0.987 0.984 0.990 40 1.000 1.000 1.000 4 90 7
0.533 0.771 0.779 10 0.529 0.812 0.775 20 0.549 0.869 0.788 40 0.571 0.867
0.803 4 180 7 0.703 0.816 0.767 10 0.696 0.901 0.782 20 0.688 0.947 0.785 40
0.720 0.953 0.815 • Notes: This table reports results for the triad pairwise
differencing estimator (TPWD) obtained after four iterations starting at the
nuclear norm regularized estimator $\widehat{\beta}^{1}(\psi)$ with
$\psi=\log(\log(T))\min(N,T)^{-1/2}/4$ and setting
$c_{NT}=\widehat{\sigma}^{1}\log(T)/\sqrt{T}$. Results are averaged across 500
Monte Carlo samples.

##  5 Empirical application: income and (waves of) democracy

Characterizing the statistical relationship between income and democracy has
been a longstanding problem in political science and economics (Lipset, 1959;
Barro, 1999). Using panel data for $N=90$ countries observed at $T=7$ points
in time over the period 1970-2000, Acemoglu, Johnson, Robinson, and Yared
(2008) documented that the statistically significant positive effect of income
on democracy disappears when including country fixed-effects 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. Bonhomme and Manresa (2015)
proposed to test this assumption by using an approximation of the GFE
estimator using alternative minimization and reporting the results for several
numbers of groups. In particular, they considered a BIC criterion to perform
model selection, unfortunately not shown to be consistent for the number of
groups when $T$ grows slowly with $N$. They argued that the true number of
groups should be less than $10$, reporting statistically significant income
effects between $0.061$ and $0.089$ (elasticity of a measure of democracy to
lagged income per capita).

In this section, I revisit their study and consistently estimate 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_{i}t}$:

$democracy_{it}=\beta_{1}democracy_{it-1}+\beta_{2}logGDPpc_{it-1}+\alpha_{g_{i}t}+v_{it}.$

The data is obtained from the balanced subsample of Acemoglu, Johnson,
Robinson, and Yared (2008).666Available here:
https://www.aeaweb.org/articles?id=10.1257/aer.98.3.808. The preliminary
estimator is a nuclear norm regularized (NN-reg) estimator with tuning
parameter set to the theoretically valid rule
$\psi=\log(\log(T))/\sqrt{16\min(N,T)}$ (results are not sensitive to
alternative choices such as 3-fold cross-validation as long as the grid is
sufficiently large and fine). The data-driven rule of the TPWD estimator is as
in Section 3.3.

Figure 1 shows the regularization path of the TPWD estimator at the first four
iterations: the estimated number of groups as a function of different values
of the thresholding parameter (data-driven choices appear as vertical lines).
Estimates values stabilize after iteration 9, and then alternate between the
results of iterations 8 and 9 where only one unit switches from Group 1 to
Group 2. The number of estimated groups moves from $2$ to $3$ at the third
iteration and remains at $\widehat{G}=3$ for iterations 8 and 9.

Figure 1: TPWD regularization path

[h]

  * •

Notes: The solid blue, dashed green, dotted red, and dash-dotted magenta lines
correspond to the 1st, 2nd, 3rd, and 4th iterations of the TPWD estimator with
its thresholding parameter varying over a grid that contains $1,000$
equispaced points between 0 and 0.176. Vertical lines represent the data-
driven rule $c_{NT}=\widehat{\sigma}^{1}\log(7)/\sqrt{7}$ at each iteration.

Table 5 reports and compares NN-reg estimates, TPWDkit estimates at iterations
$k\in\left\\{2,3,8,9\right\\}$ (results vary smoothly in between), and GFEG=g
estimates with user-specified number of groups $g\in\left\\{2,3,10\right\\}$.

Table 5: Estimation of the effect of lagged democracy and lagged income on
democracy

[h]  NN-reg TPWD1it TPWD2it TPWD8it TPWD9it GFEG=2 GFEG=3 GFEG=10 Dependent
variable: Democracy $\widehat{G}$ - 2 2 3 3 - - - Lagged Democracy
$(\beta_{1})$ 0.800 0.691 0.679 0.580 0.570 0.601 0.407 0.277 (0.044) (0.046)
(0.050) (0.051) (0.041) (0.052) (0.049) Lagged Income $(\beta_{2})$ 0.016
0.078 0.080 0.102 0.105 0.061 0.089 0.075 (0.013) (0.013) (0.013) (0.014)
(0.011) (0.011) (0.008) Cumulative Income
$\left(\frac{\beta_{2}}{1-\beta_{1}}\right)$ 0.078 0.252 0.248 0.244 0.243
0.152 0.151 0.104 (0.019) (0.019) (0.014) (0.014) (0.021) (0.013) (0.009) •
Notes: Balanced sample of Acemoglu, Johnson, Robinson, and Yared (2008) with
$N=90$ countries and $T=7$ time periods at the five-year frequency between
1970-2000. Democracy is measured as the Freedom House indicator of democracy.
This table reports results for the nuclear norm regularized (NN-reg) estimator
with regularization parameter $\psi=\log(\log(7))7^{-1/2}/4$, the triad
pairwise differencing estimator (TPWDkit) obtained after $k$ iterations
starting at NN-reg and setting $c_{NT}=\widehat{\sigma}^{1}\log(7)/\sqrt{7}$,
and Bonhomme and Manresa (2015)’s approximating Algorithm 2 for the GFE
estimator with number of groups set to $g$ (GFEG=g). When available,
analytical standard errors based on large $N,T$ approximations and clustered
at the country level are shown in parentheses.

First, while the nuclear norm is quite different from the GFE estimates, the
8th and 9th iterations (at which the sequence converges) of TPWD estimates of
the lagged democracy and lagged outcome effects are in-between the NN-reg and
the GFE with several pre-specified numbers of groups on average. After three
iterations, TPWD converges to three estimated groups and delivers estimates
closer to the GFE approximation with three groups, which confirms the
intuition of Bonhomme and Manresa (2015). The income effect is statistically
significant and similar to the GFE, but the cumulative income effect is
significantly almost twice as large.

Estimated groups differ markedly from that reported by the GFE estimator with
3 groups (permanently low, permanently high, and transitioning democracy
groups), which turn into slightly different estimates. Figures 3 plot the time
effects at iteration 8 (left panel) and 9 (right panel). On the left panel, 79
countries exhibit a high and steadily increasing trend in democracy (Group 1),
while Algeria, Cameroon, Gabon, Guinea, Iran, Mauritania, Peru, Singapore, the
Syrian Arab Republic, and Tunisia exhibit a lower level of development that is
increasing for the period 1970-1990 and decreasing on 1990-1995, and Nigeria
is surprisingly volatile in terms of democracy (Group 3). On the right panel,
Peru switches from Group 2 to Group 1, which does not affect much the results.

Figure 2: Time effects $\widehat{\alpha}_{gt}$

[h]

(a) 8 iterations (b) 9 iterations

Figure 3: 2 Figures side by side

  * •

Notes: TPWD estimates with NN-reg as a preliminary estimator and data-driven
rule $c_{NT}=\widehat{\sigma}^{1}\log(7)/\sqrt{7}$ after 8 iterations (left
panel) and 9 iterations (right panel). The solid blue, dashed red, dash-dotted
yellow lines correspond to Group 1 (“steady democracy”), Group 2 (“volatile”),
and Group 3 (“highly volatile”) respectively.

##  6 Conclusion

Grouped fixed effects models are plagued with an underlying difficult
combinatorial classification problem, rendering estimation and inference
difficult. In this paper, I propose a novel strategy for constructive
identification of all the model parameters including the number of groups. The
method therefore simultaneously solve the estimation, classification, and
model selection problem. The corresponding three-step estimator has a
polynomial computational cost and is straightforward to implement: only smooth
convex optimization and elementary arithmetic operations are required. It is
based on thresholding a suitable pairwise differencing transformation of the
regression equation and a preliminary off-the-shelf consistent estimator of
the slope. Mild conditions are given under which the proposed estimator is
uniformly consistent for the latent grouping structure and asymptotically
normal as both dimensions diverge jointly. Importantly, the number of groups
is consistently estimated without any prior knowledge, and the time dimension
can grow much more slowly than the cross-sectional dimension. The paper leaves
a few questions unanswered. For instance: could the approach be fruitful to
build a test for the grouping assumption? Could similar differencing ideas be
applied to nonlinear structural models and/or potential outcome models? I
leave these questions for further research.

Appendix

##  Appendix A Sufficient conditions for consistency of the nuclear norm
regularized estimator

Define
$\gamma^{0}\equiv(\mathbf{1}\\{g_{i}^{0}=g\\})_{i=1,\ldots,N;g=1,\ldots,G^{0}}\in\left\\{0,1\right\\}^{N\times
G^{0}}$,
$\alpha^{0}\equiv(\alpha_{gt}^{0})_{t=1,\ldots,T;g=1,\ldots,G^{0}}\in\mathcal{A}^{T\times
G^{0}}$, $x_{k}\equiv\textrm{vec}(X_{k})$, and $x\equiv(x_{1},\ldots,x_{k})$.

###### Assumption 5

  * (a)

As $N$ and $T$ tend to infinity: $\psi\to 0$ such that
$\sqrt{\min(N,T)}\psi\to\infty$.

  * (b)

Let $\mathbb{C}\equiv\left\\{A\in\mathbb{R}^{N\times
T}:\big{\|}M_{\gamma^{0}}AM_{\alpha^{0}}\big{\|}_{1}\leq
3\big{\|}A-M_{\gamma^{0}}AM_{\alpha^{0}}\big{\|}_{1}\right\\}$, where
$M_{B}\equiv I-B(B^{\prime}B)^{\dagger}B^{\prime}$, $I$ is the identity matrix
of appropriate dimensions, and † refers to the Moore-Penrose generalized
inverse. There exists $\mu>0$, independent from $N$ and $T$, such that for any
$a\in\mathbb{R}^{NT}$ with ${\rm mat}(a)\in\mathbb{C}$,
$a^{\prime}M_{x}a\geq\mu a^{\prime}a$ holds for $N,T$ sufficiently large.

  * (c)

$\big{\|}(v_{it})_{i=1,\ldots,N;t=1,\ldots,T}\big{\|}_{\infty}=O_{p}\left(\sqrt{\max(N,T)}\right)$.

  * (d)

As $N$ and $T$ tend to infinity:
$\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}x_{it}x_{it}^{\prime}\overset{p}{\to}\Sigma>0$
and $\frac{1}{\sqrt{NT}}\sum_{i=1}^{N}\sum_{t=1}^{T}v_{it}x_{it}=O_{p}(1)$.

Assumption 5(b) is a restricted eigenvalue condition common in high-
dimensional modelling (e.g., Bickel, Ritov, and Tsybakov, 2009). Sufficient
conditions for Assumption 5(c) are given in Supplementary Appendix S.2 of Moon
and Weidner (2017).

######  Proposition A.1 (Moon and Weidner (2019))

Let Assumption 5 hold. Then, as $N$ and $T$ diverge jointly to infinity,
$\big{\|}\widehat{\beta}^{1}(\psi)-\beta^{0}\big{\|}=O_{p}(\psi)$.

Proof of Proposition A.1. This follows by Moon and Weidner (2019)’s Theorem 2
after noticing the interactive fixed effects structure of model (3.1) as in
Footnote 2.

##  Appendix B Extensions

The three-step procedure underlying the TPWD estimator can be applied to
several extensions of model (1.1). I briefly discuss a few of them in this
section.

###  B.1 Linear models with heterogeneous slopes/asymptotically close groups

To allow for unit-specific effects and/or unit/time-specific slopes, one may
work with time-differenced or demeaned equations and/or use a computationally
simple first-step estimator well suited for heterogeneous/time-varying slope
(e.g., Chernozhukov, Hansen, Liao, and Zhu, 2019). Group separation conditions
need to be adjusted similarly to what is discussed in Bonhomme and Manresa
(2015) ’s Supplemental Material, but the main arguments convey.

The approach of this paper allows for some degree of asymptotic closeness in
group-specific effects. Suppose that for all
$(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}$ such that
$g\neq\widetilde{g}$:

$T^{-1}\sum_{t=1}^{T}(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})^{2}\geq\rho_{NT},$

for some deterministic sequence $\rho_{NT}\to 0$ as $\min(N,T)\to\infty$.
Then, only the rate conditions in Assumption 2 have to be adapted to recover
the true group memberships with probability approaching one. The asymptotic
behaviour of the oracle OLS estimator, however, depends on this weak factor
property.

###  B.2 Nonlinear multiplicative models for networks

Consider dyadic observations
$\\{(y_{ij},x_{ij}):(i,j)\in\\{1,\ldots,n\\}^{2},i\neq j\\}$ for $n$ agents
such that

$y_{ij}=\varphi(x_{ij};\beta_{0})u_{ij},\quad i\neq j,$ (B.1)

where $\varphi:\mathbb{R}^{K}\to\mathbb{R}^{+*}$ is a function known up to the
finite-dimensional parameter vector $\beta_{0}\in\mathbb{R}^{K}$, and
$u_{ij}\in\mathbb{R}^{+}$ is a latent disturbance. Suppose

$u_{ij}=\alpha_{i}\gamma_{j}\omega_{g_{i},h_{j}}\varepsilon_{ij},\quad i\neq
j,$ (B.2)

where $\alpha_{i}\in\mathbb{R}^{+*}$ and $\gamma_{j}\in\mathbb{R}^{+*}$ are
permanent sender (exporter) and receiver (importer) unobserved effects,
$g_{i}\in\\{1,\ldots,G_{0}\\}$ is an unobserved exporter-group membership
variable, $G_{0}$ is the number of groups of exporters (considered exogenous
and fixed), $h_{j}\in\\{1,\ldots,H_{0}\\}$ is an unobserved importer group
membership variable, $H_{0}$ is the number of groups of importers (considered
exogenous and fixed), $\omega_{g,h}\in\mathbb{R}^{+*}$ is a permanent
unobserved effect affecting group $g$ exporting to group $h$, and
$\varepsilon_{ij}\in\mathbb{R}^{+}$ is an idiosyncratic disturbance such that
$\Pr(\varepsilon_{ij}=0)<1$. Here, I treat
$\alpha_{i},\gamma_{j},\omega_{g,h},g_{i},h_{j}$ as fixed, that is I condition
on them.

Model (B.1) extends Jochmans (2017) by allowing for latent grouped
interactions on top of standard sender and receiver effects. Such grouped
effects might capture nonlinear latent measures of reciprocity between units
$i$ and $j$: trade shocks that are shared by unobserved groups of exporters or
groups of importers. Suppose interest lies in estimating $\beta_{0}$,
$\alpha\equiv(\alpha_{1},\ldots,\alpha_{n})$,
$\gamma\equiv(\gamma_{1},\ldots,\gamma_{n})$, $g\equiv(g_{1},\ldots,g_{n})$,
$h\equiv(h_{1},\ldots,h_{n})$, and
$\Omega\equiv(\omega_{g,h})_{(g,h)\in[G_{0}]\times[H_{0}]}$ under the
conditional mean restriction

$\mathbb{E}[\varepsilon_{ij}|x_{12},\ldots,x_{n(n-1)}]=1,\quad\forall i\neq
j.$ (B.3)

Because $g,h,\alpha,\gamma$ and $\Omega$ are unobserved, one needs a
normalization specified as:777The first equality is a standard choice in the
literature Jochmans (2017); Dzemski (2019). It is equivalent to
$\sum_{i=1}^{n}\nu_{i}=\sum_{j=1}^{n}\xi_{j}=0$ where
$\nu_{i}\equiv\log(\alpha_{i})$ and $\xi_{j}\equiv\log(\gamma_{j})$. The
second equality is needed because of the introduction of the group effects.
These normalization choices are arbitrary and one could alternatively assume
$\alpha_{1}=1$ and $\gamma_{1}=1$ without affecting our results.

$\prod_{i=1}^{n}\alpha_{i}=\prod_{j=1}^{n}\gamma_{j}=1.$ (B.4)

Section S.2.1 outlines identification arguments based on triad-differencing.

###  B.3 Rubin-Holland Potential Outcomes

Consider a binary treatment $D_{it}\in\\{0,1\\}$ and potential outcomes
$Y_{it}(0)$ and $Y_{it}(1)$ if individual $i$ is non-treated or treated at
period $t$ respectively. Suppose a researcher wants to learn about the
treatment effect $Y_{it}(1)-Y_{it}(0)$, but she only observes $Y_{it}(D_{it})$
which is the fundamental problem of causal inference. Under a version of
parallel trends holding within each unobserved group
$g_{i}\in\\{1\ldots,N\\}$, one strategy is to use untreated periods to
estimate the group memberships under possibly heterogeneous treatment effects
by applying directly step 2 to differenced outcomes. Because $D_{it}=0$ for
this subsample, heterogeneous treatment effects do not perturbate the
estimation of the latent group structure. With a sufficiently large number of
pre-treatment outcomes, one consistently estimate the groups and can build
counterfactual in post-treatment periods to identify a menu of average
treatment effect. The theoretical analysis of this procedure generalizing de
Chaisemartin and D’Haultfœuille (2020) and a large-scale application is work
in progress.

##  Appendix C Proofs of the results

###  C.1 Proof of Proposition 3.1

Let
$W^{0}\equiv(\mathbf{1}\\{g_{i}^{0}=g_{j}^{0}\\})_{i=1,\ldots,N;j=1,\ldots,N}$.
Equations (3.2) and (3.3) are immediate corollaries of Lemma C.1 below.

###### Lemma C.1

Let Assumptions 1-3 hold. Then, as $N$ and $T$ tend to infinity,

$\big{\|}\widehat{W}-W^{0}\big{\|}_{\max}=o_{p}(1).$ (C.1)

Proof of Lemma C.1. Let $\epsilon,K_{1}>0$. By Assumption 1, there exists
$K_{2}>0$ such that, letting
$\mathcal{E}_{1NT}\equiv\left\\{\big{\|}\widehat{\beta}^{1}-\beta^{0}\big{\|}>K_{2}r_{NT}\right\\}$,
$\Pr(\mathcal{E}_{1NT})<\epsilon$ for $N,T$ sufficiently large. Define
$Z_{1NT}(i,j)\equiv\widehat{W}_{ij}(1-W_{ij}^{0})$,
$Z_{2NT}(i,j)\equiv(1-\widehat{W}_{ij})W_{ij}^{0}$, and the probability events

$\displaystyle\mathcal{E}_{2NT}$
$\displaystyle\equiv\left\\{\min_{g\in\left\\{1,\ldots,G^{0}\right\\}}\sum_{i=1}^{N}\mathbf{1}\\{g_{i}^{0}=g\\}\geq
2\right\\}\cap\left\\{C\leq c_{NT}T^{\kappa}\leq
T^{\kappa}C\right\\}\cap\left\\{c_{NT}r_{NT}^{-1}\geq K_{1}\right\\}$
$\displaystyle\qquad\cap\left\\{r_{NT}\leq\min(\eta/K_{2},1)\right\\}\cap\left\\{c_{NT}<c/72\right\\},$

where $\eta$ is defined in Equation (C.5) and $c,C$ in Assumptions 3(c) and
2(b), and $\mathcal{E}_{NT}\equiv\mathcal{E}_{1NT}^{c}\cap\mathcal{E}_{2NT}$.
By the union bound, for $N,T$ sufficiently large,

$\displaystyle\Pr\left(\max_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\left|\widehat{W}_{ij}-W_{ij}^{0}\right|>0\right)$
$\displaystyle\quad\leq\Pr(\mathcal{E}_{NT}^{c})+\sum_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\Pr\left(\left|\widehat{W}_{ij}\neq
W_{ij}^{0}\right|,\mathcal{E}_{NT}\right)$
$\displaystyle\quad\leq\Pr(\mathcal{E}_{1NT})+\Pr(\mathcal{E}_{2NT}^{c})+\sum_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\Pr\left(\left|\widehat{W}_{ij}\neq
W_{ij}^{0}\right|,\mathcal{E}_{NT}\right)$ $\displaystyle\quad\leq
2\epsilon+\sum_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\Pr(Z_{1NT}(i,j)=1,\mathcal{E}_{NT})+\Pr(Z_{2NT}(i,j)=1,\mathcal{E}_{NT}),$
(C.2)

where the last inequality follows from
$\lim_{\min(N,T)\to\infty}\Pr\left(\mathcal{E}_{2NT}^{c}\right)=0$ by
Assumptions 2(b) and 3(e). Below, I prove that, for
$\ell\in\left\\{1,2\right\\}$, and as $N$ and $T$ tend to infinity,

$\displaystyle\max_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\Pr(Z_{\ell
NT}(i,j)=1,\mathcal{E}_{NT})=o(N^{2}T^{-\delta})\text{ for all }\delta>0.$
(C.3)

Equation (C.1) then follows by combining (C.2)-(C.3) with Assumption 2(a) and
because $\epsilon$ is unrestricted.

1\. First, I show (C.3) for $\ell=1$.888Actually, I show the stronger result
that the supremum is $o(T^{-\delta})$. Let
$(i,j)\in\left\\{1,\ldots,N\right\\}^{2}$ and $\delta>0$.

$Z_{1NT}(i,j)=\mathbf{1}\Big{\\{}\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|\leq
c_{NT}\Big{\\}}\mathbf{1}\\{g_{i}^{0}\neq g_{j}^{0}\\}.$

If $G^{0}=1$, then almost surely $g_{i}^{0}=g_{j}^{0}$ and $Z_{1NT}(i,j)=0$,
i.e., (C.3) holds. Else,

$\displaystyle\mathbf{1}\\{Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\\}$
$\displaystyle\quad=\mathbf{1}\\{\mathcal{E}_{NT}\\}\times$
$\displaystyle\quad\quad\sum_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\\{g_{i}^{0}=g\\}\mathbf{1}\\{g_{j}^{0}=\widetilde{g}\\}\mathbf{1}\Big{\\{}\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|\leq
c_{NT}\Big{\\}}.$

If $\mathbf{1}\\{\mathcal{E}_{NT}\\}\mathbf{1}\\{g_{i}^{0}\neq
g_{j}^{0}\\}=1$, there exists a pair
$\left(k^{*}(i,j,g_{i}^{0}),l^{*}(i,j,g_{j}^{0})\right)\in\mathscr{P}_{2}\left(\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}\right)$
such that $g_{k^{*}(i,j,g_{i}^{0})}^{0}=g_{i}^{0}$ and
$g_{l^{*}(i,j,g_{j}^{0})}^{0}=g_{j}^{0}$. It follows that

$\displaystyle\mathbf{1}\\{Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\\}$
$\displaystyle\quad\leq\mathbf{1}\\{\mathcal{E}_{NT}\\}\times$
$\displaystyle\quad\quad\sum_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\\{g_{i}^{0}=g\\}\mathbf{1}\\{g_{j}^{0}=\widetilde{g}\\}\mathbf{1}\Big{\\{}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{k^{*}(i,j,g_{i}^{0})t}\right|\leq
c_{NT}\Big{\\}}\times$
$\displaystyle\quad\quad\quad\mathbf{1}\Big{\\{}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{l^{*}(i,j,g_{j}^{0})t}\right|\leq
c_{NT}\Big{\\}}$
$\displaystyle\quad\leq\mathbf{1}\\{\mathcal{E}_{NT}\\}\times$
$\displaystyle\quad\quad\sum_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\\{g_{i}^{0}=g\\}\mathbf{1}\\{g_{j}^{0}=\widetilde{g}\\}\mathbf{1}\Big{\\{}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})(\widehat{v}_{k^{*}(i,j,g_{i}^{0})t}-\widehat{v}_{l^{*}(i,j,g_{j}^{0})t})\right|\leq
2c_{NT}\Big{\\}},$

where the first inequality uses the definition of the maximum and the second
inequality follows from the triangle inequality. Since there is at most one
pair $(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}$ such that
$g\neq\widetilde{g}$ and
$\mathbf{1}\\{g_{i}^{0}=g\\}\mathbf{1}\\{g_{j}^{0}=\widetilde{g}\\}=1$,
developing the product and using $\mathbf{1}\\{\left|a\right|\leq
b\\}\leq\mathbf{1}\\{a\leq b\\}$ for any
$(a,b)\in\mathbb{R}\times\mathbb{R}^{*}$ yields

$\displaystyle\mathbf{1}\\{Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\\}$
$\displaystyle\quad\leq\mathbf{1}\\{\mathcal{E}_{NT}\\}\times$
$\displaystyle\quad\quad\max_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\Big{\\{}\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}+\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)\left(v_{it}-v_{jt}+v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad\quad+\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{it}-x_{jt}+x_{k^{*}(i,j,g)t}-x_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad\quad+\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\left(v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad\quad+\frac{1}{T}\sum_{t=1}^{T}\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{it}-x_{jt}\right)\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{k^{*}(i,j,g)t}-x_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad\quad+\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{k^{*}(i,j,g)t}-x_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad\quad+\frac{1}{T}\sum_{t=1}^{T}(v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t})\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{it}-x_{jt}\right)\leq
2c_{NT}\Big{\\}}$
$\displaystyle\quad=\mathbf{1}\\{\mathcal{E}_{NT}\\}\times\max_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\Big{\\{}A_{T}(i,j,g,\widetilde{g})\leq
2c_{NT}\Big{\\}},$ (C.4)

where $A_{T}(i,j,g,\widetilde{g})$ is defined implicitely. Define

$\displaystyle B_{T}(i,j,g,\widetilde{g})$
$\displaystyle\equiv\left|A_{T}(i,j,g,\widetilde{g})-\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}\right.$
$\displaystyle\left.\quad-\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)\left(v_{it}-v_{jt}+v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)\right.$
$\displaystyle\left.\quad-\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\left(v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)\right|.$

Let $\overline{a}\equiv\sup_{a\in\mathcal{A}}\left|a\right|<\infty$ by
Assumption 3(a). It is easy to show using the Cauchy-Schwarz inequality that

$\displaystyle B_{T}(i,j,g,\widetilde{g})$
$\displaystyle\leq\big{\|}\widehat{\beta}^{1}-\beta^{0}\big{\|}\left\\{\frac{2\overline{a}}{T}\sum_{t=1}^{T}\left(\big{\|}x_{it}\big{\|}+\big{\|}x_{jt}\big{\|}+\big{\|}x_{k^{*}(i,j,g)t}\big{\|}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}\right)\right.$
$\displaystyle\left.\quad+\frac{4\big{\|}\widehat{\beta}^{1}-\beta^{0}\big{\|}}{T}\sum_{t=1}^{T}\left(\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}+\big{\|}x_{k^{*}(i,j,g)t}\big{\|}^{2}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}^{2}\right)\right.$
$\displaystyle\left.\quad+\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}}+\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{jt}^{2}}\right)\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{k^{*}(i,j,g)t}\big{\|}^{2}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}^{2}}\right.$
$\displaystyle\left.\quad+\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{k^{*}(i,j,g)t}^{2}}+\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{l^{*}(i,j,\widetilde{g})t}^{2}}\right)\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}}\right\\}.$

By Assumption 3(b), there exists $M^{*}>0$ such that
$\mathbb{E}[v_{it}^{2}]\leq M^{*}$ for all $i,t$. Let
$\widetilde{M}>\max(M,\max(M^{*},1))$, where $M$ is defined in Assumption
3(f), and $\eta>0$ such that

$\eta\leq\min\left(1,\frac{c}{24\left(2\overline{a}4\sqrt{\widetilde{M}}+8\widetilde{M}+4\sqrt{2}\widetilde{M}\right)}\right),$
(C.5)

where $c$ is defined in Assumption 3(c). By definition of $\mathcal{E}_{NT}$,
$\big{\|}\widehat{\beta}^{1}-\beta^{0}\big{\|}\leq\eta$ on $\mathcal{E}_{NT}$.
Using the Cauchy-Schwarz inequality again and $\eta\leq 1$ yields

$\displaystyle\mathbf{1}\\{\mathcal{E}_{NT}\\}B_{T}(i,j,g,\widetilde{g})$
$\displaystyle\leq\eta\left\\{2\overline{a}\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}+\big{\|}x_{k^{*}(i,j,g)t}\big{\|}^{2}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}^{2}}\right.$
$\displaystyle\left.\quad+\frac{4}{T}\sum_{t=1}^{T}\left(\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}+\big{\|}x_{k^{*}(i,j,g)t}\big{\|}^{2}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}^{2}\right)\right.$
$\displaystyle\left.\quad+\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}}+\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{jt}^{2}}\right)\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{k^{*}(i,j,g)t}\big{\|}^{2}+\big{\|}x_{l^{*}(i,j,\widetilde{g})t}\big{\|}^{2}}\right.$
$\displaystyle\left.\quad+\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{k^{*}(i,j,g)t}^{2}}+\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{l^{*}(i,j,\widetilde{g})t}^{2}}\right)\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}}\right\\}$
$\displaystyle\equiv C_{T}(i,j,g,\widetilde{g}).$

Plugging this upper bound into (C.4), I obtain

$\displaystyle\mathbf{1}\\{Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\\}$
$\displaystyle\quad\leq\max_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\mathbf{1}\Big{\\{}\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}$
$\displaystyle\quad+\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)\left(v_{it}-v_{jt}+v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)$
$\displaystyle\quad+\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\left(v_{k^{*}(i,j,g)t}-v_{l^{*}(i,j,\widetilde{g})t}\right)\leq
2c_{NT}+C_{T}(i,j,g,\widetilde{g})\Big{\\}}.$

By the Cauchy-Schwarz inequality again, and because $\widetilde{M}\geq 1$,
note the implication

$\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}\leq\widetilde{M}\implies\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}\leq\sqrt{\widetilde{M}}\leq\widetilde{M}.$

Using this result, the union bound, and some probability algebra, it follows
that

$\displaystyle\Pr\left(Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\right)$
$\displaystyle\quad\leq\sum_{\begin{subarray}{c}(g,\widetilde{g})\in\left\\{1,\ldots,G^{0}\right\\}^{2}\\\
g\neq\widetilde{g}\end{subarray}}\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\leq-\frac{c}{12}+2c_{NT}+\eta\left(2\overline{a}4\sqrt{\widetilde{M}}+8\widetilde{M}+4\sqrt{2}\widetilde{M}\right)\right)$
$\displaystyle\quad\quad+4G^{0}\left(G^{0}-1\right)\left[\sup_{g\neq\widetilde{g}}\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}\leq\frac{c}{2}\right)+\sup_{i\in\left\\{1,\ldots,N\right\\}}\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}\geq\widetilde{M}\right)\right.$
$\displaystyle\left.\quad\quad+\sup_{i\in\left\\{1,\ldots,N\right\\},g\neq\widetilde{g}}Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\right|\geq\frac{c}{12}\right)+\sup_{i\in\left\\{1,\ldots,N\right\\}}\Pr\left(\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}\geq\widetilde{M}\right)\right.$
$\displaystyle\left.\qquad+\sup_{(i,j,k)\in\mathscr{P}_{3}(\left\\{1,\ldots,N\right\\})}\Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})v_{kt}\right|\geq\frac{c}{12}\right)\right].$
(C.6)

First, I bound the terms with a supremum. By Assumption 3(c), it holds that
$\lim_{\min(N,T)\to\infty}\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}[(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})^{2}]=c_{g,\widetilde{g}}>c$.
So for $N,T$ large enough, I have

$\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\left[\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}\right]\geq\frac{2c}{3}.$

Applying Lemma B.5 in Bonhomme and Manresa (2015) to
$z_{t}=(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})^{2}-\mathbb{E}[(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})^{2}]$,
which satisfies appropriate mixing and tail conditions by Assumption 3(b) and
(d), and taking $z=c/6$ yields, as $N,T$ tends to infinity,

$\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)^{2}\leq\frac{c}{2}\right)=o(T^{-\delta}),$
(C.7)

uniformly across $g$ and $\widetilde{g}$. Similarly, applying Lemma B.5 to
$z_{t}=v_{it}^{2}-\mathbb{E}[v_{it}]^{2}$ and taking $z=\widetilde{M}-M^{*}$
yields

$\Pr\left(\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}\geq\widetilde{M}\right)=o(T^{-\delta}),$
(C.8)

uniformly across units $i$. Note that $\left\\{v_{it}^{2}\right\\}_{t}$ is
strongly mixing as $\left\\{v_{it}\right\\}$ is strongly mixing by Assumption
3(d). By Assumption 3(d), the process
$\left\\{(\alpha_{\widetilde{g}t}^{0}-\alpha_{gt}^{0})v_{it}\right\\}_{t}$ has
zero mean, and is strongly mixing with faster-than-polynomial decay rate.
Moreover, for all $i,t$ and $m>0$,

$\Pr\left(\left|\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\right|>m\right)\leq\Pr\left(\left|v_{it}\right|>\frac{m}{2\overline{a}}\right),$

so $\left\\{(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})v_{it}\right\\}_{t}$
also satisfies the tail condition of Assumption 3(b), albeit with a different
constant $b^{\prime}>0$ instead of $b>0$. Lastly, applying Lemma B.5 from
Bonhomme and Manresa (2015) again with
$z_{t}=(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0})v_{it}$ and taking
$z=c/12$ yields

$\Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\right|\geq\frac{c}{12}\right)=o(T^{-\delta})$
(C.9)

uniformly across $i$, $g$, and $\widetilde{g}$. An analogous reasoning yields

$\sup_{(i,j,k)\in\mathscr{P}_{3}(\left\\{1,\ldots,N\right\\})}\Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})v_{kt}\right|\geq\frac{c}{12}\right)=o(T^{-\delta}).$
(C.10)

Finally, because $c_{NT}\leq c/72$ on $\mathcal{E}_{NT}$, a similar reasoning
yields

$\displaystyle\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\leq-\frac{c}{12}+2c_{NT}+\eta\left(2\overline{a}4\sqrt{\widetilde{M}}+8\widetilde{M}+4\sqrt{2}\widetilde{M}\right)\right)$
$\displaystyle\quad\leq\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\left(\alpha_{gt}^{0}-\alpha_{\widetilde{g}t}^{0}\right)v_{it}\leq-\frac{c}{72}\right)$
$\displaystyle\quad=o(T^{-\delta}),$ (C.11)

uniformly across $g,\widetilde{g}$, where I have used the value of $\eta$
given in (C.5). Combining (C.6)-(C.11) and using Assumption 3(f) yields

$\sup_{(i,j)\in\left\\{1,\ldots,N\right\\}^{2}}\Pr\left(Z_{1NT}(i,j)=1,\mathcal{E}_{NT}\right)=G^{0}(1-G^{0})\times
o_{p}(T^{-\delta})=o_{p}(T^{-\delta}),$

i.e., (C.3) for $\ell=1$ holds.

2\. Second, I show (C.3) for $\ell=2$. I now have

$\displaystyle\mathbf{1}\\{Z_{2NT}(i,j)=1,\mathcal{E}_{NT}\\}$
$\displaystyle\quad=\mathbf{1}\\{\mathcal{E}_{NT}\\}\mathbf{1}\Big{\\{}\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(\widehat{v}_{it}-\widehat{v}_{jt})\widehat{v}_{kt}\right|>c_{NT}\Big{\\}}\mathbf{1}\\{g_{i}^{0}=g_{j}^{0}\\}$
$\displaystyle\quad\leq\mathbf{1}\\{\mathcal{E}_{NT}\\}\mathbf{1}\Big{\\{}\max_{k\in\left\\{1,\ldots,N\right\\}\backslash\left\\{i,j\right\\}}\left|\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})v_{kt}+\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\alpha_{g_{k}^{0}t}^{0}\right.$
$\displaystyle\left.\quad\quad+\frac{1}{T}\sum_{t=1}^{T}\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{it}-x_{jt}\right)\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}x_{kt}\right.$
$\displaystyle\left.\quad\quad+\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}x_{kt}+\frac{1}{T}\sum_{t=1}^{T}\alpha_{g_{k}^{0}t}^{0}\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}(x_{it}-x_{jt})\right.$
$\displaystyle\left.\quad\quad+\frac{1}{T}\sum_{t=1}^{T}v_{kt}\left(\beta^{0}-\widehat{\beta}^{1}\right)^{\prime}\left(x_{it}-x_{jt}\right)\right|>c_{NT}\Big{\\}}.$

By the union bound, the triangle inequality, and the Cauchy-Schwarz
inequality,

$\displaystyle\Pr\left(Z_{2NT}(i,j)=1,\mathcal{E}_{NT}\right)$
$\displaystyle\quad\leq(N-2)\sup_{(i,j,k)\in\left\\{1,\ldots,N\right\\}^{3}}\left\\{\Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})v_{kt}\right|>\frac{c_{NT}}{10},\mathcal{E}_{NT}\right)\right.$
$\displaystyle\left.\qquad+\Pr\left(\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{it}\big{\|}^{2}+\big{\|}x_{jt}\big{\|}^{2}+\big{\|}x_{kt}\big{\|}^{2}>\frac{c_{NT}}{10\times
4K_{2}^{2}r_{NT}^{2}},\mathcal{E}_{NT}\right)\right.$
$\displaystyle\qquad\left.+4\Pr\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}}\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{kt}\big{\|}^{2}}>\frac{c_{NT}}{10K_{2}r_{NT}},\mathcal{E}_{NT}\right)\right.$
$\displaystyle\qquad\left.+2\Pr\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}}>\frac{c_{NT}}{10K_{2}r_{NT}\times\overline{a}}\right)\right.$
$\displaystyle\left.\qquad+2\Pr\left(\sqrt{\frac{1}{T}\sum_{t=1}^{T}v_{it}^{2}}\sqrt{\frac{1}{T}\sum_{t=1}^{T}\big{\|}x_{kt}\big{\|}^{2}}>\frac{c_{NT}}{10K_{2}r_{NT}\times\overline{a}},\mathcal{E}_{NT}\right)\right\\}.$

Under the strong mixing and tail conditions given by Assumptions 3(b) and
3(d), and because $K_{1}\leq c_{NT}/r_{NT}\leq c_{NT}/r_{NT}^{2}$ on
$\mathcal{E}_{NT}$, for $K_{1}$ sufficiently large, all noninitial
probabilities in the above expression can be shown to be $o(T^{-\delta})$ for
all $\delta>0$, uniformly accross $(i,j,k)$, by similar arguments as in Step
1. For the first probability, a close inspection of the proof of Lemma B.5 in
Bonhomme and Manresa (2015) reveals that, by taking
$z_{t}=(v_{it}-v_{jt})v_{kt}$ and $z=c_{NT}/6$, and because $c_{NT}\geq
CT^{-\kappa}$ on $\mathcal{E}_{NT}$,

$\displaystyle\Pr\left(\left|\frac{1}{T}\sum_{t=1}^{T}(v_{it}-v_{jt})v_{kt}\right|\geq\frac{c_{NT}}{10}\right)$
$\displaystyle\leq
4\left(1+\frac{T^{1/2-2\kappa}}{C_{1}}\right)^{-(1/2)T^{1/2}}$
$\displaystyle\qquad+C_{2}T^{\kappa}\exp\left(-C_{3}\left(T^{(1/2-\kappa)/C_{4}}\right)\right),$
(C.12)

where $C_{1},C_{2},C_{3}$, and $C_{4}$ are positive constants that do not
depend on $i,j,k$. Since $\kappa<1/2$, the upper bound is $o_{p}(T^{-\delta})$
for all $\delta>0$. This shows (C.3) for $\ell=2$. $\square$

###  C.2 Proof of Corollary 3.2

Let $\widetilde{\beta}$ and
$(\widetilde{\alpha}_{11},\ldots,\widetilde{\alpha}_{G^{0}T})^{\prime}$ denote
the infeasible oracle estimators computed from a pooled OLS regression of
$y_{it}$ on $x_{it}$ and the interactions of group and time indicators
$\mathbf{1}\\{g_{i}^{0}=1\\},\ldots,\mathbf{1}\\{g_{i}^{0}=G^{0}\\}$,
$\mathbf{1}\\{t=1\\},\ldots,\mathbf{1}\\{t=T\\}$. By the same reasoning as in
section S.A.1. in Bonhomme and Manresa (2015)’s Supplemental Material, it
holds

$\sqrt{NT}(\widetilde{\beta}-\beta^{0})\stackrel{{\scriptstyle
d}}{{\to}}\mathcal{N}\left(0,\Sigma_{\beta}^{-1}\Omega_{\beta}\Sigma_{\beta}^{-1}\right),$
(C.13)

and, for all
$(g,t)\in\left\\{1,\ldots,G^{0}\right\\}\times\left\\{1,\ldots,T\right\\}$,

$\sqrt{N}(\widetilde{\alpha}_{gt}-\alpha_{gt}^{0})\stackrel{{\scriptstyle
d}}{{\to}}\mathcal{N}\left(0,\frac{\omega_{gt}}{\pi_{g}^{2}}\right).$ (C.14)

Without loss of generality, I assume that the labelling of predicted groups
matches the true group labeling. By Proposition 3.1, for all
$(g,t)\in\left\\{1,\ldots,G^{0}\right\\}\times\left\\{1,\ldots,T\right\\}$,

$\displaystyle\Pr\left(\left\\{\widehat{\alpha}_{gt}\neq\widetilde{\alpha}_{gt}\right\\}\cup\left\\{\widehat{\beta}\neq\widetilde{\beta}\right\\}\right)$
$\displaystyle\leq\Pr\left(\widehat{G}\neq
G^{0}\right)+\Pr\left(\max_{i\in\left\\{1,\ldots,N\right\\}}\;\left|\widehat{g}_{i}-g_{i}^{0}\right|>0\right)$
$\displaystyle=o(1)+o(1)$ $\displaystyle=o(1).$

Then, Eq. (3.5) follows from

$\displaystyle\left|\Pr\left(\sqrt{N}(\widehat{\alpha}_{gt}-\alpha_{gt}^{0})\leq
a\right)-\Pr\left(\sqrt{N}(\widetilde{\alpha}_{gt}-\alpha_{gt}^{0})\leq
a\right)\right|$
$\displaystyle\quad\leq\left|\Pr\left(\sqrt{N}(\widehat{\alpha}_{gt}-\alpha_{gt}^{0})\leq
a,\sqrt{N}(\widetilde{\alpha}_{gt}-\alpha_{gt}^{0})>a\right)\right|$
$\displaystyle\qquad+\left|\Pr\left(\sqrt{N}(\widehat{\alpha}_{gt}-\alpha_{gt}^{0})>a,\sqrt{N}(\widetilde{\alpha}_{gt}-\alpha_{gt}^{0})\leq
a\right)\right|$
$\displaystyle\quad\leq\quad\Pr\left(\widehat{\alpha}_{gt}\neq\widetilde{\alpha}_{gt}\right)+\Pr\left(\widehat{\alpha}_{gt}\neq\widetilde{\alpha}_{gt}\right)=o(1).$

for any $a>0$. Eq. (3.4) follows from a similar argument.

## References

  * (1)
  * Acemoglu, Johnson, Robinson, and Yared (2008) Acemoglu, D., S. Johnson, J. A. Robinson, and P. Yared (2008): “Income and Democracy,” _American Economic Review_ , 98(3), 808–42.
  * Ando and Bai (2022) Ando, T., and J. Bai (2022): “Large-scale generalized linear longitudinal data models with grouped patterns of unobserved heterogeneity,” Discussion paper.
  * Arellano and Bond (1991) Arellano, M., and S. Bond (1991): “Some Tests of Specification for Panel Data: Monte Carlo Evidence and an Application to Employment Equations,” _The Review of Economic Studies_ , 58(2), 277–297.
  * Armstrong, Weidner, and Zeleneev (2022) Armstrong, T. B., M. Weidner, and A. Zeleneev (2022): “Robust Estimation and Inference in Panels with Interactive Fixed Effects,” .
  * Armstrong, Weidner, and Zeleneev (2023) (2023): “Robust Estimation and Inference in Panels with Interactive Fixed Effects,” .
  * Auerbach (2022) Auerbach, E. (2022): “Identification and Estimation of a Partially Linear Regression Model Using Network Data,” _Econometrica_ , 90(1), 347–365.
  * Bai (2003) Bai, J. (2003): “Inferential Theory for Factor Models of Large Dimensions,” _Econometrica_ , 71(1), 135–171.
  * Bai (2009) (2009): “Panel Data Models With Interactive Fixed Effects,” _Econometrica_ , 77(4), 1229–1279.
  * Barro (1999) Barro, R. J. (1999): “Determinants of Democracy,” _Journal of Political Economy_ , 107(S6), S158–S183.
  * Beyhum and Gautier (2019) Beyhum, J., and E. Gautier (2019): “Square-root nuclear norm penalized estimator for panel data models with approximately low-rank unobserved heterogeneity,” .
  * Beyhum and Gautier (2023) (2023): “Factor and Factor Loading Augmented Estimators for Panel Regression With Possibly Nonstrong Factors,” _Journal of Business & Economic Statistics_, 41(1), 270–281.
  * Bickel, Ritov, and Tsybakov (2009) Bickel, P. J., Y. Ritov, and A. B. Tsybakov (2009): “Simultaneous analysis of Lasso and Dantzig selector,” _The Annals of Statistics_ , 37(4), 1705 – 1732.
  * Bonhomme, Lamadon, and Manresa (2019) Bonhomme, S., T. Lamadon, and E. Manresa (2019): “A Distributional Framework for Matched Employer Employee Data,” _Econometrica_ , 87(3), 699–739.
  * Bonhomme, Lamadon, and Manresa (2022) (2022): “Discretizing Unobserved Heterogeneity,” _Econometrica_ , 90(2), 625–643.
  * Bonhomme and Manresa (2015) Bonhomme, S., and E. Manresa (2015): “Grouped Patterns of Heterogeneity in Panel Data,” _Econometrica_ , 83(3), 1147–1184.
  * Brownlees, Guðmundsson, and Lugosi (2022) Brownlees, C., G. S. Guðmundsson, and G. Lugosi (2022): “Community Detection in Partial Correlation Network Models,” _Journal of Business & Economic Statistics_, 40(1), 216–226.
  * Charbonneau (2017) Charbonneau, K. B. (2017): “Multiple fixed effects in binary response panel data models,” _The Econometrics Journal_ , 20(3), S1–S13.
  * Chen (2019) Chen, J. (2019): “Estimating latent group structure in time-varying coefficient panel data models,” _The Econometrics Journal_ , 22(3), 223–240.
  * Cheng, Schorfheide, and Shao (2019) Cheng, X., F. Schorfheide, and P. Shao (2019): “Clustering for multi-dimensional heterogeneity,” Discussion paper, Working paper.
  * Chernozhukov, Hansen, Liao, and Zhu (2019) Chernozhukov, V., C. Hansen, Y. Liao, and Y. Zhu (2019): “Inference for Heterogeneous Effects using Low-Rank Estimation of Factor Slopes,” .
  * Chetverikov and Manresa (2021) Chetverikov, D., and E. Manresa (2021): “Spectral and Post-Spectral Estimators for Grouped Panel Data Models,” Discussion paper.
  * de Chaisemartin and D’Haultfœuille (2020) de Chaisemartin, C., and X. D’Haultfœuille (2020): “Two-Way Fixed Effects Estimators with Heterogeneous Treatment Effects,” _American Economic Review_ , 110(9), 2964–96.
  * Dzemski (2019) Dzemski, A. (2019): “An Empirical Model of Dyadic Link Formation in a Network with Unobserved Heterogeneity,” _The Review of Economics and Statistics_ , 101(5), 763–776.
  * Dzemski and Okui (2018) Dzemski, A., and R. Okui (2018): “Confidence set for group membership,” .
  * Freeman and Weidner (2023) Freeman, H., and M. Weidner (2023): “Linear panel regressions with two-way unobserved heterogeneity,” _Journal of Econometrics_ , 237(1), 105498.
  * Graf and Luschgy (2002) Graf, S., and H. Luschgy (2002): “Rates of convergence for the empirical quantization error,” _The Annals of Probability_ , 30(2), 874 – 897.
  * Graham (2017) Graham, B. S. (2017): “An Econometric Model of Network Formation With Degree Heterogeneity,” _Econometrica_ , 85(4), 1033–1063.
  * Gu and Volgushev (2019) Gu, J., and S. Volgushev (2019): “Panel Data Quantile Regression With Grouped Fixed Effects,” _Journal of Econometrics_ , 213(1), 68 – 91, Annals: In Honor of Roger Koenker.
  * Hahn and Moon (2010) Hahn, J., and H. R. Moon (2010): “Panel Data Models With Finite Number of Multiple Equilibria,” _Econometric Theory_ , 26(3), 863–881.
  * Higgins (2022) Higgins, A. (2022): “Fixed T Estimation of Linear Panel Data Models with Interactive Fixed Effects,” .
  * Hocking, Joulin, Bach, and Vert (2011) Hocking, T., A. Joulin, F. Bach, and J. Vert (2011): “Clusterpath: An algorithm for clustering using convex fusion penalties,” in _Proceedings of the 28th International Conference on Machine Learning, ICML 2011_ , Proceedings of the 28th International Conference on Machine Learning, ICML 2011, pp. 745–752. 28th International Conference on Machine Learning, ICML 2011 ; Conference date: 28-06-2011 Through 02-07-2011.
  * Honoré and Powell (1994) Honoré, B. E., and J. L. Powell (1994): “Pairwise difference estimators of censored and truncated regression models,” _Journal of Econometrics_ , 64(1), 241–278.
  * Janys and Siflinger (2023) Janys, L., and B. Siflinger (2023): “Mental Health and Abortions among Young Women: Time-Varying Unobserved Heterogeneity, Health Behaviors, and Risky Decisions,” _Journal of Econometrics, Forthcoming_.
  * Jochmans (2017) Jochmans, K. (2017): “Two-Way Models for Gravity,” _The Review of Economics and Statistics_ , 99(3), 478–485.
  * Krasnokutskaya, Song, and Tang (2022) Krasnokutskaya, E., K. Song, and X. Tang (2022): “Estimating unobserved individual heterogeneity using pairwise comparisons,” _Journal of Econometrics_ , 226(2), 477–497.
  * Lewis, Melcangi, Pilossoph, and Toner-Rodgers (2022) Lewis, D. J., D. Melcangi, L. Pilossoph, and A. Toner-Rodgers (2022): “Approximating grouped fixed effects estimation via fuzzy clustering regression,” _Journal of Applied Econometrics, Forthcoming_.
  * Lipset (1959) Lipset, S. M. (1959): “Some Social Requisites of Democracy: Economic Development and Political Legitimacy,” _The American Political Science Review_ , 53(1), 69–105.
  * Lovász (2012) Lovász, L. (2012): _Large Networks and Graph Limits._ , vol. 60 of _Colloquium Publications_. American Mathematical Society.
  * Lumsdaine, Okui, and Wang (2023) Lumsdaine, R. L., R. Okui, and W. Wang (2023): “Estimation of panel group structure models with structural breaks in group memberships and coefficients,” _Journal of Econometrics_ , 233(1), 45–65.
  * Mammen, Wilke, and Zapp (2022) Mammen, E., R. Wilke, and K. M. Zapp (2022): “Estimation of group structures in panel models with individual fixed effects,” ZEW Discussion Papers 22-023, ZEW - Leibniz Centre for European Economic Research.
  * Mehrabani (2023) Mehrabani, A. (2023): “Estimation and identification of latent group structures in panel data,” _Journal of Econometrics_ , 235(2), 1464–1482.
  * Moon and Weidner (2017) Moon, H., and M. Weidner (2017): “Dynamic Linear Panel Regression Models With Interactive Fixed Effects,” _Econometric Theory_ , 33(1), 158–195.
  * Moon and Weidner (2015) Moon, H. R., and M. Weidner (2015): “Linear Regression for Panel With Unknown Number of Factors as Interactive Fixed Effects,” _Econometrica_ , 83(4), 1543–1579.
  * Moon and Weidner (2019) (2019): “Nuclear Norm Regularized Estimation of Panel Regression Models,” .
  * Mugnier (2022a) Mugnier, M. (2022a): “Supplemental Material for A Simple and Computationally Trivial Estimator for Grouped Fixed Effects Models,” https://martinmugnier.github.io/research.
  * Mugnier (2022b) (2022b): “Unobserved Clusters of Time-Varying Heterogeneity in Nonlinear Panel Data Models,” Discussion paper.
  * Ng, Jordan, and Weiss (2002) Ng, A., M. Jordan, and Y. Weiss (2002): “On Spectral Clustering: Analysis and an algorithm,” in _Advances in Neural Information Processing Systems_ , ed. by T. Dietterich, S. Becker, and Z. Ghahramani, vol. 14. MIT Press.
  * Pesaran (2006) Pesaran, M. H. (2006): “Estimation and Inference in Large Heterogeneous Panels with a Multifactor Error Structure,” _Econometrica_ , 74(4), 967–1012.
  * Pollard (1981) Pollard, D. (1981): “Strong Consistency of $K$-Means Clustering,” _The Annals of Statistics_ , 9(1), 135 – 140.
  * Pollard (1982) (1982): “A Central Limit Theorem for $k$-Means Clustering,” _The Annals of Probability_ , 10(4), 919 – 926.
  * Su, Shi, and Phillips (2016) Su, L., Z. Shi, and P. C. B. Phillips (2016): “Identifying Latent Structures in Panel Data,” _Econometrica_ , 84(6), 2215–2264.
  * Tan and Witten (2015) Tan, K. M., and D. Witten (2015): “Statistical properties of convex clustering,” _Electronic Journal of Statistics_ , 9(2), 2324 – 2347.
  * Vogt and Linton (2016) Vogt, M., and O. Linton (2016): “Classification of Non-Parametric Regression Functions in Longitudinal Data Models,” _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 79(1), 5–27.
  * von Luxburg (2007) von Luxburg, U. (2007): “A Tutorial on Spectral Clustering,” .
  * Wang and Su (2021) Wang, W., and L. Su (2021): “Identifying Latent Group Structures in Nonlinear Panels,” _Journal of Econometrics_ , 220(2), 272–295, Annals Issue: Celebrating 40 Years of Panel Data Analysis: Past, Present and Future.
  * Wooldridge (2010) Wooldridge, J. (2010): _Econometric Analysis of Cross Section and Panel Data, second edition_ , Econometric Analysis of Cross Section and Panel Data. MIT Press.
  * Yu, Gu, and Volgushev (2022) Yu, L., J. Gu, and S. Volgushev (2022): “Group structure estimation for panel data – a general approach,” .
  * Zeleneev (2020) Zeleneev, A. (2020): “Identification and Estimation of Network Models with Nonparametric Unobserved Heterogeneity,” Discussion paper.
  * Zhang, Levina, and Zhu (2017) Zhang, Y., E. Levina, and J. Zhu (2017): “Estimating network edge probabilities by neighbourhood smoothing,” _Biometrika_ , 104(4), 771–783.