Subspace Clustering for Panel Data with Interactive Effects

Subspace Clustering for Panel Data with Interactive Effects

\jyear{2020}

\startabstract{

} \makechaptertitle

Introduction

Panel data, also known as longitudinal data, contain multi-dimensional observations obtained over multiple periods for a sample of individuals. There is evidence to show that unobservable heterogeneity among different individuals in multiple dimensions exists in panel data and hence, suitable statistical models are needed. Among the models for panel data, the interactive fixed-effect model uses the interactive effects that combined individual effects and time effects to reveal the effects of the common factors, which can capture the unobserved information. Panel data model with interactive effects has been widely studied in the literature Pesaran2006,Bai2009.

Due to the number of parameters in the standard fixed-effects model is the same as the number of individuals, the estimation of the fixed-effects may be inaccurate. Therefore, in order to study the individual heterogeneity, the number of parameters in the panel data model needs to be reduced. The grouped panel data model is an effective method to solve this incidental-parameter problem. In a grouped panel data model, individuals in the same group are assumed to have the same effect, which is called group effect, and these group effects can reflect the individual heterogeneity. Grouped panel data models have been studied in the past decade. For example, Lin2012 studied the linear panel data models with time-varying grouped heterogeneity by using the $K$-means clustering algorithm. Bonhomme2015 developed the “grouped fixed-effects" (GFE) estimation method to estimate the model parameters and derived the statistical properties of the GFE estimators when the sample sizes of cross-section ($N$) and length of the time series ($T$) are large. Ando2016 studied the grouped panel data models with unobserved group factor structures, and estimated the model parameters by minimizing the sum of least-squared errors with a shrinkage penalty. Ando2016 also proved the consistency and asymptotic normality of the estimators under large $N$ and large $T$. SuPhillips2016 proposed the C-Lasso method to estimate the parameters in heterogeneous linear panel data models, where the slope parameters are heterogeneous across groups but homogeneous within a group with unknown group membership. Su_Ju2018 considered the penalized principal component estimation method by extending the C-Lasso method to deal with panel data with interactive fixed effects. Su_Ju2018 also assumed the individual slope parameters are heterogeneous across groups but homogeneous within a group.

In those existing studies, the unobserved heterogeneous parts are assumed to live in the same space and they have the ball-shaped (spherical Gaussian) clusters. In many practical applications, a data object often has multiple attributes and many of which may be live in some low dimensional subspaces. For example, for disease detection in newborns, various tests (e.g., blood test and heart rate test) are performed on the newborns and the level of those factors are measured. Each newborn is associated with a vector containing the values of the factors and one can further construct a newborn-factor level matrix in which each row contains the factor levels of a different newborn. Pediatricians want to cluster groups of newborns based on the disease those newborns suffer from. Usually, each disease is correlated with a specific set of factors, which implies that points corresponding to newborns suffering from a given disease in a lower-dimensional subspace Kriegel2009. Therefore, the clustering of newborns based on their specific diseases together with the identification of the relevant factors associated with each disease can be formulated as a subspace clustering problem. The $K$-means algorithm clusters data from around cluster centers to the clustered data in the entire data space and estimates the cluster centers by minimizing the sum of squared distances from the data points to their nearest cluster centers. Therefore, traditional clustering methods such as the $K$-means algorithm may not be meaningful in these cases.

In this paper, we characterize the unobserved effects of the $i$-th unit at time $t$ as $\eta_{it}=\mbox{\boldmath $f$}_t'\mbox{\boldmath $\lambda$}_i$, where $\mbox{\boldmath $f$}_t$ and $\mbox{\boldmath $\lambda$}_i$ are $r \times 1$ dimensional vectors and $\mbox{\boldmath $\lambda$}_i$ $i = 1, \cdots, N$ live in some low-dimensional subspaces. A key feature of the low dimensional subspaces is that it can deal with more general data types, and the factor loadings live in the low dimensional subspaces can encapsulate the main direction of variation. Terada2014 showed that subspace clustering is a more general clustering method which includes the conventional $K$-means clustering method as a special case. A novel approach called least squares subspace clustering (LSSC) is proposed to simultaneously estimate the model parameters and cluster the individual effects by using the least-squares criterion with the subspace clustering principle. In the clustering of individual effects, we treat each individual effect as a vector in a high-dimensional space and then cluster these vectors into some low-dimensional subspaces. The main contributions of this article are highlighted as (i) the grouping is done according to limiting the data points to a specific subspace instead of according to the distance between points, which can better reflect the underlying structure of the data and can be applied to a more general data structure. Furthermore, the consistency of the clustering procedure are proved; (ii) the proposed model allows the covariates to be correlated with the factor structure, and Monte Carlo simulation results show that the LSSC method performs well even when $T$ is small; (iii) the proposed methods allow different groups to share common factors, and different groups of factors can be correlated. The model structure can well capture the spatial structure of individual effects and reflect the cluster structure in real data; (iv) we proposed a model selection criterion to choose the subspace dimensions consistently, which makes the proposed model and methods more general and flexible.

The rest of this paper is organized as follows. In Section 2, we describe the model and states some constraint conditions. In Section 3, we propose an algorithm for estimating the model parameters and subspace clustering simultaneously. In Section 4, we derive the consistency of the subspace clustering and study the asymptotic properties of the estimators. In Section 5, Monte Carlo simulation studies with different settings are used to illustrate the performance of the proposed method and to compare with the GFE estimation method and the method proposed in Bai2009. In Section 6, we discuss some further considerations when the number of subspaces for factors, the dimension of factors and the dimension of subspaces are unknown. For illustrative purposes, the proposed approach is applied to study the relationship between income and democracy in Section 7. Finally, some concluding remarks are provided in Section 8. All the proofs of the theoretical results and additional simulation results are provided in the Supplementary Materials.

Model Description

