Spectral and post-spectral estimators for grouped panel data models

Spectral and Post-Spectral Estimators for Grouped Panel Data Models

\address[D. Chetverikov]{ Department of Economics, UCLA, Bunche Hall, 8283, 315 Portola Plaza, Los Angeles, CA 90095, USA. } \email{[email removed]}

\address[E. Manresa]{ Department of Economics, New York University, 19 West 4th Street, New York, NY 10003, USA. } \email{[email removed]}

Introduction

Consider a grouped panel data model

where $i$ denotes cross-sectional units, $t$ denotes time periods, $y_{it}\in\mathbb R$ is an observable dependent variable, $x_{it} = (x_{it1},\dots,x_{itd})'\in\mathbb R^d$ is a corresponding vector of observable covariates, $g_i\in\{1,\dots,G\}$ is an unobservable group-membership variable, $v_{i t}\in\mathbb R$ is an unobservable zero-mean noise random variable, $\beta = (\beta_1,\dots,\beta_d)'\in\mathbb R^d$ is a vector of parameters of interest, and $(\alpha_{1 t},\dots,\alpha_{G t})' \in\mathbb R^G$ is a vector of unobservable group-specific time effects. Here, we assume that the noise and covariates are uncorrelated,

where $0_d = (0,\dots,0)'\in\mathbb R^d$, but group-specific time effects and group-membership variables can be arbitrarily correlated with covariates. Also, throughout the paper, we assume that the number of groups $G$ is known (consistently estimated).

The model (ref) was originally introduced by Bonhomme and Manresa in BM15, who also developed a so-called grouped fixed-effect estimator of the vector of parameters $\beta$ in this model. This estimator has attractive theoretical properties but is computationally difficult. It is therefore of interest to see if there exist alternative estimators that would be easier to compute. In this paper, we answer this question affirmatively, under certain additional assumptions to be specified in the next paragraph, and propose an estimator of $\beta$, which we call the {\em post-spectral estimator}, that also has nice theoretical properties but, in contrast to the grouped fixed-effect estimator, is computationally simple.

Like in the previous papers on grouped panel data models, we consider large $(N,T)$-asymptotics, i.e. we assume that $T\to\infty$, potentially very slowly, as $N\to\infty$, since otherwise $\beta$ is in general not identified. In contrast to the previous papers, however, we impose a special structure on the data-generating process for the covariates $x_{it}$. In particular, we assume that for some $M\geq 1$,

where $(\alpha_{1t}^m,\dots,\alpha_{Gt}^m)'\in\mathbb R^G$ for $m=1,\dots,M$ are group-specific time effects, $\rho_{im}\in \mathbb R^d$ for $m=1,\dots,M$ are individual-specific vectors of coefficients, and $z_{it}$ is a zero-mean component of $x_{it}$ that is independent of group-specific time effects, group-membership variable $g_{i}$, and vectors of coefficients $\rho_{i1},\dots\rho_{iM}$. Here, no quantity on the right-hand side of (ref) is observed, except for the number of time effects $M$, which we assume to be known (consistently estimated). Also, without loss of generality, we assume that $\alpha_{\gamma t}^1 = \alpha_{\gamma t}$ for all $\gamma = 1,\dots,G$ and $t = 1,\dots,T$. We believe that this factor-analytic model for the covariates $x_{it}$ is rather flexible as it allows for individual-specific correlations between covariates and group-specific time effects. In Appendix (ref), we also provide an example in terms of agricultural production functions and environmental economics to motivate equation ((ref)).

Our post-spectral estimator consists of three steps. In the first step, we carry out preliminary estimation of $\beta$. To do so, we prove that as long as the data-generating process is given by equations (ref) and (ref), there exists a convex quadratic function $f\colon \mathbb R^d\to\mathbb R$ such that (i) its unique minimum is achieved by $\beta$ and (ii) for each value of $b\in\mathbb R^d$, the value of $f(b)$ can be consistently estimated by the sum of $2G M+2$ largest in absolute value eigenvalues of a certain matrix. We then demonstrate that this function and its consistent estimator can be used to construct an estimator of $\beta$ that is both consistent and computationally simple. This estimator, which we call the {\em spectral estimator}, may have slow rate of convergence if $T$ is growing rather slowly, and so we proceed to the second and the third steps. In the second step, using the preliminary spectral estimator of $\beta$ obtained in the first step, we carry out classification of units $i=1,\dots,N$ into groups $\gamma = 1,\dots,G$. Importantly, our classification algorithm, which is a version of spectral clustering method (VW04, M18, V18, LZZ21), is fast and does not require solving any non-convex optimization problems. In the third step, we obtain the {\em post-spectral estimator} of $\beta$ by performing OLS-type estimation pooling all units within the same group together.

We prove that our post-spectral estimator is generally consistent and has particularly attractive properties under the assumption of well-separated groups, which means that the vectors $(\alpha_{\gamma 1},\dots,\alpha_{\gamma T})'$, $\gamma = 1,\dots,G$, are not too close to each other, and which was also used in BM15.\footnote{We make no assumptions on the distance between vectors $(\alpha_{\gamma 1}^m,\dots,\alpha_{\gamma T}^m)'$ for $m\geq 2$.} Specifically, we show that under this assumption, the classifier constructed in the second step is consistent in the sense that with probability approaching one, any two units are getting classified into the same group if and only if they belong to the same group, and so the post-spectral estimator of $\beta$ is asymptotically equivalent to the pooled-OLS estimator with known group memberships (i.e., the oracle estimator). In turn, the latter is $\sqrt{NT}$ consistent and admits the standard OLS inference. Under the assumption of well-separated groups, our post-spectral estimator thus can be used for testing hypotheses and for constructing confidence intervals for $\beta$ using standard panel-data methods, ignoring the preliminary estimation and classification steps. Inference without the assumption of well-separated groups, however, remains an open (and challenging) question for future work. Theoretical properties of our post-spectral estimator are thus comparable to those of the grouped fixed-effect estimator, with the caveat that we impose the special structure on the data-generating process for the covariates $x_{it}$ given in (ref).

To compare computational properties of the post-spectral estimator to those of the grouped fixed-effect estimator, we note that the latter, as well as many related estimators (AB16,AB17,CLW19,CSS19,ZWZ17), are built around the K-means optimization problem. This optimization problem is known to be NP-hard (see ADHP09), and so it is rather unlikely that there exists a fast algorithm for finding its solutions. Any proposed fast implementation of the aforementioned estimators therefore is likely to fail occasionally. For example, the grouped fixed-effect estimator is defined as the (global) minimizer of the sum of squared residuals over all parameter values and over all partitions of units into groups, and the main algorithm to calculate this minimizer in BM15 proceeds by initializing randomly selected values of parameters $\beta$ and $\{\alpha_{\gamma t}\}_{\gamma,t=1}^{G,T}$ and then alternating between two steps: (1) optimization over the values of group memberships $g_1,\dots,g_N$ given the values of parameters $\beta$ and $\{\alpha_{\gamma t}\}_{\gamma,t=1}^{G,T}$, and (2) optimization over the values of parameters $\beta$ and $\{\alpha_{\gamma t}\}_{\gamma,t=1}^{G,T}$ given the values of group memberships $g_1,\dots,g_N$. Since this procedure at best converges to a local minimum, it is repeated over many different initial values of parameters to find the global minimum, which corresponds to the grouped fixed-effect estimator. As BM15 notes, however, “a prohibitive number of initial values may be needed to obtain reliable solutions.” In addition, it seems never possible to say whether the global minimum has been found, as this in general would require minimizing the sum of squared residuals over parameter values for {\em each} partition of units into groups, and the number of these partitions, $G^N$, is tremendously large even in small samples. BM15 also proposed a few other algorithms to calculate the grouped fixed-effect estimator that tend to perform better in simulations but they are all subject to the same critique: if they are fast, they must fail occasionally. In contrast, our post-spectral estimator is easy to compute and does not suffer from the potential failure problem.

There is also growing literature on non-linear panel data models with group structure of individual-level parameters (HM10, LN12, SSP16, BML17, SJ18, WPS18, GV19, GS20) originated by Hahn and Moon in HM10. This literature is conceptually related to the grouped panel data model (ref) but estimation techniques developed in this literature are very different from those considered here because all aforementioned papers assume that the individual-level parameters are time-independent, and so preliminary consistent estimation of these parameters is possible by performing estimation separately for each unit. The latter is not possible in the model (ref) because individual effects $\alpha_{g_it}$ are varying over time, which creates one of the key challenges in estimating this model.

We note also that the grouped panel data model (ref) is a special case of a panel data model with interactive fixed effects, corresponding to factor loadings with finite support in the latter model. The methods developed for estimating panel data models with interactive fixed effects can therefore be used to estimate $\beta$ in (ref) as well. To the best of our knowledge, however, most of these methods are either computationally difficult or require conditions that are substantially different from those used in our paper. For example, the estimator in B09 is based on a solution to a non-convex optimization problem, and the estimators in HNR88, ALS01, ALS13 require certain IV-type conditions. Like in our approach, the estimators in P06 also require restricting the data-generating process for the covariates $x_{it}$ but the nature of imposed restrictions is very different. In particular, P06 either imposes a certain rank condition, which can only be satisfied if the dimensionality of $x_{it}$ is sufficiently large, or requires the factor loadings in the equation for covariates to be independent of the factor loadings in the equation for the dependent variable, which in our model would correspond to assuming that $g_i$ in (ref) is different and independent of $g_i$ in (ref),\footnote{P06 claims that his estimators are consistent without requiring independence of the factors but a counter-example is given in WU13, which proves that if the rank condition is not satisfied, then independence is essentially a necessary condition for consistency of the estimators in P06.} thus leading to a random, rather than fixed, effect model. In fact, the only exception in this literature we are aware of is MW19, who developed a computationally relatively simple method and used conditions that are similar to (actually somewhat weaker than) those used in our paper. Their estimator can potentially replace our preliminary spectral estimator. However, the convergence rate of their estimator is only $(T\wedge N)^{-1/2}$, which can be particularly slow if $T$ is growing much slower than $N$, a case of special interest in the grouped panel data model. In contrast, the rate of convergence of our spectral estimator is $(T\wedge N)^{-1}$, which seems much more acceptable for the preliminary estimation. Indeed, we find via simulations reported below that our spectral estimator leads to much better results in reasonably large samples. Finally, AWZ22 has recently developed a method for debiasing estimators with the slow convergence rate $(T\wedge N)^{-1/2}$ yielding estimators with the fast convergence rate $(T\wedge N)^{-1}$. An advantage of our estimator here is that we obtain the fast rate in just one step, instead of two steps, which again may be preferable when $T$ is relatively small, so that the original estimator to be debiased is rather imprecise and is hard to debias.

In addition, we study three extensions of the model (ref). First, we consider a dynamic version of the model, where lagged values of $y_{it}$ appear on the right-hand side of (ref). We demonstrate that our spectral and post-spectral estimators work for this model too, as long as the number of factors $M$ is appropriately modified. We allow for both pre-determined and exogenous covariates in this model. Second, we consider a high-dimensional version of the model (ref), where the number of covariates $d$ is large, potentially much larger than $NT$, but the vector of coefficients $\beta$ is sparse in the sense that it has relatively few non-zero components. We demonstrate how to modify the spectral estimator via $\ell^1$-penalization to obtain a computationally simple and consistent estimator of $\beta$ in this case. Third, we consider an interactive fixed-effect panel data model, where $\alpha_{g_it}$ in (ref) is replaced by $\kappa_i'\phi_t$, with $\phi_t$ being a vector of factors and $\kappa_i$ being a vector of factor loadings. We demonstrate that the spectral estimator, with appropriately modified parameters $G$ and $M$, is consistent in this model with the convergence rate $(T\wedge N)^{-1}$. Being computationally simple, our spectral estimator thus can serve as an alternative to existing estimators in the literature on interactive fixed-effect panel data models. Note, however, that in all three extensions, we maintain a version of (ref).