Let $k$ be the number of subspaces which is unknown but fixed, and $G = \{g_{1}, g_{2}, \ldots, g_{N}\}$ be the grouping of the cross-sectional units into the $k$ subspaces, where the subspace membership variable $g_{i} = j$ denotes the $i$-th unit belongs to the $j$-th subspace with $g_{i} \in \{1, 2, \ldots, k\}$. We further let $N_{j}$ be the number of cross-sectional units within the $j$-th subspace and the total number of units is $N = \sum_{j=1}^{k} N_{j}$. We consider the following panel data model with subspace factor structure

where $y_{it}$ is the respond variable of the $i$-th unit observed at time $t$, $\mbox{\boldmath $x$}_{it}$ is a $p \times 1$ observable vector, $\mbox{\boldmath $\beta$}$ is a $p \times 1$ unknown regression coefficient vector. $\mbox{\boldmath $\lambda$}_{g_{i},i} = (\lambda_{g_{i},i}^{1}, \ldots, \lambda_{g_{i},i}^{r})^{'}$ is a $r \times 1$ factor loading vector that represents the unobserved unit/individual effect for the $i$-th individual. $\mbox{\boldmath $f$}_{g_{i}, t}$ is a $r \times 1$ vector of unobservable subspace-specific pervasive factors that affect the units only in $g_{i}$-th subspace, and $\varepsilon_{i,t}$ is the unit-specific error.

In this paper, we assume that all the factor loadings or individual effects are inside a $r$-dimensional space, which contains some low-dimensional subspaces and these subspaces contain all the individual effects $\mbox{\boldmath $\lambda$}_{i},\;i=1, 2, \ldots, N$. The mathematical notation can be expressed as $\mbox{\boldmath $\lambda$}_{i:g_{i} = j} \subset S_{j}$ and $\displaystyle\bigcup_{j=1}^{k}S_{j} \subset \mathbb{R}^{r}$, where $S_{j}$ is the $j$-th subspace. The covariate $\mbox{\boldmath $x$}_{it}$ can be correlated to $\mbox{\boldmath $\lambda$}_{g_i,i}$ alone or to $\mbox{\boldmath $f$}_{g_i,t}$ alone, or it can be correlated to both $\mbox{\boldmath $\lambda$}_{g_i,i}$ and $\mbox{\boldmath $f$}_{g_i,t}$ simultaneously. Here, $\mbox{\boldmath $x$}_{it}$ can be a nonlinear function of $\mbox{\boldmath $\lambda$}_{g_i,i}$ and $\mbox{\boldmath $f$}_{g_i,t}$. Stacking the observations over $t$, we have $\mbox{\boldmath $F$} = (\mbox{\boldmath $f$}_{1}, \mbox{\boldmath $f$}_{2}, \ldots, \mbox{\boldmath $f$}_{T})$. Furthermore, if we let $\mbox{\boldmath $F$}_{j}$ be the vector of factors for the $j$-th subspace, then we have $\mbox{\boldmath $F$}_{j} = (\mbox{\boldmath $f$}_{g_{i}=j,1}, \mbox{\boldmath $f$}_{g_{i}=j,2}, \ldots, \mbox{\boldmath $f$}_{g_{i}=j,T})^{'}$. Similarly, we have $\mbox{\boldmath $\Lambda$}_{j}=(\mbox{\boldmath $\lambda$}_{j,1},\mbox{\boldmath $\lambda$}_{j,2},\cdots,\mbox{\boldmath $\lambda$}_{j,N_j})^{'}$. We also consider the constraints $\mbox{\boldmath $F$}_{j}^{'}\mbox{\boldmath $F$}_{j}/T = I_{r}(j=1,2, \ldots,k)$, $\;\mbox{\boldmath $\Lambda$}_{j}^{'}\mbox{\boldmath $\Lambda$}_{j}(j=1,2,\ldots,k)$ being diagonal for the issue of identifiability as described in Bai2009, Ando2016 and Stock2002. We aim to estimate the parameters $\mbox{\boldmath $\beta$}$, $\mbox{\boldmath $\Lambda$}_{j}$, $\mbox{\boldmath $F$}_{j}$, $j = 1, 2, \ldots,k$ for each subspace, identify the subspace membership $G=\{g_{1}, g_{2}, \ldots,g_{N}\}$ and the bases $B_{1},\ldots,B_{k}$ of the orthogonal spaces of $S_{1},\ldots,S_{k}$ (denoted as $S_{1}^{\bot},\ldots,S_{k}^{\bot}$) simultaneously.

Parameter Estimation and Clustering

In this section, we propose a method for the estimation of model parameters and subspace clustering simultaneously. The proposed method is a natural generalization and combination of both the least-squares estimation method and the $K$-means clustering algorithm. Thus, the proposed parameter estimation and clustering is innovative compared to those existing approaches because most of the existing statistical methods are based on the $K$-means and least-squares method to cluster and estimate simultaneously.

Estimation Procedure

For a given number of subspaces $k$, the objective function is

where $\mbox{\boldmath $y$}_{i}=(y_{i1},y_{i2},\cdots,y_{iT})$, $\mbox{\boldmath $x$}_{i}=(\mbox{\boldmath $x$}_{i1},\mbox{\boldmath $x$}_{i2},\cdots,\mbox{\boldmath $x$}_{iT})$. The estimator of the vector of model parameters $\{\hat{\mbox{\boldmath $\beta$}},\hat{\mbox{\boldmath $F$}}_{1},\ldots,\hat{\mbox{\boldmath $F$}}_{k},\hat{\mbox{\boldmath $\Lambda$}}_{1},\ldots,\hat{\mbox{\boldmath $\Lambda$}}_{k},\hat{B}_{1},\ldots,\hat{B}_{k},\hat{G} \}$ is defined as

subject to the constraints $\mbox{\boldmath $F$}_{j}^{'}\mbox{\boldmath $F$}_{j}/T = I_{r} (j = 1, \ldots, k)$, $\mbox{\boldmath $\Lambda$}_{j}^{'}\mbox{\boldmath $\Lambda$}_{j} (j = 1, \ldots, k)$ being diagonal, where $\mbox{\boldmath $\Lambda$}_{j}^{'}=(\mbox{\boldmath $\lambda$}_{j,1}$, $\ldots,\mbox{\boldmath $\lambda$}_{j,N_{j}})$ is the $r\times N_{j}$ factor loading matrix $(j=1,\ldots,k)$ for the subspace-specific factors and they live in $k$ different subspaces embedding in the $r$-dimensional space. In Eq. ((ref)), $span\{B_{j}\}$ represents the subspace spanned by the basis $B_{j}$, and the $span\{B_{j}\}^{\bot}$ is the orthogonal subspace of $span\{B_{j}\}$. These constraints and assumptions are needed to ensure the model is identifiable. Here, we aim to estimate the model parameters and to cluster the $\mbox{\boldmath $\lambda$}_{i} \in \mathbb{R}^{r}, i = 1, 2, \ldots, N$ into the $k$ different subspaces simultaneously. Different from the existing classification methods, we approach this challenging problem from a subspace clustering point-of-view. The major idea is to divide the space $\mathbb{R}^{r}$ into several subspaces and project $\mbox{\boldmath $\lambda$}_{i}$ $(i = 1, 2, \ldots, N)$ into the nearest subspace for the classification.

The constrained minimization of the objective function $Q(\mbox{\boldmath $\beta$}, \mbox{\boldmath $F$}_{1}, \ldots, \mbox{\boldmath $F$}_{k}, \mbox{\boldmath $\Lambda$}_{1}, \ldots, {\mbox{\boldmath $\Lambda$}}_{k}, {B}_{1}, \ldots, {B}_{k}, G )$ in Eq. ((ref)) can be obtained by the following iterative algorithm:

Although the least squares objective function is not globally convex Bai2009, from the results of the Monte Carlo simulation studies, the proposed algorithm is robust to the starting value under large $N$ and large $T$. In the numerical experiments, we propose using the least squares method to get the initial value $\mbox{\boldmath $\beta$}^{(0)}$ by ignoring the unobserved group factor structures. We observe in the Monte Carlo simulation studies (Section 5) that the initial value $\mbox{\boldmath $\beta$}^{(0)}$ has good convergence property. In practice, if one has concern that the algorithm may converge to a local optimizer, we suggest using different random starting values, and then select the solution that yields the lowest value of the objective function if those solutions are different.

To evaluate the complexity of the above iterative procedure, we consider the complexity of Steps 2--4 in the above algorithm. In Step 2, we need to obtain the top $r$ singular values of the $T \times N$ data matrix in order to obtain the values $\mbox{\boldmath $F$}$ and $\mbox{\boldmath $\Lambda$}$, which requires a complexity of $O(\delta_{NT}^2\tilde{\delta}_{NT})+O(\delta_{NT}^3)$ where $\delta_{NT}=\min[N,T],\tilde{\delta}_{NT}=\max[N,T]$. Step 3 of the algorithm requires sorting the $\mbox{\boldmath $\lambda$}_i$ according to their projection distances to these subspaces; the computational cost is $O(Nkr)$. In Step 4, the least squares estimation requires a complexity of $O(p^2N)$ because of the inverse operation. Thus, the complexity of the above iterative process is $O(\delta_{NT}^2\tilde{\delta}_{NT})+O(\delta_{NT}^3)+O(Nkr)+O(p^2N)=O(\delta_{NT}^2\tilde{\delta}_{NT})$.

Subspace Clustering for Factor Loadings

In this section, we present the procedure of the subspace segmenting {for the factor loadings $\mbox{\boldmath $\lambda$}_{i},\;i=1,\ldots,N$. We assume that the $i$-th factor loading $\mbox{\boldmath $\lambda$}_{i}\in \mathbb{R}^{r}$ ($i=1,\ldots,N$) is inside $k$ different subspaces, in which the dimension of these $k$ subspaces are $d_{1}, d_{2}, \ldots, d_{k}$, where $\;0 < d_{j} < r, \;j=1,2, \ldots, k$. For the $j$-th subspace $S_{j} \subset \mathbb{R}^{r}$ with dimension $d_{j}$, a basis $B_{j} = [\mbox{\boldmath $b$}_{j1},\ldots,\mbox{\boldmath $b$}_{j,r-d_{j}}]\in \mathbb{R}^{r\times (r-d_{j})}$ is selected for its orthogonal complement $S_{j}^{\perp}$.} Using these notations, we can obtain the following equation for the $j$-th subspace $S_{j}$ and the $i$-th factor loading $\mbox{\boldmath $\lambda$}_i$:

Since $\mbox{\boldmath $\lambda$}_i \in \mathbb{R}^{r}$ belongs to $\cup_{j=1}^{k}S_{j}$ if and only if $( \mbox{\boldmath $\lambda$}_i \in S_{1}) \vee \ldots \vee (\mbox{\boldmath $\lambda$}_i\in S_{k})$, where the notation $\bigvee$ represents the “or" operator. This condition is equivalent to

where the notation $\bigwedge$ represents the “and" operator, and $\sigma$ is a particular choice of the normal vector $\mbox{\boldmath $b$}_{j\sigma(j)}$ from the $j$-th basis $B_{j}$. Note that the right-hand side of Eq. ((ref)) is obtained by exchanging the “and" and “or" operators using De Morgan's laws. Since

which is a homogeneous polynomial of degree $k$ in $r$ variables, we can write each of the polynomials as

where $\mbox{\boldmath $c$}_{k}$ is the vector of polynomial coefficients, $c_{k_{1},k_{2},\ldots,k_{r}}$ is the polynomial coefficient, $\mbox{\boldmath $\lambda$}_i = (\lambda_{i1},\cdots,\lambda_{ir})$, $\mbox{\boldmath $v$}_{k}: \mathbb{R}^{r}\rightarrow \mathbb{R}^{M_{k}(r)}$ is the Veronese map of degree $k$ Fischler1981 which is also known as the polynomial embedding in machine learning defined as $\mbox{\boldmath $v$}_{k}: [\lambda_{i1}, \ldots, \lambda_{ir}]^{T} \mapsto [\ldots,\mbox{\boldmath $\lambda$}_i^{I},\ldots]^{T}$ with $I$ being chosen in the degree-lexicographic order, $\mbox{\boldmath $\lambda$}_i^{I}=\lambda_{i1}^{k_{1}}\lambda_{i2}^{k_{2}}\ldots \lambda_{ir}^{k_{r}}$ and the dimension $M_{k}(r)=C_{k+r-1}^{r-1}$.

Since the polynomial in Eq. ((ref)) can be satisfied by all the factor loadings $\mbox{\boldmath $\lambda$}_{i}, i = 1, 2, \ldots, N$, we can then use these factor loadings to obtain the subspaces. Although the polynomial equations in Eq. ((ref)) are nonlinear in each point $\mbox{\boldmath $\lambda$}_i$, these polynomials are actually linear in the vector of polynomial coefficients $\mbox{\boldmath $c$}_{k}$. Indeed, since each polynomial $p_{k \sigma}(\mbox{\boldmath $\lambda$}_i)=\mbox{\boldmath $c$}_{k}^{T}\mbox{\boldmath $v$}_{k}(\mbox{\boldmath $\lambda$}_i)$ must be satisfied by every data point, we can obtain $\mbox{\boldmath $c$}_{k}^{T}\mbox{\boldmath $v$}_{k}(\mbox{\boldmath $\lambda$}_{i})=0$ for all $i = 1, 2, \ldots, N$.

Suppose that $I_{k}$ is the space of the vector of polynomial coefficients $\mbox{\boldmath $c$}_{k}$ of all the homogeneous polynomial that vanishes in the $k$ subspaces, then the vector of polynomial coefficients of the factorizable polynomial defined in Eq. ((ref)) span a (possibly proper) subspace in $I_{k}$ as $span_{\sigma}\{ p_{k\sigma} \}\subseteq I_{k}$. As every vector $\mbox{\boldmath $c$}_{k}$ in $I_{k}$ represents a polynomial that vanishes on all the data points (on the subspaces), the vector $\mbox{\boldmath $c$}_{k}$ must satisfy the system of linear equations

where $V_{k}(r)\in \mathbb{R}^{M_{k}(r)\times N}$ is the embedded data matrix.Hance, we have the relationship $I_{k} \subseteq null(V_{k}(r))$.

Remark 1: The zero set of each vanishing polynomial $p_{k}(\mbox{\boldmath $\lambda$}_i),\;i=1, 2, \ldots, N$ is a surface in $\mathbb{R}^{r}$, therefore, the derivative of $p_{k}(\mbox{\boldmath $\lambda$}_i)$ at $\mbox{\boldmath $\lambda$}_{i} \in S_{j}$, denoted as $Dp_{k}(\mbox{\boldmath $\lambda$}_i)$, gives a vector normal to the surface. Since a union of subspaces is locally flat, i.e., in a neighborhood of $\mbox{\boldmath $\lambda$}_{i}$ the surface is merely the surface $S_{j}$, then the derivative at $\mbox{\boldmath $\lambda$}_{i}$ lies in the orthogonal complement $S_{j}^{\perp}$ of $S_{j}$. By evaluating the derivatives of all the polynomials in $I_{k}$ at the same point $\mbox{\boldmath $\lambda$}_{i}$, we obtain a set of normal vectors that span the orthogonal complement of $S_{j}$.

Following the results in Theorem 3 of Vidal2005, we can obtain a set of polynomial $p_k(\mbox{\boldmath $\lambda$}_i),i=1\cdots,N$ with coefficients equal to the eigenvectors in the null space of $V_{k}(r)$. By evaluating the derivatives $Dp_k(\mbox{\boldmath $\lambda$}_i)$ at each $\mbox{\boldmath $\lambda$}_i$, $i = 1, \cdots,N$, we can obtain a set of vectors orthogonal to the subspace that the points lie in. Note that the generalized principal component analysis method Vidal2005,Vidal2016 relies on reliable samples per subspace to segment the dataset, however, in the presence of noise, the sample may not be reliable. Here, for each sample, we assume that the sample could be obtained from all the candidate co-dimension classes, and the sample is voted by the dominate vectors of $Dp_k(\mbox{\boldmath $\lambda$}_i)$ as a basis. Finally, the base associated with the highest vote will be used as the normal vectors perpendicular to the subspaces as suggested by Yang2005. After obtaining the orthogonal bases of those subspaces, we can assign $\mbox{\boldmath $\lambda$}_{i}$ to the subspace $j^{*}$, where $j^{*} = \arg \min_{j =1,\ldots,k}|| B_{j}^{T}\mbox{\boldmath $\lambda$}_{i}||$.

Asymptotic Properties

In this section, we characterize the asymptotic properties of the estimators as $N$ and $T$ tend to infinity. We prove that the estimated clustering converges to the corresponding true subspaces under some conditions. In particular, we follow the method in Pollard1981 to establish the consistency of subspace clustering.

Consistency of clustering procedure

The proposed clustering procedure prescribes a criterion for partitioning a set of points into $k$ subspaces. To divide the factor loadings $\mbox{\boldmath $\lambda$}_{1}, \ldots, \mbox{\boldmath $\lambda$}_{N}$ in $\mathbb{R}^{r}$, we first choose $k$ ($k$ is fixed) cluster subspaces $S_{1},S_{2},\ldots,S_{k}$ with dimensions $d_{1}, d_{2}, \ldots, d_{k}$, respectively, that minimize

where $\mbox{\boldmath $\lambda$}_{1}, \ldots, \mbox{\boldmath $\lambda$}_{N}$ can be viewed as $N$ vectors of the sample points and $\Delta( \mbox{\boldmath $\lambda$}_{i},S_{j} )$ is the angle between the vector $\mbox{\boldmath $\lambda$}_{i}$ and the subspace $S_{j}$ which is a value in $[0, \pi/2]$. Let $|| \mbox{\boldmath $\lambda$}_{i} ||=1$ and $\mbox{\boldmath $b$}_{j1},\ldots,\mbox{\boldmath $b$}_{jm}$ be an orthonormal basis of $S_{j}$ and $\theta$ be the angle between $\mbox{\boldmath $\lambda$}_{i}$ and $S_{j}$, then $\sin(\theta)=|| \mbox{\boldmath $\lambda$}_{i}-\sum_{\ell=1}^{m}\mbox{\boldmath $b$}_{j \ell}(\mbox{\boldmath $b$}_{j \ell}^{T}\mbox{\boldmath $\lambda$}_{i}) ||=\sqrt{1-\sum_{\ell=1}^{m}(\mbox{\boldmath $b$}_{ j\ell}^{T}\mbox{\boldmath $\lambda$}_{i})^{2}}$. Thus, clustering by angles is equivalent to clustering by projecting the $\mbox{\boldmath $\lambda$}_{i}$ to the nearest subspace. Here, the function $\phi$ must satisfy some regularity conditions that $\phi$ is continuous and nondecreasing, with $\phi(0) = 0$. For any subspaces $\mbox{\boldmath $x$}$ and $\mbox{\boldmath $y$}$, $\Delta(\mbox{\boldmath $x$}, \mbox{\boldmath $y$}) \in [0,\pi/2]$, therefore, $\phi(\Delta( \mbox{\boldmath $x$},\mbox{\boldmath $y$} ))$ must be in a compact set.

Since $\max\limits_{1\leq i \leq N}||\hat{\mbox{\boldmath $\lambda$}}_i-\mbox{\boldmath $\lambda$}_i||=o_{p}(1)$, we can show that the empirical distribution function converge uniformly to the true distribution function by the strong law of large numbers and the Glivenko-Cantelli theorem, i.e., $$\lim\limits_{N\rightarrow \infty} \sup\limits_{\mbox{\boldmath $x$}\in \mathbb{R}^{r}}|\hat{P}_{N}(\mbox{\boldmath $x$})-P(\mbox{\boldmath $x$})|=0,$$ where $\hat{P}_{N}(\mbox{\boldmath $x$})=\frac{1}{N}\sum\limits_{i=1}^{N}I_{\{ \hat{\mbox{\boldmath $\lambda$}}_{i}\leq \mbox{\boldmath $x$} \}}$, with $I_{A} = 1$ if $A$ is true and 0 otherwise, is the empirical distribution function and $P(\mbox{\boldmath $x$})$ is the true distribution function. Therefore, clustering for the estimate $\hat{\mbox{\boldmath $\lambda$}}_{i}$ obtained by our proposed method is equivalent to clustering for the true value $\mbox{\boldmath $\lambda$}_{i}$ ($i = 1, 2, \ldots, N$). Since $\phi(\Delta( \mbox{\boldmath $\lambda$},S ))$ is an increasing function of the angle deviation which can be used in defining a within cluster sum of angle deviations, the criterion considered here minimizes the within-cluster sum of angle deviations.

We assume that $\{\mbox{\boldmath $\lambda$}_{1},\ldots,\mbox{\boldmath $\lambda$}_{N}\}$ is a sample of independent observations on some probability measure $P$. Here, we consider the empirical measure

where $\mathcal{S}$ is a set of subspaces. For a fixed set of subspaces $\mathcal{S}$, we can obtain

where $\mbox{\boldmath $\lambda$}_{1},\ldots,\mbox{\boldmath $\lambda$}_{N}$ are the $N$ vectors of the sample points which can be clustered by minimizing the within cluster sum of angle deviations. Let $\mathcal{S}_{N}$ be the set of subspaces that minimizes $W(\cdot, {\hat P}_{N})$ (i.e., the set of optimal clustered subspaces based on the samples) and $\bar{\mathcal{S}}$ be the set of subspaces that minimizes $W(\cdot, P)$. Provided that $\bar{\mathcal{S}}$ can be uniquely determined, it is expected that $\mathcal{S}_{N}$ should lie close to $\bar{\mathcal{S}}$.

For a probability measure $Q$ on $\mathbb{R}^{r}$ and a finite set of subspaces $\mathcal{S}$ of $\mathbb{R}^{r}$, we define

and

For a given value of $k$, the set of optimal clustered subspaces based on the samples $\mathcal{S}_{N} = \mathcal{S}_{N}(k)$ is chosen to satisfy $\Phi(\mathcal{S}_{N}, {\hat P}_{N})=m_{k}({\hat P}_{N})$ and the set of optimal population clustered subspaces $\bar{\mathcal{S}}=\bar{\mathcal{S}}(k)$ is chosen to satisfy $\Phi(\bar{\mathcal{S}},P)=m_{k}(P)$.

To define the distance measures, we have the following assumption:\\ Assumption A. Suppose that $\int \phi(|| x ||)P(dx)<\infty$ and that $m_j(P)>m_k(P)$ for $j=1,\cdots,k-1$.

Our aim here is to prove a consistency result for the cluster subspaces that $\mathcal{S}_{N} \xrightarrow{a.s.} \bar{\mathcal{S}}$. To show $\mathcal{S}_{N}\xrightarrow{a.s.}\bar{\mathcal{S}}$, we first consider the subspace distance defined in LiWang2006:

Definition 1. The symmetric distance between any $m$ dimensional subspace $U$ and $n$-dimensional subspace $\tilde{U}$ is defined as

where $(\mbox{\boldmath $u$}_{1},\ldots,\mbox{\boldmath $u$}_{m})$ and $(\tilde{\mbox{\boldmath $u$}}_{1},\ldots,\tilde{\mbox{\boldmath $u$}}_{n})$ are the bases of $U$ and $\tilde{U}$ respectively. Note that this subspace distance satisfies the triangle inequality

where $W$ is any non-empty subspace.

Remark 2: The angle $\Delta(\cdot, \cdot)$ used in Eq. ((ref)), $$\Delta(U,\tilde{U})=\sqrt{\min(m,n)-\sum_{i=1}^{m}\sum_{j=1}^{n}(\tilde{\mbox{\boldmath $u$}}_{j}^{'}\mbox{\boldmath $u$}_{i})^{2}},$$ is different from the distance defined in Definition 1 for measuring the subspaces distance. In fact, the angle $\Delta(\cdot, \cdot)$ projects the low subspace onto the high-dimensional subspace, while the distance $D(\cdot,\cdot)$ projects the high-dimensional subspace onto the low subspace. In order to avoid different dimension subspaces being treated as the same subspace, $D(\cdot, \cdot)$ is used to characterize the difference between two subspaces. For example, consider a two-dimensional plane and a line which is parallel to this plane as two subspaces, if we use the distance $\Delta(\cdot, \cdot)$ to characterize the difference between these two subspaces, it is likely to get the result that these two subspaces are treated as the same subspace because the angel between these two subspaces is $0$. However, using the distance $D(\cdot, \cdot)$ can avoid this issue. Here, we further define a distance measure similar to the Hasudroff distance:\\ Definition 2. Let $\mathcal{X}$ and $\mathcal{Y}$ be two non-empty compact subsets that contain multiple subspaces. We define their Hausdorff distance $D_{H}(\mathcal{X}, \mathcal{Y})$ by

where $x \in \mathcal{X}$ is a subspace rather than a point, and the $D(\cdot,\cdot)$ is the distance between two subspaces defined in Definition 1.

By Definition 2, we have $D_{H}(\mathcal{X},\mathcal{Y}) < \delta$ if and only if every subspace of $\mathcal{X}$ is within the distance $\delta$ of at least one of the subspaces of $\mathcal{Y}$, and vice versa. Suppose $\mathcal{X}$ contains exactly $k$ distinct subspaces, and that $\delta$ is chosen to be a value less than half of the minimum distance between the subspaces of $\mathcal{X}$. Then, if $\mathcal{Y}$ is any set of $k$ or fewer subspaces for which $D_{H}(\mathcal{X},\mathcal{Y})< \delta$, the $\mathcal{Y}$ must contain exactly $k$ distinct subspaces. Therefore, the almost sure convergence of $\mathcal{X}_{N}$ in the above sense of distance could be translated into the almost sure convergence of subspaces with a suitable labeling. By definition, for any two subspaces $S_{1}$ and $S_{2}$, if $D(S_{1},S_{2}) < \delta$, then $\Delta( S_{1},S_{2})<\delta$. We will provide the following theorem for the consistency of clustering procedure. Since the conclusion of the theorem is in terms of almost sure convergence, there might be aberrant null sets of subspaces $S$'s for which the convergence does not hold. In order to estimate the parameters in the model presented in Eq. ((ref)) and to prove the consistency of the estimators, similar to Bonhomme2015 and Ando2016, we add the following assumption, Assumption B, that each group must have a certain proportion individuals. This assumption also guarantees the uniqueness because the null set situation is excluded.

Assumption B. All units are divided into a finite number of subspaces $k$, each of them containing $N_j$ units such that $0<\underline{a}<N_j/N<\bar{a}<1$.

For notational simplicity and clarity, we assume that $\lambda_i, i = 1, \cdots, N$ is known in the following theorem.

The proof of Theorem 1 is presented in the Supplementary Materials S2.

Consistency of the estimators

In this subsection, we discuss the asymptotic properties of the proposed estimators. Recall that the proposed estimators can be obtained by minimizing the objective function defined in Eq. ((ref)) subject to the constraints $g_{i}=\arg\min\limits_{j \in\{1,\ldots,k\}} || B_{j}^{T}\mbox{\boldmath $\lambda$}_{i}||$, $\mbox{\boldmath $F$}_{j}^{'}\mbox{\boldmath $F$}_{j}/T=I_{r}(j=1,\ldots,k),\mbox{\boldmath $\Lambda$}_{j}^{'}\mbox{\boldmath $\Lambda$}_{j}(j=1,\ldots,k)$ being diagonal. While the consistency of the subspace clustering procedure has been discussed in Section 4.1, we have the following theorems (Theorems 2, 3 and 4) to show the property of the estimators when $T$ and $N$ are large. To show the property of the estimators, Assumptions C, D, E and F, are needed as presented in Bai2009 and Ando2016 and we present Assumptions C--F in the Supplementary Materials S1.

The proofs of Theorems 2--4 are presented in the Supplementary Materials S3.

Monte Carlo Simulation Studies

In this section, Monte Carlo simulation studies with different settings are used to illustrate the proposed methodologies and to study the finite sample properties of the proposed methods. We assume that $\mbox{\boldmath $\lambda$}_{i}$ ($i = 1, 2, \ldots, N$) comes from different subspaces with dimension $d_{1}, d_{2}, \ldots d_{k}$ and $\mbox{\boldmath $\lambda$}_{i}$ can follow different probability distributions over different subspaces. Furthermore, we consider that $\mbox{\boldmath $\lambda$}_{i}, i=1,2,\ldots, N$, can have moderate noise. The simulation results are based on $100$ simulations.

Setting 1

We consider the situation that there are three different subspaces in the $\mathbb{R}^{3}$ space with known dimensions $d_{1}$, $d_{2}$ and $d_{3}$, i.e., $k = 3$ and $r = 3$. The bases of the three subspaces with dimensions $d_{1}$, $d_{2}$ and $d_{3}$ are represented as $\mbox{\boldmath $\alpha$}_{1}$, $\mbox{\boldmath $\alpha$}_{2}$ and $\mbox{\boldmath $\alpha$}_{3}$, respectively. Let $N_{\ell_1 \times \ell_2}(\mu, \sigma^{2})$ be a $\ell_1 \times \ell_2$ matrix whose elements are random variables that are independent and identically distributed as normal with mean $\mu$ and variance $\sigma^{2}$. Under this setting, we generate the panel data $y_{it}$ in the $j$-th subspace ($j = 1, 2, 3$; $i = 1, 2, \ldots, N_{j}$; $t = 1, 2, \ldots, T$) based on the panel data model in Eq. ((ref)) with $N_{1} = N_{2} = N_{3} = N = 100$, $T = 6$ and the following scheme:

In the simulation study, we compare the performance of the proposed LSSC method with the GFE method Bonhomme2015 and the estimation method proposed by Bai2009 (BAI) in terms of the biases and the root mean squared errors (RMSEs). The simulated biases and RMSEs of the estimators obtained from the GFE, BAI and LSSC estimation methods for Setting 1 are presented in Table (ref). From Table (ref), we observe that the proposed LSSC method has better performance than the GFE and BAI's methods in terms of biases and RMSEs. It is noteworthy that although the theoretical proofs of the asymptotic properties require $T$ to be large, the simulation results show that the proposed method performs well even when $T$ is small.

To compare the performance of the clustering methods, we present the simulated average misclassification rates of the clustering methods based on GFE and the proposed LSSC in Table (ref). From Table (ref), we can see that the simulated misclassified rates of the LSSC method are lower than the corresponding misclassified rates of the GFE method.

To verify the consistency of the proposed LSSC method, in addition to $N_{1} = N_{2} = N_{3} = N = 100$, we consider different sample sizes $N = 50$, 200, 300 and 500 in order to study the effect of the sample size on the biases and RMSEs. We also consider different values of time period $T = 5$, 10, 30, 50, 100 to study the effect of the time period on the biases and RMSEs of parameter $\mbox{\boldmath $\beta$}$. The results are presented in Tables (ref) and (ref). From Tables (ref) and (ref), we can see that the biases and RMSEs of the LSSC estimators decrease as the sample size $N$ or time period $T$ increases, which verifies Theorem 2.

Setting 2

The method proposed by Ando2016 is effective when the regressors are not correlated with factors and factor loadings under large $N$ and large $T$, but it does not performs well when the regressors and the factors are correlated, such as the set-up in Setting 1. To further compare the proposed method with the method proposed by Ando2016, we consider the following settings:

Model Selection and Possible Extensions

In the previous sections, we assume that the number of subspaces and the dimension of subspaces are known. In this section, we extend the procedure to a more general setting in which the dimension of factors is unknown and discuss the situations that the number of subspaces for the factors and the dimension of these subspaces are unknown.

Determine the number of subspaces for factors

One of the critical aspects of cluster analysis is to determine the number of subspaces empirically based on the observed data. For experimental data, however, there is no real “true" number of subspaces, but only a choice of the suitable value of $k$ which can provide stable and replicable results with a good fit to the data. In fact, the problem of estimating the subspaces number is a challenging model selection problem. Here, we are not intended to give a detailed review of all the existing methods for obtaining the number of subspaces for factor, but we aim to provide a feasible solution based on the work by Liu_Ma2013.

Liu_Ma2013 proposed a novel objective function named low-rank representation (LRR), which seeks the lowest rank representation among all the candidates that can represent the samples as linear combinations of the bases in a given dictionary. The computational procedure of LRR is to solve a nuclear norm regularized problem Fazel2002, which is a convex optimization problem that can be solved in polynomial time. The estimate of the number of subspaces can be obtained as Liu_Ma2013

where $\tau$ is a cut-off threshold, $\sigma_{i}$ denotes a singular value of the normalized Laplacian matrix of the affinity matrix of data, $\text{int}[a]$ is the nearest integer of a real number $a$ and $f_{\tau}$ is a summation function which counts different values regarding that $\sigma_{i}<\tau$ defined as $f_{\tau}(\sigma)=1$ if $\sigma \geq \tau$ and $f_{\tau}(\sigma)=\log_{2}(1+\frac{\sigma^{2}}{\tau^{2}})$ if $\sigma < \tau$, where $0<\tau<1$ is a parameter. Specifically, we can use the following steps to obtain the number of subspaces:

To verify the performance of the above algorithm in obtaining the number of subspaces, we assume that the numbers of subspaces in Settings 1 and 2 are unknown and apply the above algorithm to obtain the number of subspaces based on each simulated data set. Based on the simulation study, we find that the above algorithm can obtain ${\hat k} = 3$ correctly for all the 100 simulated data sets. For future research, evaluating the performance of the above algorithm for obtaining the number of subspaces under different settings (e.g., different sample sizes, different number of subspaces, etc.) is of interest.

Determine the dimension of factors and the dimension of subspaces

Determination of the dimension of factors (dimension of ambient space) is an interesting research topic. The dimension of factors can be specified based on the particular practical problem using professional or expert knowledge. When professional or expert knowledge about the dimension of factors is not available, Bai_Ng2019 developed a regularization criterion to determine the number of factors, and this criterion is more stable when the nominal number of factors is inflated by the presence of weak factors or large measurement noise. To choose the dimension of factors $r \in [0, rmax]$, the expression of the criterion is $$\bar{r} = \min\limits_{r = 0, \cdots, rmax} \log\left[ 1-\sum\limits_{j=1}^r(D_{jj}-\gamma)^2 \right] + kg(N,T),$$ where $g_{N,T}=\frac{N+T}{NT}\log(\frac{NT}{N+T})$, $D_{jj}$ is the $j$-th singular value of the scaled observable data, and $\gamma$ is a constant threshold. Through Monte Carlo simulation, we found that this criterion performs well in determining the number of factors when the factors and factor loadings have the subspace structure.

When the dimensions of the subspaces are unknown, determining the dimension for each subspace is still an open and challenging problem. In this section, we suggest to obtain the solution of the optimal model selection as

where $SSR(\hat{Z})$ represents the mean squared errors under the subspaces set $\hat{Z}$ (i.e., a measure of the data fidelity), $\tau$ is the error tolerance, $d_j$ is the dimension of the $j$-th subspace and $N_j$ is the number of individuals in $j$-th subspace, $k$ is the number of subspaces, $\hat{\sigma}^2$ is the estimated variance, and $\sum_{i=1}^{k}d^2_i \frac{T+N_i}{NT}\log(TN_i)$ is the penalty term which measures the model complexity under the subspaces set $\hat{Z}$. The proposed criterion can be viewed as a tradeoff between how well the model fits the data and the model complexity. It can be shown that the penalty function $d^2_j \frac{T+N_j}{NT}\log(TN_j)\to 0$ and $\min\{N,T\}d^2_j \frac{T+N_j}{NT}\log(TN_j)\to \infty$ as $T,N\to \infty$ and $T/N$ converges to constant.

To examine the proposed method for model selection, we simulated the panel data from the models with number of factors $r = 3$ and 4 and then obtain the solution of the optimal model selection in Eq. (ref) with different number of units in each subspace and different time period $T$. We consider that there is no covariate, i.e., $\beta = 0$ and we use the error tolerance $\tau=\frac{10 (d_1^3 + \cdots + d_{k-1}^3)}{\min\{N,T\}}$. The simulated percentages of identifying the correct dimension of subspaces (based on 1000 simulations for each setting) are presented in Tables (ref) and (ref) for $r = 3$ and 4 with three and four subspaces, respectively. From the simulation results in Tables (ref) and (ref), the proposed model selection method performs reasonably well in the case of hyperplane, i.e., these subspaces have the same dimensions. Compared with the case of the hyperplane, when the dimensions of the subspaces are not all the same, the simulated percentages of identifying the correct dimension can be lower to about 80%.

Real Data Application

In this section, we illustrate the proposed methodologies by using the real data provided by Bonhomme2015 and studying the linkage between income growth and democracy across different countries. Following Bonhomme2015, we use the linear dynamic model to identify the group membership and the linkage between income growth and democracy across countries, i.e.,

where $democracy_{it}$ is the democracy index (measured by the Freedom House indicator with values in between 0 (the lowest) and 1 (the highest)) for the $i$-th country at time $t$, $GDPpc_{it}$ is the GDP per capita of the $i$-th country at time period $t$, and $\mbox{\boldmath $\lambda$}_{g_i,i}$ and $\mbox{\boldmath $f$}_{g_i,t}$ are the unobservable grouped factor loadings and factors, respectively. Here, the dependent variable is the country's democracy index and the explanatory variables are the first-order lagged democracy index and the income of a country measured by the logarithm of GDP per capita.

The data set contains a balanced panel of 90 countries and 7 periods at a five-year interval over 1970--2000. First, using the information criteria suggested in Bai_Ng2019 to estimate the number of factors, we obtain the dimension of factor space as $r = 5$. Then, the number of subspaces is estimated as $k = 3$ based on Eq. ((ref)). The results are consistent with those presented in SuPhillips2016. Next, we use the criterion in Eq. ((ref)) to select the optimal model, and the results show that the optimal model have the dimensions $d_1 = d_2 = d_3 = 4$. Finally, we use BAI, GFE and LSSC methods to obtain the parameter estimates as $(\hat{\theta}_1, \hat{\theta}_2)$ and corresponding fitting errors (defined as $\hat{SSR}=\frac{1}{NT}\sum\limits_{i=1}^N\sum\limits_{t=1}^T(democracy_{it}-\theta_1 democracy_{i(t-1)} -\theta_2log GDPpc_{i(t-1)} - \mbox{\boldmath $\lambda$}_{g_i,i} \mbox{\boldmath $f$}_{g_i,t})^2$). The estimated results are presented in Table (ref). From Table (ref), we can see that all these estimates imply the effect of income on democracy is positive, but the LSSC method has the smallest fitting error.

In order to visualize the group membership obtained by the proposed method, we put these grouped countries on a world map in Figure (ref) in which the countries in the same group are represented in the same color. The detailed lists of grouped countries are presented as followings:

From these groupings, it can be seen that most of the early developed countries are distributed in the first group, which has a certain relationship with the economic and political structure. It includes the United States and Canada, most of the countries in continental Europe, coastal countries of South America and Australia. Most of the countries in the second group are developing countries with rapid economic development in Asia, Africa and South America, which includes China, Brazil and South Africa. Japan and South Korea also belong to the second group because they are both countries with high-speed economic development during this period and similar culture and policies. Most of the countries in the third group have slower development and relatively backward economies and policies during this period.

Concluding Remarks

In this paper, we consider a panel data model that allows the covariates and the unobservable latent variables to be correlated. We propose a subspace clustering method for factor loadings of the panel data model that captures the grouped unobserved heterogeneity. The common regression parameters, grouped unobservable factor structure and group membership can be estimated simultaneously with the proposed method. The asymptotic results show that the subspace clustering and the estimators are consistent. The Monte Carlo simulation results show that the proposed methodologies outperform the existing methods under different settings. Under the model considered in this paper, we propose a consistent model selection criterion to determine a suitable subspace dimension. We also discuss some possible future research directions in determining the number of subspaces and factor dimension when these values are unknown. These issues are under investigation and we hope to report the results in a future paper.