The rest of the paper is organized as follows. In the next section, we discuss details of implementation of our spectral and post-spectral estimators. In Section (ref), we state their asymptotic properties. In Section (ref), we provide the extensions of the baseline model. In Section (ref), we discuss results of a small-scale Monte Carlo simulation study that shed some light on finite-sample properties of our estimators. In Section (ref), we present main proofs. In Appendix (ref), we collect some technical lemmas that are useful for the proofs of our main results. In Appendix (ref), we present remaining proofs. In Appendix (ref), we describe a method for calculating eigenvalues of large matrices, which may be needed for implementing our estimators. In Appendix (ref), we provide some details on the assumptions of the dynamic model extension. In Appendix (ref), we discuss an example motivating equation (ref).

Estimation

Our proposed estimation procedure consists of three steps. The first step is preliminary consistent estimation of $\beta$, which is based on the spectral analysis of certain matrices and gives the spectral estimator. The second step is classification of units into groups. The third step is pooled-OLS estimation of $\beta$ on classified units, which gives the post-spectral estimator. Under the assumption of well-separated groups, the post-spectral estimator will be $\sqrt{NT}$-consistent and asymptotically normal with mean zero, making inference based on this estimator straightforward.

Spectral Estimator

For all $b\in\mathbb R^d$, let $A^b$ be an $N\times N$ matrix whose $(i,j)$-th element is

Since $A^b$ is an $N\times N$ symmetric matrix, it has $N$ real eigenvalues. Let $\lambda^b_1,\dots,\lambda^b_{2GM+2}$ be its $2G M+2$ largest in absolute value eigenvalues. We will show below that under mild conditions,

where $\Sigma$ is a $d\times d$ symmetric positive definite matrix, $S$ is a $d\times 1$ vector, $L$ is a scalar, and, importantly, $\beta$ is the unique minimizer of the function $b\mapsto f(b) = b'\Sigma b + S'b + L$. We therefore define our spectral estimator as

where $\hat\Sigma$, $\hat S$, and $\hat L$ are estimators of $\Sigma$, $S$, and $L$, respectively, to be constructed below. By the first-order conditions, the estimator $\tilde\beta$ can be equivalently defined as $$ \tilde\beta = -\hat\Sigma^{-1}\hat S/2. $$ We will prove in the next section that $\tilde\beta\to_P \beta$.

Next, we discuss estimators $\hat\Sigma$, $\hat S$, and $\hat L$. For brevity of notation, denote

so that by (ref),

Also, for all $k=1,\dots,d$, let $e_k = (0,\dots,0,1,0,\dots,0)'$ be the $d\times 1$ vector with 1 in the $k$-th position and 0 in all other positions, and let $0_d = (0,\dots,0)'$ be the $d\times 1$ vector with 0 in all positions. Since (ref) implies $\hat f(0_d) = L + o_P(1)$, we set $ \hat L = \hat f(0_d). $ Further, since $\hat f(e_k) - \hat f(-e_k) \to_P 2S_k$, we set $$ \hat S_k = \frac{\hat f(e_k) - \hat f(-e_k)}{2},\quad\text{for all }k=1,\dots,d, $$ and $\hat S = (\hat S_1,\dots,\hat S_d)'$. Finally, since $\hat f(e_k) + \hat f(-e_k) = 2(\Sigma_{k k} + L) + o_P(1)$, we set $$ \hat \Sigma_{k k} = \frac{\hat f(e_k) + \hat f(-e_k)}{2} - \hat L,\quad\text{for all }k=1,\dots,d, $$ and since $\hat f(e_{k} + e_{l}) = \Sigma_{k k} + \Sigma_{l l} + 2\Sigma_{k l} + S_{k} + S_{l} + L + o_p(1)$, we set $$ \hat \Sigma_{k l} = \hat \Sigma_{l k} = \frac{\hat f(e_{k} + e_{l}) - \hat\Sigma_{k k} - \hat\Sigma_{l l} - \hat S_{k} - \hat S_{l} - \hat L}{2}, $$ for all $k,l=1,\dots,d, \ k> l$, and let $\hat \Sigma$ be the matrix whose $(k,l)$-th component is $\hat \Sigma_{k l}$.\footnote{Note that the quality of the estimators $\hat\Sigma$, $\hat S$, and $\hat L$ could potentially be improved by exploiting additional values of the vector $b$ but we leave the question of optimal estimation for future work.}

Under result (ref), the estimators $\hat\Sigma$ and $\hat S$ are consistent for $\Sigma$ and $S$, respectively, and so $\tilde\beta = -\hat\Sigma^{-1}\hat S/2 \to_P -\Sigma^{-1}S/2 = \beta$, as long as $\Sigma$ is invertible, which is the case under mild conditions. The bulk of our derivations in the next section will thus be related to proving (ref).

Before we move on, however, we note that the $N\times N$ matrices $A^b$ may be rather large, and the reader might wonder how much time it takes to calculate their eigenvalues $\lambda^b_1,\dots,\lambda^b_{2GM+2}$. Fortunately, there exists a class of fast randomized algorithms that allow to calculate these eigenvalues arbitrarily well; see Appendix (ref) for details.

Classifier

To classify units into groups, we will use a version of the spectral clustering method.\footnote{Note that the special structure on the data-generating process for the covariates $x_{it}$ given in (ref) was used for the construction of the spectral estimator only. This structure has no role for our classifier and for the post-spectral estimator described below.} For reasons to be explained in Remark (ref) in the next section, we will also rely on sample splitting. To this end, let $h_1,\dots,h_N$ be i.i.d. random variables that are independent of the data and that are taking values $0$ and $1$, each with probability 1/2. We split all cross-sectional units $i=1,\dots,N$ into two subsamples, $\mathcal I_0 = \{i =1,\dots,N\colon h_i=1\}$ and $\mathcal I_1 = \{i =1,\dots,N\colon h_i=0\}$. Further, for $i=1,\dots,N$, denote $y_i = (y_{i1},\dots,y_{iT})'$ and $x_i = (x_{i1},\dots,x_{iT})'$. Also, for $h=0,1$, let $\tilde\beta^h$ be the spectral estimator calculated using the subsample $\mathcal I_h$ and let $\hat B^h$ be a $T\times T$ matrix given by

Since $\hat B^h$ is a $T\times T$ symmetric positive definite matrix, it has $T$ non-negative eigenvalues and $T$ corresponding orthonormal eigenvectors. Let $\hat F_h$ be a $T\times G$ matrix whose columns are orthonormal eigenvectors corresponding to $G$ largest eigenvalues of the matrix $\hat B^h$. Moreover, for all $i=1,\dots,N$, let $\hat A_i$ be a $T\times1$ vector defined by

Intuitively, the vectors $\hat A_1,\dots,\hat A_N$ estimate the vectors $\alpha_{g_1},\dots,\alpha_{g_N}$, where we denoted $\alpha_{\gamma} = (\alpha_{\gamma 1},\dots,\alpha_{\gamma T})'$ for all $\gamma=1,\dots,G$. We therefore classify units $i=1,\dots,N$ into $G$ groups using these vectors. To do so, fix a tuning parameter $\lambda > 0$, to be chosen below, and consider the following algorithm:

{\bf Classification Algorithm.}

{\em Step 1:} set $\mathcal A_1 = \{1\}$, $m = 1$, and $i = 1$;

{\em Step 2:} replace $i$ by $i + 1$;

{\em Step 3:} if $i = N + 1$, stop the algorithm;

{\em Step 4:} set $\mathcal C_i = \{\gamma=1,\dots,m\colon \|\hat A_i - |\mathcal A_{\gamma}|^{-1}\sum_{l\in\mathcal A_{\gamma}}\hat A_l\| \leq \lambda\}$;

{\em Step 5:} if $\mathcal C_i$ is empty, replace $m$ by $m+1$, set $\mathcal A_m = \{i\}$, and go to Step 2;

{\em Step 6:} if $\mathcal C_i$ is not empty, replace $\mathcal A_{\gamma}$ by $\mathcal A_{\gamma}\cup \{i\}$ for $\gamma=\min\mathcal C_i$ and go to Step 2.

This algorithm creates $m$ groups $\mathcal A_1,\dots,\mathcal A_m$, with the number of groups $m$ depending on $\lambda$, so that $m = m(\lambda)$. Clearly, $\lambda\mapsto m(\lambda)$ is a right-continuous function, and so $$ \hat\lambda = \min\Big\{\lambda>0\colon m(\lambda)\leq G\Big\} $$ is well-defined. (In practice, $\hat\lambda$ can be calculated using the values of $m(\lambda)$ on a fine grid.) We classify units $i=1,\dots,N$ into $G$ groups using this algorithm with $\lambda = \hat\lambda$. The result of the algorithm is then $m(\hat\lambda)\leq G$ groups $\mathcal A_1,\dots,\mathcal A_{m(\hat\lambda)}$, and for all $i=1,\dots,N$, there exists a unique $\gamma = \gamma(i)\in\{1,\dots,m(\hat\lambda)\}$ such that $i\in\mathcal A_{\gamma}$. We set $$ \hat g_i = \gamma(i),\quad\text{for all }i=1,\dots,N, $$ and $\hat g = (\hat g_1,\dots,\hat g_N)'$. Note that this classifier can occasionally lead to less than $G$ groups, which happens when $m(\hat\lambda) < G$, but we will show in the next section that it is consistent in the sense that

under the assumption of well-separated groups.

Post-Spectral Estimator

Once we have classified units into groups, estimation of $\beta$ is straightforward. In particular, we rely upon a pooled-OLS estimator:

where $\mathcal B$ is a parameter space for the vector $\beta$, and $\mathcal A_{G,T}$ is a parameter space for the matrix $\{\alpha_{\gamma t}\}_{\gamma,t=1}^{G,T}$. We refer to $\hat\beta$ as the post-spectral estimator for grouped panel data models. We will show in the next section that under the assumption of well-separated groups, this estimator is asymptotically equivalent to the estimator based on correct classification,

and thus to the grouped fixed-effect estimator of BM15. Hence, the standard OLS inference ignoring group classification applies.

Asymptotic Theory

In this section, we derive asymptotic properties of the procedures described above. For convenience, we do so in three separate subsections: spectral estimator, classifier, and post-spectral estimator.

Throughout the rest of the paper, we assume that membership variables $g_i$, group-specific time effects $\alpha_{\gamma t}^m$ and individual specific vectors of coefficients $\rho_{im}$ are non-stochastic, i.e. our analysis is conditional on these random quantities. Also, given that we set $\alpha_{\gamma t}^1 = \alpha_{\gamma t}$, group-specific time effects $\alpha_{\gamma t}$ are non-stochastic as well. Moreover, we assume that the units $i$ are independent.

Spectral Estimator

Let $\mathcal S^T$ denote the unit sphere in $\mathbb R^T$, i.e. $\mathcal S^T = \{u\in\mathbb R^T\colon \|u\|=1\}$. Also, for any random variable $w$, let $\|w\|_{\psi_2}$ denote the sub-Gaussian norm of $w$, i.e. $$ \|w\|_{\psi_2} = \inf\Big\{\epsilon > 0\colon \mathbb E[\exp(w^2/\epsilon^2)]\leq 2\Big\}; $$ see Section 2.5.2 in V18 on properties of the sub-Gaussian norm.\footnote{Sub-Gaussian norm also often appears in the literature under the name of Orlicz norm.} Intuitively, a random variable has a finite sub-Gaussian norm if the tails of its distribution are not heavier than tails of the Gaussian distribution. For example, every bounded random variable has a finite sub-Gaussian norm. To prove consistency and to derive the rate of convergence of the spectral estimator $\tilde\beta$, we will use the following assumptions.

By Hoeffding's inequality (Proposition 2.6.1 in V18), Assumption (ref)(i) holds if the random variables $v_{it}$ have finite sub-Gaussian norm and are independent across $t$. More generally, due to numerous versions of Hoeffding's inequality for time series data (e.g., see D09,V02), Assumption (ref)(i) holds as long as the dependence of the random variables $v_{it}$ across $t$ is not too strong. Assumption (ref)(ii) is similar to Assumption (ref)(i) but imposes the integrability and time series dependence restrictions on $z_{it}$ instead of $v_{it}$. We admit that the assumption of random variables having finite sub-Gaussian norm may be somewhat strong but we emphasize that a version of Theorem (ref) below, with slower rates, can be derived under weaker integrability assumptions. We have chosen to work with Assumption (ref) in order to minimize technicalities of our analysis.

Since $\mathbb E[v_{it}] = 0$, $\mathbb E[v_{it} x_{it}] = 0_d$, and we assume that $\alpha_{g_it}^m$ and $\rho_{im}$ are non-stochastic, it follows from (ref) that $\mathbb E[v_{it}z_{it}] = 0_d$. Hence, Assumption (ref)(i) is a quantitative law of large numbers for the random vector $(NT)^{-1}\sum_{i=1}^N\sum_{t=1}^T v_{it}z_{it}$. Similarly, Assumption (ref)(ii) is a quantitative law of large numbers for the random matrix $(NT)^{-1}\sum_{i=1}^N\sum_{t=1}^T z_{it}z_{it}'$. Assumption (ref)(ii) also imposes the constraint that the probability limit of this matrix is positive-definite, which is an identification condition.

Classifier

For all $\gamma = 1,\dots,G$, let $N_\gamma = 1\{g_i = \gamma\}$ be the number of units $i$ within group $\gamma$. To prove consistency of the classifier $\hat g$, we will use the following assumptions.

Assumption (ref) is self-explanatory. Assumption (ref) means that the groups are well-separated in the sense that the vectors of group-specific time effects, $(\alpha_{\gamma 1},\dots,\alpha_{\gamma T})'$ for $\gamma = 1,\dots,G$, are not too close to each other. We will use this assumption to prove consistency of the classifier $\hat g$ and to derive asymptotic normality of the post-spectral estimator $\hat\beta$ but we will not use it to prove consistency of the post-spectral estimator $\hat\beta$. Note, however, that if groups are not well-separated, it is possible that the spectral estimator actually outperforms the post-spectral one. Also, note that Assumption (ref) is sufficient for consistency of the classifier $\hat g$ but by no means necessary. In particular, using more sophisticated arguments as in LZZ21 and stronger conditions on the noise variables $v_{it}$ (i.e. isotropic Gaussianity, which means that the random variables $v_{it}$ are i.i.d. centered Gaussian), one can replace Assumption (ref) by a much weaker condition $\sum_{t=1}^T(\alpha_{\gamma_1t} - \alpha_{\gamma_2t})^2\geq c \log N$ for all $\gamma_1,\gamma_2=1,\dots,G$ such that $\gamma_1\neq\gamma_2$ and a suitable constant $c>0$.

Assumption (ref) requires that each group $\gamma = 1,\dots,G$ constitutes a non-trivial fraction of all units. If we were to assume random group assignment, where each unit is assigned to group $\gamma$ with probability $\mathrm{p}_{\gamma}>0$, so that $\sum_{\gamma=1}^G\mathrm p_{\gamma}=1$, and units are assigned independently, this assumption would be satisfied with probability approaching one as $N\to\infty$.

Assumption (ref) specifies how fast $T$ and $N$ are required to grow relative to each other in the asymptotics. The most important observation here is that we allow $T$ to be much smaller than $N$, which is the main case of interest for grouped panel data models; see BM15.

Post-Spectral Estimator

In this subsection, we present two results on our post-spectral estimator $\hat\beta$. First, we show that this estimator is generally consistent. Second, we show that under the assumption of well-separated groups, Assumption (ref), this estimator is $\sqrt{NT}$-consistent and has simple asymptotic distribution.

For all $\nu = (\nu_1,\dots,\nu_N)'\in\{1,\dots,G\}^N$ and $\gamma_1,\gamma_2 = 1,\dots,G$, denote $$ \mathcal I(\nu,\gamma_1,\gamma_2) = \Big\{i=1,\dots,N\colon \nu_i =\gamma_1\text{ and }g_i = \gamma_2\Big\} $$ and $$ \bar x_{\nu,\gamma_1,\gamma_2,t} = \frac{1}{|\mathcal I(\nu,\gamma_1,\gamma_2)|}\sum_{i\in\mathcal I(\nu,\gamma_1,\gamma_2)}x_{it},\quad\text{for all }t=1,\dots,T. $$ To derive consistency of the post-spectral estimator $\hat\beta$, we will use the following conditions:

Assumption (ref) requires that covariates $x_{it}$ have sufficient within-group variation over time and across units. This assumption was used in BM15 as well. Assumption (ref) is a standard compactness condition used in the statistical analysis of non-linear models.

Finally, we prove $\sqrt{NT}$-consistency and asymptotic normality of the post-spectral estimator $\hat\beta$. To do so, denote $$ \bar x^{\gamma, t} = \frac{1}{N_{\gamma}}\sum_{i\colon g_i = \gamma}x_{it},\quad\text{for all }\gamma = 1,\dots,G, \ t = 1,\dots, T $$ and $$ \check x_{it} = x_{it} - \bar x^{g_i, t},\quad\text{for all }i=1,\dots,N, \ t = 1,\dots,T. $$ We will use the following condition:

This assumption is similar to the corresponding assumptions in BM15.

Extensions

In this section, we consider three extensions of the model we studied above. The first extension is concerned with a dynamic version of the model, i.e. a model that allows for lagged values of $y_{it}$ on the right-hand side of equation (ref). The second extension is concerned with a high-dimensional version of the model, i.e. a model with high-dimensional $\beta$. The third extension is concerned with an interactive fixed-effect model.

Dynamic Model

Consider a dynamic grouped panel data model

where $x_{it}$ is a $d\times 1$ vector of pre-determined covariates. As before, assume also that (ref) is satisfied. This model is different from the model we studied above as we now allow for the lagged dependent variable on the right-hand side of (ref). In this section, we explain what changes one has to carry out in the spectral and post-spectral estimators to estimate parameters $\theta$ and $\beta$ of this model. We will assume throughout that $|\theta|<1$ since otherwise the extension seems difficult. The results below can be easily extended to allow for additional lagged values of $y_{it}$ on the right-hand side of (ref), i.e. $y_{it-2}$, $y_{it-3}$, etc.

To motivate our approach, note that iterating (ref) and substituting (ref) yields

where we denoted $$ \rho_{im}^y = 1\{m=1\} + \rho_{im}'\beta,\quad\text{for all }i=1,\dots,N , m=1,\dots,M, $$ $$ \alpha_{\gamma t}^{m,y} = \sum_{r=0}^{t-2}\theta^r \alpha_{\gamma t-r-1}^m,\quad\text{for all }t=1,\dots,T, \ \gamma = 1,\dots,G, \ m=1,\dots,M, $$ $$ z_{it}^y = \theta^{t-1}y_{i0} + \sum_{r=0}^{t-2}\theta^r(z_{it-r-1}'\beta + v_{it-r-1}),\quad\text{for all }i=1,\dots,N, \ t=1,\dots,T, $$ where the sum $\sum_{r=0}^{-1}$ is treated as zero. Hence, we can write $$ y_{it} = \mathring{x}_{it}'\mathring\beta + \alpha_{g_it}+v_{it}, $$ where we denoted $\mathring x_{it} = (y_{it-1},x_{it}')'$ and $\mathring\beta = (\theta,\beta')'$ and the vector of covariates $\mathring x_{it}$ satisfies $$ \mathring x_{it} = \sum_{m=1}^{2M}\mathring\rho_{im}\mathring\alpha_{g_it}^m + \mathring z_{it},\quad\text{for all }i=1,\dots,N, \ t=1,\dots,T, $$ where $$ \mathring\rho_{im} = (0,\rho_{im}')'1\{m\leq M\} + (\rho_{im-M}^y,0_d')'1\{m>M\}, $$ $$ \mathring \alpha_{g_it}^m = \alpha_{g_it}^m1\{m\leq M\} + \alpha_{g_it}^{m-M,y}1\{m>M\}, \text{ and }\mathring z_{it} = (z_{it}^y,z_{it}')'. $$ Thus, the dynamic model considered here reduces to the model studied in Sections (ref) and (ref) with $x_{it}$, $\beta$, $\rho_{im}$, $\alpha_{\gamma t}^m$, $z_{it}$, and $M$ replaced by $\mathring x_{it}$, $\mathring \beta$, $\mathring\rho_{im}$, $\mathring\alpha_{\gamma t}^m$, $\mathring z_{it}$, and $2M$, respectively. Therefore, the parameters of the dynamic model, $\mathring\beta = (\theta,\beta')'$, can be estimated by the spectral and post-spectral estimators as in Section (ref) with $M$, $x_{it}$, $d$, and $\mathcal B$ replaced by $2M$, $\mathring x_{it}$, $d+1$, and $\mathring{\mathcal B}$, respectively, where $\mathring{\mathcal B}$ is a parameter space for $\mathring\beta$, e.g. $\mathring{\mathcal B} = [-1+\delta,1-\delta]\times\mathcal B$ for some $\delta\in(0,1)$.

To state the formal results, let Assumptions 4.1--4.12 be the same as Assumptions 3.1--3.12 with $z_{it}$, $x_{it}$, $d$, $\Sigma$, and $\mathcal B$ replaced by $\mathring z_{it}$, $\mathring x_{it}$, $d+1$, $\mathring\Sigma$, and $\mathring{\mathcal B}$, respectively.

High-Dimensional Model

Consider a high-dimensional grouped panel data model $$ y_{it} = x_{it}'\beta + \alpha_{g_it} +v_{it},\quad\text{for all }i=1,\dots,N, \ t=1,\dots,T, $$ where $x_{it}$ is a $d\times1$ vector of covariates and $d$ can be large, potentially much larger than $NT$, but the vector of coefficients $\beta$ is sparse, i.e. $s = \|\beta\|_0 = \sum_{k=1}^d 1\{\beta_k\neq 0\} $ is relatively small (in the sense to be made precise later). As before, assume also that (ref) is satisfied. Estimating this model requires introducing penalized versions of the spectral and post-spectral estimators. For brevity, however, we focus here on the penalized spectral estimator only, as deriving results for the penalized post-spectral estimator requires taking care of a lot of technicalities but does not seem to bring any new insight.\footnote{Under consistent classification (ref), the $\ell^1$-penalized post-spectral estimator will be similar to the usual Lasso estimator, with comparable properties.}

To define the penalized spectral estimator, let $\hat S$ and $\hat\Sigma$ be the same $d\times1$ vector and $d\times d$ matrix as those appearing in Section (ref). Note that calculating these quantities requires only $O(d^2)$ operations, and thus is computationally rather simple. We then define the penalized spectral estimator as $$ \hat\beta_{\lambda} = \arg\min_{b\in\mathbb R^d}\Big\{ b'\hat\Sigma b + \hat S'b +\lambda\|b\|_1 \Big\}, $$ where $\lambda>0$ is a penalty parameter and $\|b\|_1 = \sum_{k=1}^d |b_k|$ denotes the $\ell^1$-norm of the vector $b=(b_1,\dots,b_d)'$. The optimization problem here is convex and can be carried out using standard software.

To analyze this estimator, we are going to rely on the triangular array asymptotics, where the model, as well as the dimension $d$ of the vector of covariates, are allowed to depend on $N$ and $T$ but for brevity of notation, we keep this dependence implicit. Also, we now have to modify Assumption (ref). Indeed, Assumption (ref)(i) is impossible to satisfy if $d$ is growing together with $N$ and $T$ because of the $\ell^2$-norm appearing in the assumption. Assumption (ref)(ii) also has to be modified as the concept $O_P(\cdot)$ is not well-defined when applied to matrices of growing dimensions. However, since the required modification are only used in this subsection, we spell them directly in the statement of the theorem. In addition, to state the formal result, for any matrix $A$, we will use $\|A\|_{\infty}$ to denote its $\ell^{\infty}$-norm, i.e. the maximum of absolute values of components of $A$.

Interactive Fixed-Effect Model

Consider an interactive fixed-effect panel data model

where $\phi_t$ is a $J\times 1$ vector of factors and $\kappa_i$ is a $J\times 1$ vector of factor loadings. B09 developed an OLS-type interactive fixed-effect estimator of $\beta$ in this model. The proposed estimator, however, is a solution to a non-convex optimization problem, and like in the case of the grouped fixed-effect estimator, it may be difficult to find (ensure that we found) the global solution to this optimization problem. To fix this issue, MW19 developed an estimator of $\beta$ based on the nuclear-norm minimization. Their estimator solves a convex optimization problem and so is computationally rather simple. However, the convergence rate of their estimator is only $(T\wedge N)^{-1/2}$, which can be rather slow, especially if $T\ll N$ or $N\ll T$. It is therefore of interest to look for alternative estimators of $\beta$ in this model.

In this subsection, we demonstrate that a trivial modification of our spectral estimator can be used to estimate $\beta$ in this model with the rate $(T\wedge N)^{-1}$ as long as we assume that the covariates $x_{it}$ satisfy the factor model as well, namely

where $\phi_t$ is a $J\times 1$ vector of factors and $\omega_i$ is a $J\times d$ matrix of factor loadings. Note here that it is without loss of generality to assume that the same factors $\phi_t$ appear both in (ref) and in (ref) as we can always merge the factors from two equations.

Let $\tilde\beta_J$ be the same spectral estimator $\tilde\beta$ as that defined in Section (ref) with $J$ replacing $GM$. We then have the following result.

Monte Carlo Simulation Study

In this section, we present results of a Monte Carlo simulation study. We compare the performance in finite samples of the spectral (S) and post-spectral (P-S) estimators to several natural alternatives. Specifically, we consider the grouped fixed effect (GFE) estimator proposed in BM15, the least squares (LS) estimator proposed in B09, and the penalized nuclear norm (Pen NN) estimator proposed in MW19. The comparison with LS and Pen NN estimators is pertinent because equation ((ref)) can be seen as a particular case of an interactive fixed effects model where the factors are the group trends, $\alpha_{gt}$, and the loadings are restricted to be of the form $(0,\ldots,1,\ldots,0)$, with 1 in the $g_i^{th}$ position.

Computation of the GFE estimator is not trivial. Even with $N=100$ and $G=2$, it is hard to be sure that any given algorithm produces the global minimum of the GFE optimization problem since this would require separately considering $G^N$ partitions of units $i=1,\dots,N$ into $G$ groups. Having in mind this issue, we calculate the GFE estimator using Algorithm 1 in BM15 with 500 random initial values, which is in a nutshell the Lloyd algorithm with covariates. Given a classification, this algorithm obtains the regression coefficients and the group-specific trends by minimizing OLS, and given the parameters, units are classified into groups according to the smallest individual-specific residual. We generate the initial values from the standard normal distribution in $\mathbb R^{2 + G\times T}$, where the first two and the last $G\times T$ components correspond to $\beta = (\beta_1,\beta_2)'$ and $\{\alpha_{\gamma t}\}_{\gamma,t=1}^{G,T}$, respectively.

We also consider an infeasible GFE estimator, which uses only one initial value but this value is chosen to be the oracle estimator, i.e. the pooled OLS estimator based on the true partition of units into groups. Using this value substantially increases the chances to find the global minimum of the GFE optimization problem as we expect the solution to be near the oracle estimator. In what follows, we refer to this infeasible GFE estimator as I-GFE.

Computation of the LS estimator can also be problematic, as its objective function is prone to local minima as well. To deal with this problem, like in the case of the GFE estimator, we use 500 random initial values and choose the one that gives the best value of the LS criterion function.

Computation of the Pen-NN estimator is relatively straightforward: it minimizes a convex objective function. To make the comparison with the other estimators fair, however, we choose a penalty parameter for this estimator so that the total number of generated groups is equal to the true value, $G$, instead of using the data-driven procedure proposed by the authors. Also, instead of presenting results for the plain Pen-NN estimator, following the original paper MW19, we actually present results for the LS estimator that uses the Pen-NN estimator as initial value, as adding an extra step of least squares minimization is supposed to make the estimator more precise. With some abuse of notation, we still refer to this estimator as Pen-NN.\footnote{MW19 proposed several other estimators as well but we have chosen to present results for the Pen-NN estimator only, as it dominates the other estimators in our simulations.}

Finally, for all data-generating processes, we also report results for the oracle estimator, which is the pooled OLS estimator based on the true partition of units into groups. Although we do not discuss these results, an interested reader can use them as a benchmark for the results on the other estimators.

Next, we discuss the data-generating processes. We generate the data $(x_{it}',y_{it})'$, $i=1,\dots,N$ and $t=1,\dots,T$, according to equations (ref) and (ref), where we set $d=2$, so that the model contains two covariates. Depending on the experiment, we set $N=100$, $200$, or $400$ and $T=20$, $50$, or $100$. We also set $G=2$ or $7$ and $M=1$ or $2$. The random variables $v_{it}$ and components of the random vectors $z_{it}$ are generated from a truncated standard normal distribution, $Z1\{|Z|\leq C\}$, where $Z\sim N(0,1)$ and $C=20$, independently of each other and across indices. The random variables $\alpha_{\gamma t}^m$ are generated from a truncated normal distribution as well but with different variance values, $Z1\{|Z|\leq C\}$, where $Z\sim N(0,\sigma^2)$, independently across indices and of all other random variables. Depending on the experiment, we set $\sigma^2$ equal to either 1 or 4. We have checked that results are robust with respect to different specifications of the truncation constant $C$. Further, we set the number of observations $i$ within each group $\gamma$ to be the same, except for the last group, which contains more observations if $G$ does not divide $N$. We have also checked with uneven units across groups and results are unchanged.

Throughout all experiments, we set $\beta = (\beta_1,\beta_2)'$, where $\beta_1=-1$ and $\beta_2 = .8$ but we note that the particular choice of these values does not seem to matter much for the simulation results. Also, when $M=1$, we set $\rho_{i1} = (\rho_{i11},\rho_{i12})'$ with $\rho_{i11} = \varrho + Z_{i1}1\{|Z_{i1}|\leq C\}$ and $\rho_{i12} = \varrho + Z_{i2}1\{|Z_{i2}|\leq C\}$, where $Z_{i1}$ and $Z_{i2}$ are independent of each other, across $i$, and of all other random variables in the data-generating process, and where $\varrho = 3$. Here, $\varrho$ can be understood as a measure of endogeneity in the model as it govern the correlation between covariates and grouped-specific time effects. When $M=2$, we set $\rho_{i11} = \varrho + Z_{i11}1\{|Z_{i11}|\leq C\}$, $\rho_{i21} = 1 + Z_{i21}1\{|Z_{i21}|\leq C\}$, and $\rho_{i12} = \varrho + Z_{i12}1\{|Z_{i12}|\leq C\}$ and $\rho_{i22} = Z_{i22}1\{|Z_{i22}|\leq C\}$, where random variables $Z_{imj}$ for $m=1,2$ have standard normal distribution and are independent of each other, across indices, and of all other random variables in the data-generating process.

We present results in Tables 1--4 at the end of Supplementary Materials. For each combination of $N$ and $T$, these tables give the mean absolute error (MAE) for S, P-S, LS, Pen-NN, I-GFE, GFE, and oracle estimators. In addition, the tables give the fraction of misclassified units based on our classification algorithm (which is based on the spectral estimator, and which is the second step in constructing our post-spectral estimator) and the fraction of misclassified units in the construction of the GFE estimator. Each table presents results in the upper panel for $G =2$ and in the lower panel for $G = 7$. Tables 1 and 2 correspond to $M=1$ and Tables 3 and 4 correspond to $M=2$. In addition, Tables 1 and 3 correspond to $\sigma^2=1$, and Tables 2 and 4 correspond to $\sigma^2=4$. When increasing the variance $\sigma^2$, we are effectively making the groups more separated, and hence group classification improves. Both the MAE and the fraction of misclassified units are calculated as the average over 50 simulations.

We now describe the results. We start with the behavior of the S and P-S estimators. In most cases, the P-S estimator outperforms the S estimator. Hence, it seems that doing the "post" step delivers a better estimator. Also, the P-S estimator tends to be comparable to the oracle estimator across all tables as long as $N$ and $T$ increase sufficiently. In addition, while the P-S estimator seems to have more trouble when the separation of the groups is less, i.e. in Tables 1 and 3, it still quickly achieves the oracle for moderate values of $N$ and $T$. Moreover, when $G=2$, the P-S estimator seems to achieve the oracle faster than when $G = 7$, but even for small values of $N$ and $T$, the P-S estimator is a robust reliable estimator.

Let us now turn to the comparison with the grouped fixed effect estimators. Not surprisingly, the infeasible I-GFE estimator is often comparable to the oracle. However, the performance of the feasible GFE estimator varies substantially depending on the number and separation of the groups. For instance, in the upper panel of Table 1 ($G=2)$, where separation of groups is moderate, and for moderate values of $T$, the GFE estimator is close to the I-GFE estimator, indicating that the minimum might have been achieved. However, in the lower panel of Table 1 ($G = 7$), it is unlikely that 500 initial values cover a significant part of the space of all possible partitions of $N$ units into $G$ groups, and the performance of the GFE estimator deteriorates. In fact, missclassification for the GFE estimator is on the level of $70\%$ for moderate values of $T$ when $G = 7$. In this case, the P-S estimator has a comparative advantage relative to the GFE estimator. This same pattern is reproduced in Table 2. However, since the groups are more separated, and the r-squared of the model is higher, the GFE estimator improves. Interestingly, the behavior of the GFE estimator when $M=2$ improves relative to $M=1$.

As far the LS estimator is concerned, when the number of groups increases, the finite sample performance of this estimator deteriorates, although it improves as both $N$ and $T$ increase. However, it is only for large values of $T$ that the behavior is comparable to the oracle. At least in this exercise, the P-S estimator outperforms the LS estimator.

Finally, the pen NN estimator is, in general, dominated by the other estimators for our data-generating processes. However, like with the LS estimator, its performance improves as we increase $N$ and $T$. For example, in Table 4, when $G=2$, the estimator is actually comparable to the oracle for large values of $N$ and $T$.

Proofs of Theorems (ref), (ref), and (ref)