Blocked Clusterwise Regression

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

41,712 characters · 16 sections · 39 citation commands

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

abstractA recent literature in econometrics models unobserved cross-sectional heterogeneity in panel data by assigning each cross-sectional unit a one-dimensional, discrete latent type. Such models have been shown to allow estimation and inference by regression clustering methods. This paper is motivated by the finding that the clustered heterogeneity models studied in this literature can be badly misspecified, even when the panel has significant discrete cross-sectional structure. To address this issue, we generalize previous approaches to discrete unobserved heterogeneity by allowing each unit to have multiple, imperfectly-correlated latent variables that describe its response-type to different covariates. We give inference results for a k-means style estimator of our model and develop information criteria to jointly select the number clusters for each latent variable. Monte Carlo simulations confirm our theoretical results and give intuition about the finite-sample performance of estimation and model selection. We also contribute to the theory of clustering with an over-specified number of clusters and derive new convergence rates for this setting. Our results suggest that over-fitting can be severe in k-means style estimators when the number of clusters is over-specified.

Introduction

We often believe that there may be significant cross-sectional heterogeneity in the structural relationship between observed covariates and response. Even with panel data, however, estimating distinct regression coefficients for each cross-sectional unit can be noisy or infeasible when the time dimension is small. Clusterwise regression methods (e.g. LinNg2012, BonhommeManresa2015), which model individual heterogeneity as a function of a one-dimensional discrete latent type, have recently become popular as a viable compromise between the common parameter assumption and full heterogeneity. However, as we show, this discretization of heterogeneity can be badly misspecified even when the panel has significant cross-sectional structure. By introducing multiple, imperfectly correlated latent types, we can relieve this issue and significantly enrich the set of panel structures that can be handled by clustering methods. In particular, our approach is motivated by a class of data-generating processes where units are clustered along multiple latent dimensions or “response-types” to distinct blocks of the covariate vector. We motivate this generalization with several examples from finance and production function estimation. The main contribution of this paper is to modify existing clustering methods for use in this larger family of models and show that they can likewise be used to perform inference on regression parameters. \\

Following BonhommeManresa2015 (henceforth BM), we establish consistency and asymptotic normality for a k-means style estimator in our setting. The estimation algorithm is iterative and alternates between (1) solving a least-squares problem to estimate cluster parameters for each block and (2) updating latent types for each cross-sectional unit based on a unit-wise predictive criterion. As in BM, our proof proceeds by establishing asymptotic equivalence with the oracle estimator where each unit's latent types are known. We extend the approach in AndoBai2016 to give a $C_p$ style information criterion to choose the number of clusters (types) for all the latent variables simultaneously. \\

In general, the true number of clusters in a given data set is unknown. Thus, the behavior of estimators with a misspecified number of clusters is important both for model selection theory as well as for our understanding the finite-sample properties of clustering estimators. Here, we make some contributions to the theory of models with an over-specified number of clusters, improving the convergence rates given in LiuOverspecifiedGroups for the linear regression setting. In contrast to the well-specified case, difficulty obtaining the “fast rate” $O_p(\frac{1}{NT})$ when we over-specify the number of clusters suggests that over-fitting may be severe when the number of clusters is over-specified. We conjecture that $\sqrt{T}$-consistency may be optimal for over-specified models.

Motivating Example - Production Function Estimation

Consider panel data on firms' production levels and factor usage. We are interested in estimating the firm-specific production functions

equation[equation omitted — 161 chars of source]

where $y_{it}$ is a measure of output and $L_{it}, K_{it}, M_{it}$ are labor, capital and materials (all in logs), and $\text{Elec}_{it}$ is a measure of electricity usage. Suppose that the heterogeneity in factor elasticities can be well approximated by \[ \theta_{i\ell} \in \{\theta^{low}_{\ell}, \theta^{mid}_{\ell}, \theta^{high}_{\ell} \} \quad 1 \leq \ell \leq 4 \]

Ignoring endogeneity in input choice, we consider estimation of (ref) with clusterwise regression. The problem with this approach is readily apparent - although each firm can only have one of $3 \cdot 4 = 12$ elasticity types, $\theta_i$ can take up to $3^4 = 81$ distinct values. Thus, estimating this model with clusterwise regression would require $k=81$ clusters to be well-specified. For a panel of $200$ firms, this would lead to estimation with approximately $N/81 \leq 3$ firms in each regression, in spite of significant cross-sectional homogeneity. However, with $k=81$ clusters the model is also significantly over-parameterized. For instance, there will be $27$ distinct clusters with each level of labor elasticity coefficient. \\

The problem is that current clustering models assume limited heterogeneity in the individual parameter vectors $\theta_i$. In our example, however, cross-sectional heterogeneity takes the form of a few discrete elasticity levels for each input factor, while the support of $\theta_i$ itself is large. This suggests a model with multiple latent heterogeneity types. For instance \[ \theta_i = (\theta_1(c_{i1}), \theta_2(c_{i2}), \theta_3(c_{i3}), \theta_4(c_{i4})) \] with latent type $c_{i\ell}$ for $1 \leq \ell \leq 4$ controlling the elasticity level of factor $\ell$.

Related Literature and Outline

Early contributions to the econometric literature on clustering include Sun2005 and BuchinskyHotz2005. Linear panel data models with discrete unobservable heterogeneity have recently been studied in LinNg2012, BonhommeManresa2015, Su2016, Su2016Homogeneity, DzemskiOkui2018. Our asymptotic normality results for the well-specified case closely follow the analysis pioneered in BonhommeManresa2015. AndoBai2016 extends clustering methods to linear factor models and gives an information criterion for choosing the number of clusters. We develop a similar $C_p$-style criterion in our setting. Outside of the linear case, Zhang2019 and Chen2019 study clustered linear conditional quantile regression, and BonhommeManresa2019Discretizing considers discrete latent types as an approximation to continuous unobserved heterogeneity. LiuOverspecifiedGroups studies clustering in M-estimation with an over-specified number of groups. We build on their techniques and significantly sharpen their rate results for the linear case. In contemporaneous work, Cheng2019 consider a clustering model with two latent types in a GMM setting. By contrast, we allow for $B > 1$ latent types in a linear model with individual fixed effects. Clusterwise regression was initially proposed in Spath1979 as “Algorithm 39 - Clusterwise Linear Regression.”\\

Further afield, this paper is related to a number of literatures in statistics and computer science, such as the literature on clustering functional data, e.g. Wasserman2005CATS, Yamamoto2014, LintonVogt2019, and subspace clustering, e.g. Candes2012. In statistics, related methods include homogeneity pursuit, proposed in Ke2015. See also Ke2016 and Lian2019. In the Bayesian literature, clusterwise regression is also known as multilevel regression; see GelmanHill2007. \\

The paper is organized as follows - we introduce our model and estimator in section (ref). Asymptotic properties of the estimator and consistency of model selection are given in section (ref). Section (ref) discusses models with an over-specified number of clusters. Monte Carlo simulations are given in section (ref), and proofs in section (ref). Supplementary appendix (ref) collects technical lemmas and other ancillary discussions.

Model and Estimation

Model

Let $y_{it}$ and $x_{it}$ denote repsonse and covariates for $t = 1, \dots, T$ time periods and $i = 1, \dots, N$ cross-sectional observations. The covariate vector $x_{it} \in \mathbb{R}^p$ is divided into $1 \leq \ell \leq B$ blocks $x_{it}^{\ell}$, where $x_{it} = (x_{it}^1, \dots, x_{it}^{B})$, and $B$ denotes the total number of blocks. We let $k = (k_1, \dots, k_{B})$, where $k_{\ell}$ denotes the number of distinct latent types (clusters) associated with the $\ell^{th}$ block. Possible cluster assignments are denoted $c = (c_1, \dots, c_{B}) \in \prod_{\ell} [k_{\ell}] \equiv \mathcal{C}$. For instance, a unit in cluster $1$ in the first block and cluster $3$ in the second block would have $c=(1, 3)$. \\

Each cross-sectional unit belongs to exactly one cluster for each block. We let $\gamma: [N] \to \prod_{\ell} [k_{\ell}]$ denote an assignment of cross-sectional units to cluster vectors, so that $\gamma(i) = c_i$. The set of all possible cluster assignments is denoted $\Gamma$. In the main specification, we assume that the response $y_{it}$ is given by

equation[equation omitted — 67 chars of source]

with the $\ell^{th}$ block parameter selected by the latent variable $c_{i\ell}$

equation[equation omitted — 111 chars of source]

Thus, each covariate grouping $\ell$ is associated with $k_{\ell}$ parameter sub-vectors with $\theta_{\ell}(c_{\ell}) \in \mathbb{R}^{d_{\ell}}$. Our goal is to jointly estimate the true parameter $\theta^0$ and true cluster assignments $\gamma^0(i) = (c^0_{i1}, \dots, c^0_{iB})$ of each cross-sectional unit. In appendix (ref), we also give results for the model with individual fixed effects

equation[equation omitted — 61 chars of source]

\quad \\ Relationship with Clusterwise Regression: The model above nests clusterwise regression (as in LinNg2012) when $B=1$. For $B > 1$, it is statistically equivalent to clusterwise regression when the conditional pdf $\mathbb{P}(c^0_{i(-\ell)} | c^0_{i\ell})$ is degenerate (perfectly correlated types). Similarly, clusterwise regression ($B=1$) with exponentially many clusters $k = \prod_{\ell=1}^{B} k_{\ell}$ and exponentially many constraints nests our model. For instance, let $p=B$ and assume $\theta_{\ell} \in \{\pm 1\}$ for each $\ell$. Then our model would be equivalent to a clusterwise regression model with $2^p$ clusters and $p2^{p-1}$ linear equality constraints. \\

Example - Exchange Rates: Consider financial market data with $y_{it}$ the exchange rate against USD of country $i$ at time $t$, $p^{oil}_{t}$ the price of crude oil, $b_{it}$ a measure of the country's business cycle and $r_{t}$ the US discount rate, we model \[ y_{it} = \theta_{i1} p^{oil}_{t} + \theta_{i2} b_{it} + \theta_{i3} r_{t} + e_{it} \] Due to differences in national industry composition, the magnitude and composition of foreign trade, financial openness and so on, we may expect heterogenous marginal responses $\theta_{i \ell}$ of $y_{it}$ to each of the factors above. As in the introduction, we may model $\theta_{i\ell} \in \{\theta_{\ell}(1), \dots, \theta_{\ell}(k_{\ell})\}$, corresponding to $k_{\ell}$ different sensitivity levels to factor $\ell$. However, we don't expect these unobserved types to be perfectly correlated across factors. For instance, we might expect both Venezuela and China to have large $\theta_{i1}$ but very different $\theta_{i3}$. A factor error structure could be accommodated using the techniques in AndoBai2016.

Estimator

We define our estimator of the parameter $\theta^0$ and cluster assignment $\gamma^0$ as

equation[equation omitted — 204 chars of source]

We let $\widehat{Q}(\theta, \gamma)$ denote the sample risk in (ref). There are many algorithms available for the least squares partitioning problem above\footnote{See, for instance, the discussion in BM Appendix S1.}. One benchmark approach, known as Lloyd's algorithm (Lloyd1982) in the setting of k-means clustering, takes a coordinate ascent approach to the problem in (ref), alternately updating the parameters $\theta$ and assignments $\gamma$ until convergence. \\

Lloyd's Algorithm - Fix a division of the covariate vector $x_{it} = (x_{it}^1, \dots, x_{it}^{B})$ into blocks with $x_{it}^{\ell} \in \mathbb{R}^{d_{\ell}}$ and fix the number of clusters $k = (k_1, \dots, k_{B})$ in each block. Our approach is a modification of Lloyd's algorithm for k-means clustering. We perform coordinate ascent on the sample objective $\widehat{Q}(\theta, \gamma)$ by alternating parameter updates and cluster assignment updates until convergence. \\

(1) Randomly initialize parameters $\theta^1$ and cluster assignments $\gamma^1$. \\ (2) Given $\theta^s$, set $\gamma^{s + 1} = \operatorname*{argmin}_{\gamma \in \Gamma} \widehat{Q}(\theta^s, \gamma)$. \\ (3) Given $\gamma^{s + 1}$, update $\theta^s \to \theta^{s + 1}$. \\ (4) Repeat (2) and (3) until convergence. \\

Since problem (ref) is not generally convex, we repeat steps (1) through (4) from different random initializations (in parallel), and take the estimate that achieves the lowest sample risk $\widehat{Q}$. See appendix (ref) for more discussion of the implementation of this algorithm and related computational issues.

Asymptotic Properties

In this section, we investigate the asymptotic properties of the estimator $(\widehat{\theta}, \widehat{\gamma})$ defined above as $N, T \to \infty$. In what follows, we assume the data is generated from the model (ref) with $(\theta^0, \gamma^0)$ the true slope parameters and cluster assignment function. We let $\| \cdot \|$ denote the usual Euclidean norm.

Consistency

assumptionWe make the following assumptions \begin{enumerate}[label={(\alph*)}, ref ={\ref*{assumptions:consistency}.(\alph*)}, itemindent=.5pt, itemsep=.5pt] • For each $\ell \in [B]$, the parameter space $\Theta_{\ell} \subset \mathbb{R}^{d_{\ell} \times k_{\ell}}$ is compact • $\frac{1}{NT^2} \sum_{i=1}^N \sum_{t=1}^T \sum_{s=1}^T e_{it} e_{is} x_{it}' x_{is} \overset{p}{\to} 0$ • $\|\theta^0_{\ell a} - \theta^0_{\ell b} \| \equiv d(\ell, a, b) > 0$ for each pair of clusters $a, b \in [k_{\ell}]$ • Define $M(c, c', \gamma) \equiv \frac{1}{NT} \sum_{i, t} x_{it} x_{it}' \mathds{1}(c_i^0 = c')\mathds{1}(c_i = c)$ and let $\rho(c, c', \gamma) \equiv \lambda_{min}(M(c, c', \gamma))$. Then there exists $\delta > 0$ such that $\inf_{c', \gamma} \max_{c} \rho(c, c', \gamma) \geq \delta - o_P(1)$ as $N, T \to \infty$. \end{enumerate}

Assumption (ref) is the usual parameter space compactness condition. Assumption (ref) can be seen as limiting the time-series dependence of errors and covariates, averaged over cross-sectional units. Condition (ref) ensures that the clusters within each grouping are non-identical. The final assumption (ref) is the analogue in our setting of assumption S2(a). in BM. This condition is used to ensure curvature of the sample risk function $\widehat{Q}$. If there is a common parameter ($B = 1$, $k = 1$), this is the usual non-collinearity condition for pooled panel regression. See section (ref) in the appendix for further discussion, as well as section S4.2 in the supplementary appendix of BM.\\

Cluster Label Ambiguity. The minimizer $\operatorname*{argmin}_{\gamma \in \Gamma, \theta \in \Theta} \widehat{Q}(\theta, \gamma)$ is only unique up to permutations of the labels in $\mathcal{C}$ and their associated parameter vectors in $\theta$. Thus, the $c \in \mathcal{C}$ used to label estimated clusters in each block is arbitrary, and to resolve this ambiguity we need to fix a correspondence $\sigma_{\ell}: [k_{\ell}] \to [k_{\ell}]$ between true and estimated cluster parameters for each $\ell$\footnote{ Note that, for finite $T$, it can be the case that $\widehat{\theta}(\widehat{c}_i) \not = \widehat{\theta}(\widehat{c}_j)$, but $c_i^0 = c_j^0$, so the estimates $\widehat{\theta}(\widehat{c}_i)$ do not in general induce a well-defined estimator of any fixed cluster parameter $\theta^0_{\ell a}$.}. Let

equation[equation omitted — 170 chars of source]

and define the estimator of the true parameter $\theta^0_{\ell a}$ to be $\widehat{\theta}_{\ell \sigma({a})}$. Note that this is infeasible without access to the true parameters $\theta^0$.

lemmaUnder the assumptions in (ref), $\mathbb{P}(\sigma_{\ell} \, \text{invertible}) \to 1$ as $N, T \to \infty$

By the lemma, we can relabel the estimated clusters $\widehat{\theta}_{\ell \sigma(a)} \to \widehat{\theta}_{\ell a}$, and this is well-defined w.h.p as $N, T \to \infty$.

thmUnder the assumptions in (ref), for all groupings $\ell$ and $a \in [k_{\ell}]$, we have \[ \|\theta^0_{\ell a} - \widehat{\theta}_{\ell a} \| = o_P(1) \] equivalently \[ \min_{x \in \mathcal{C}} \|\widehat{\theta}(x) - \theta(c)\| = o_p(1) \quad \forall c \in \mathcal{C} \] as $N, T \to \infty$.

See section (ref) for the proof of the theorem and lemma.

Asymptotic Equivalence

In this section, we establish asymptotic equivalence of $\widehat{\theta}$ to the infeasible oracle estimator with known clusters. We need the following assumptions in addition to those already stated in (ref).

assumptionMake the following assumptions \begin{enumerate}[label={(\alph*)}, ref ={\ref*{assumptions:inference}.(\alph*)}, itemindent=.5pt, itemsep=.5pt] • $\frac{1}{NT^2} \sum_{i} \sum_{t, s} \|x_{it} \|^2 \|x_{is}\|^2 = O_P(1)$ • Define $M_{NT}^c = \frac{1}{NT} \sum_{i, t} \mathds{1}(c_i^0 = c) x_{it}x_{it}'$. Then there exists $\underline{\rho} > 0$ such that for all $a > 0$, this sequence of matrices satisfies $\min_{c \in \mathcal{C}} \lambda_{min}(M_{NT}^c) \overset{p}{\to} \underline{\rho}$ as $N, T \to \infty$. • There exist constants $b > 0$ and $d_1 > 0$ and sequence $\alpha(t) \leq e^{-bt^{d_1}}$ such that for all $i \in [N]$ $\{x_{it}\}_t$ and $\{x_{it} e_{it}\}_t$ are strongly mixing with coefficients $\alpha(t)$. • There exist constants $f>0$ and $d_2>0$ such that for all $i \in [N]$ and all $z > 0$, for all components $x_{it}^j$, $x_{it}^{j'}$ of the vector $x_{it}$ we have $\mathbb{P}(|x_{it}^j x_{it}^{j'} - E(x_{it}^j x_{it}^{j'}) | > z)$ and $\mathbb{P}( |e_{it} x_{it}^j - E e_{it} x_{it}^j | > z)$ are bounded above by $e^{1 - (z/f)^{d_2}}$. • The uniform limits $\max_{i \in [N]} \frac{1}{T} \sum_t E[e_{it} x_{it}] \to 0$ and $\min_{i \in [N]} \frac{1}{T} \sum_t \mathbb{E}(x_{it}'(\theta(c) - \theta(c')))^2 \to d(c, c')$ hold as $T \to \infty$, and $d(c, c') \geq d_{min} > 0$ for $c \not = c'$. • There exists $M' > 0$ such that for all $a > 0$ \[ \max_{i \in [N]} \mathbb{P} \left (\frac{1}{T} \sum_t \|x_{it} \|^2 > M' \right ) = o(T^{-a}) \] \end{enumerate}

We will show that $\widehat{\theta}$ is asymptotically equivalent to the infeasible oracle estimator where true cluster membership $c_i^0$ is known for all $i$. Define the problem

align[align omitted — 233 chars of source]

The following theorem shows that $\widehat{\theta}$ and $\Tilde{\theta}$ are asymptotically equivalent.

thmLet the assumptions in (ref) and (ref) hold. Then for any $a > 0$ and as $N, T \to \infty$, we have \begin{equation} \widehat{\theta} = \Tilde{\theta} + o_P(T^{-a}) \end{equation} Moreover, individual cluster estimates satisfy \begin{equation} \mathbb{P} \left (\exists i \in [N] \, s.t. \, \widehat{c}_i \not = c_i^0 \right ) = o(1) + o(NT^{-a}) \end{equation}

See appendix (ref) for the proof. Because of this theorem, for asymptotic sequences with $N$ growing at a sub-polynomial rate relative to $T$, it suffices to characterize the asymptotic distribution of the estimator $\Tilde{\theta}$.

Inference

Notation - To aid the exposition, we start with a few definitions. For $A \in \mathbb{R}^{p \times q}$, let $\operatorname{vec}(A) \equiv ((A^1)', \dots, (A^q)')' \in \mathbb{R}^{pq}$. Thinking of $\theta = \{\theta_1, \dots, \theta_{B}\}$ as a collection of matrices $\theta_{\ell} \in \mathbb{R}^{d_{\ell} \times k_{\ell}}$, we denote $\operatorname{vec}(\theta) = (\operatorname{vec}(\theta_1)', \dots, \operatorname{vec}(\theta_{B})')' \in \mathbb{R}^{d_{\theta}}$, where $d_{\theta} \equiv \sum_{\ell} k_{\ell} d_{\ell}$ is the total dimension of $\operatorname{vec}(\theta)$. For $1 \leq \ell \leq B$ and $a \in [k_{\ell}]$, we use the block index convention that $\operatorname{vec}(\theta)_{\ell a}$ refers to the $d_{\ell}$ dimensional sub-vector in the $a^{th}$ position of the $\ell^{th}$ block. Using the notation above, for $1 \leq \ell, s \leq B$ and $a \in [k_{\ell}]$, $b \in [k_s]$ define $\widehat{M} \in \mathbb{R}^{d_{\theta} \times d_{\theta}}$ and $v \in \mathbb{R}^{d_{\theta}}$ by

align[align omitted — 388 chars of source]

propThe solution $\Tilde{\theta}$ to problem (ref) satisfies \begin{equation} \widehat{M} \operatorname{vec}(\Tilde{\theta}) = v \end{equation}

The proof follows by taking the first order conditions of (ref) and rearranging. Note that the first order conditions $\nabla_{\theta_{\ell a}} \Tilde{Q}(\Tilde{\theta}) = 0$ can potentially vary with all the other parameters $\theta_{sb}$ in the model (for $s \not = \ell$). Therefore, in contrast to the $B = 1$ case considered in the existing literature, the estimator $\Tilde{\theta}$ is not equivalent to simply running $k$ separate regressions over the partition of the cross-sectional units. \\

Consider the following assumptions that allow us to characterize the asymptotic distribution of the infeasible $\Tilde{\theta}$.

assumptionWe make the following assumptions \begin{enumerate}[label={(\alph*)}, ref ={\ref*{assumptions:clt}.(\alph*)}, itemindent=.5pt, itemsep=.5pt] • There is a matrix $\Omega \succ 0$ such that for all $\ell, s \in [B]$ and $1 \leq a \leq k_{\ell}$, $1 \leq b \leq k_s$ \begin{equation} \frac{1}{NT} \sum_{i, j = 1}^N \sum_{t, t' = 1}^T E[e_{it} e_{jt'} \mathds{1}(c_{i \ell}^0 = a) \mathds{1}(c_{js}^0 = b) x_{it\ell} x_{jt's}'] \to \Omega_{\ell a, s b} \end{equation} as $N, T \to \infty$. • $\mathbb{E}[e_{it} \mathds{1}(c_{i \ell}^0 = a) x_{it \ell}] = 0$ for all $\ell, a$. • $\widehat{M} \overset{p}{\to} M$ as $N, T \to \infty$, with $M \succ 0$. • $\sqrt{NT} w \overset{d}{\to} \mathcal{N}(0, \Omega)$ as $N, T \to \infty$. \end{enumerate}

thmSuppose that the assumptions in (ref) are satisfied. Also suppose there is some $r > 0$ such that $\sqrt{N}T^{-r} = o(1) $ as $N, T \to \infty$. Then we have \begin{equation} \sqrt{NT}(\operatorname{vec}(\widehat{\theta} - \theta^0)) \overset{d}{\to} \mathcal{N}(0, M ^{-1} \Omega M) \end{equation}

The proof of this theorem is given in appendix (ref). \\

Consider the case where cross-sectional units are independent, then under assumption (ref), the terms in (ref) with $i \not = j$ vanish. In this case, we propose the HAC estimators

align[align omitted — 337 chars of source]

where $\widehat{M}$ is as in equation (ref). Variance estimators of this form were originally proposed in Arellano1987, and their asymptotic theory for $N, T \to \infty$ jointly was first analyzed in Hansen2007Robust. For further discussion on adapting the results of Hansen2007Robust to our setting, see appendix (ref).

Model Selection

In this section we let $k^0 = (k_1^0, \dots, k_{B}^0)$ denote the true number of clusters and develop a Cp-like criterion to estimate $k^0$. We suppose that prior information can be used to bound the true number of clusters from above $k^0 \leq k_{max}$. So far, we have defined the sample risk $\widehat{Q}(\theta, \gamma)$ with domain the true parameter space i.e. $\theta \in \prod_{\ell} \mathbb{R}^{d_{\ell} \times k^0_{\ell}}$ and $\gamma: [N] \to \prod_{\ell} [k^0_{\ell}]$. However, note that $\widehat{Q} = \frac{1}{NT} \sum_{i, t} (y_{it} - x_{it}'\theta(\gamma(i)))^2$ only varies through $\theta(\gamma(i)) \in \mathbb{R}^p$. Thus, we can extend the domain of $\widehat{Q}$ to models with $k \not = k^0$, since $\theta(\gamma(i)) \in \mathbb{R}^p$ for any conformable $(\theta, \gamma)$.\footnote{Formally, let $\widehat{Q}: \bigcup_{k \geq 0} \left \{\prod_{\ell} \mathbb{R}^{d_{\ell} \times k_{\ell}} \times [\, [N] \to \prod_{\ell} [k^0_{\ell}]\, ]\right \} \to \mathbb{R}_{\geq 0}$ with $\widehat{Q}(\theta, \gamma) = \frac{1}{NT} \sum_{i, t} (y_{it} - x_{it}'\theta(\gamma(i)))^2$}\\

We slightly strengthen some assumptions above

assumptionImpose the following assumptions \begin{enumerate}[label={(\alph*)}, ref ={\ref*{assumptions:mo}.(\alph*)}, itemindent=.5pt, itemsep=.5pt] • For all $c \in \mathcal{C}^{k_0}$, $\frac{1}{NT} \sum_{i, t} e_{it} x_{it} \mathds{1}(c_i^0 = c) = O_p(1/\sqrt{NT})$ • With $\rho(c, c', \gamma)$ defined as in assumption (ref), there exists $\delta > 0$ such that for all $k \geq k^0$ we have \[ \rho^k_{NT} \equiv \min_{c' \in \mathcal{C}^{k_0}} \min_{\gamma \in \Gamma^k} \max_{c \in \mathcal{C}^k} \rho(c, c', \gamma) \geq \delta - o_p(1) \] • As $N, T \to \infty$ \[ \inf_{j \in [N]} \lambda_{min} \left (\frac{1}{T} \sum_t E[x_{jt} x_{jt}'] \right ) \to \underline{\lambda} > 0 \] • For some $0 < \epsilon < \frac{1}{2} \wedge \frac{d_1 d_2}{d_1 + d_2}$ we have $\log N = o(T^{\epsilon})$ as $N, T \to \infty$ , where $d_1, d_2$ are the mixing and tail parameters defined in assumptions (ref) and (ref) \end{enumerate}

We can think of assumption (ref) as stating that a CLT holds for $e_{it} x_{it} \mathds{1}(c_i^0 = c)$ for each $c \in \mathcal{C}^{k_0}$. This will be easiest to satisfy when $E[e_{it} x_{it} \mathds{1}(c_i^0 = c)] = 0$, a stronger form of unconfoundedness. Assumption (ref) is the extension of assumption (ref) in our consistency proof to the case of models with misspecified number of clusters. See section (ref) of the supplementary appendix for a discussion of this condition. In the stationary case with identically distributed cross-sectional units, having $E[x_{it} x_{it}']$ of full rank is sufficient for assumption (ref). Finally, assumption (ref) requires that $\log N$ is sub-polynomial in $T$ as $N, T \to \infty$. \\

Information Criterion - Let $(\widehat{\theta}^k, \widehat{\gamma}^k)$ be the minimizer of $\widehat{Q}$ with $(k_i)_i$ clusters and denote $\widehat{Q}(k) = \widehat{Q}(\widehat{\theta}^k, \widehat{\gamma}^k)$. Then we define the $C_p$ criterion

align[align omitted — 165 chars of source]

We have the following result on consistency of model selection

thmSuppose the assumptions in (ref) hold. Let $f(N, T)$ be such that $f(N, T) \to 0$ and for some $\epsilon$ as in assumption (ref), $f(N, T) T^{1-3\epsilon} \to \infty$ as $N, T \to \infty$. Then \[ \mathbb{P}(\widehat{k} = k^0) \to 1 \] as $N, T \to \infty$

For the proof, see appendix (ref).

remark[Choice of $f$] While any function $f(N, T)$ satisfying the conditions of the theorem will give asymptotically consistent model selection, the choice of $f$ will significantly affects finite sample performance. To put $f(N, T)$ on the same scale as $\widehat{Q}(k)$, we use $f(N, T) = \widehat{\sigma}^2 g(N, T)$ in our simulations, where $\widehat{\sigma}^2$ is a consistent estimate of the long run variance $\lim_{N, T} \frac{1}{NT} \sum_{i, t} E[e_{it}^2]$. By lemma (ref) in the appendix, $\widehat{\sigma}^2 \equiv \widehat{Q}(k_{max})$ is such a consistent estimator. We find good performance with $g(N, T) = \frac{\log T}{T}$ in our simulations. Alternatively, e.g. $g(N, T) = \frac{\log T}{T^{1-\epsilon'}}$ for small $\epsilon'$ can be used to be technically consistent with the theory.

Overspecification of $k$

In this section, we report new results on the performance of k-means style estimators with an over-specified number of clusters. The work in this section builds on and sharpens the results in LiuOverspecifiedGroups for the case of linear regression. Our proof of model selection consistency in Theorem (ref) heavily relies on the following result.

thmSuppose the assumptions in (ref) hold. Then \begin{align} \sup_{i \in [N]} \|\widehat{\theta}^k(\widehat{c}^k_i) - \theta^0(c_i^0)\|^2 &= o_p(T^{-1 + 4 \epsilon}) \\ \min_{x \in \mathcal{C}^k} \|\widehat{\theta}^k(x) - \theta^0(c)\|^2 &= o_p(T^{-1 + 3 \epsilon}) \\ \frac{1}{N} \sum_i (\widehat{\theta}^k(\widehat{c}_i^k) - \theta^0(c_i^0))^2 &= o_p(T^{-1 + 3 \epsilon}) \end{align}

proofFollows from Proposition (ref), Proposition (ref), and Corollary (ref) in the appendix.

remarkThe preceding result can be compared with LiuOverspecifiedGroups Theorem 1 and Lemma 5.16, which give $o_p(1)$, $o_p(1)$, and $o_p(T^{\frac{-1}{2(1 + d)}})$ rates (respectively) for each of the losses above, with $d = \frac{d_1 d_2}{d_1 + d_2}$. Note that their setting is a more general model of clustered M-estimation. Our rate improvements come from (1) optimizing the Fuk-Nagaev inequality in Rio2011 for our purposes (see lemma (ref)) and (2) an inductive strategy that allows us to “boost” $O_p(T^{-1/4})$ rates arbitrarily close to $O_p(T^{-1/2})$. For a description of this approach, see lemma (ref) as well as the propositions and corollary referenced above.

remarkThe rate established in equation (ref) above is used to bound the magnitude of over-fitting for estimators with $k > k^0$, as in the second part of corollary (ref). A result of this form is necessary to determine the complexity penalty in (ref). In particular, in contrast to the result in LiuOverspecifiedGroups, theorem (ref) gives feasible rates for $f(T)$ that do not depend on mixing parameters and tail bounds of $e_{it}$ and $x_{it}$, which may be difficult or impossible to estimate.

remarkDifficulty obtaining the fast rate $\widehat{Q}(k) - \widehat{Q}^0 = O_p(\frac{1}{NT})$ for $k > k^0$ suggests that over-fitting may be severe under over-specification of $k$. Difficulty obtaining $\sqrt{NT}$-consistency of $\widehat{\theta}^k$ when $k > k^0$ suggests a type of incidental parameter problem. In fact, in the linear case it is known\footnote{See the example in BM, appendix S3.1.} that under $N \to \infty$, finite $T$ asymptotics, estimators with $k > k^0$ can suffer a bias of order up to $\frac{1}{\sqrt{T}}$.

Monte Carlo Simulations

In this section, we describe the results of our Monte Carlo simulations. All tables are reported in section (ref) of the appendix. Throughout, we denote

enumerate[itemsep=0.4pt] • $\text{Param. MSE} = \frac{1}{N} \sum_{i=1}^N \|\widehat{\theta}(\widehat{c}_i) - \theta^0(c_i^0)\|^2$ • $\text{Function MSE} = \frac{1}{NT} \sum_{i=1}^N \sum_{t=1}^T \|\widehat{\theta}(\widehat{c}_i)'x_{it} - \theta^0(c_i^0)'x_{it}\|^2$ • $\text{Cluster Loss} = \frac{1}{N} \sum_{i=1}^N \mathds{1}(\widehat{c}_i \not = c_i^0)$

We use two specifications for the joint distribution $(x_{it}, y_{it})$. (1) Specifications labeled AR(1) take $e_{it} \sim \text{AR}(1)$ and $x_{it} \sim \text{VAR}(1)$. The AR process $e_{it}$ has normal innovations, and $x_{it}$ has multivariate normal innovations with constant, diagonal covariance matrix. The respective autocorrelation parameters are $\rho_e = 0.3$, $\rho_x = 0.5$. (2) Specifications labeled HK use a heteroskedastic design inspired by Hansen2007Robust. With $x_{it}$ as above, we use $e_{it} = \rho e_{it-1} + v_{it} \cdot \sqrt{\frac{1}{2} + \frac{\|x_{it}\|^2}{2p}}$, with independent normal innovations $v_{it}$. Innovation variances are normalized so that $\mathrm{Var}(e_{it}) =1$ for all designs and all $(i, t)$. All simulations use $500$ independent samples. See appendix (ref) for additional details on the computational specification.

Estimator Performance

design[Cluster Separation] We let $p=4$, $k = (2, 2)$, $B=2$ and parameters \[ \theta_1^0 = \begin{pmatrix} 1 & \cos\alpha \\ 0 & \sin\alpha \end{pmatrix} \quad \theta_2^0 = \begin{pmatrix} 0 & -\sin\alpha \\ 1 & \cos\alpha \end{pmatrix} \] where the columns of $\theta_{\ell}$ list parameters of block $1 \leq \ell \leq B$. Thus, as $\alpha \to 0$, the cluster parameters in each block rotate towards each other. Cluster estimation accuracy radically worsens for small $\alpha$. Coverage is around $80-90\%$ for well-separated clusters. As $\alpha \to 0$, our confidence intervals do not account for variation due to cluster estimation, and coverage is poor. Parameter loss is inverse U-shaped in cluster separation $\alpha$. For small $\alpha$, classification of $c_i^0$ becomes worse, giving large losses on some units. For $\alpha$ near $0$, misclassification contributes less to parameter loss since the cluster centers are very close. Results are shown in Table (ref).

design[Sample Size $(N, T)$] We use the specification in the simulation above with $\alpha=\frac{\pi}{2}$. Cluster loss is quite insensitive to $N$, in line with the theory. Increasing $T$ has a much larger effect than $N$ on parameter loss and coverage. For $T=5$, we find coverage actually decreases with larger $N$, which could be an example of the over-fitting issue discussed in section (ref). Results are shown in Table (ref).

design[Number of Clusters] Again with $p=4$ and $B=2$, we let $k=(k_1, k_2)$ vary. We define clusters $\theta_{1a} = (\cos(\frac{2 \pi}{5} \cdot a), \sin(\frac{2 \pi}{5} \cdot a))'$ for $1 \leq a \leq k_1$ and similarly for the second block. All performance measures decrease as the number of clusters increase. Results are shown in Table (ref).

design[Misspecification] In this simulation, we repeat the design above using $B=1$ and $k = k_1 \cdot k_2$, the minimal number of clusters for consistent estimation using the single latent type assumption ($B=1$) considered in the literature. As expected, there is a significant power loss. Results are shown in Table (ref).

design[Block Dimension Imbalance] We let $p=12$, $B=2$ and $(k_1, k_2) = (2, 2)$. We vary the grouping of covariates, taking the first block to be $(x_{it1}, \dots, x_{itm})$ for $m \in \{1, \dots, 6\}$. Coverage-small denotes the average coverage for parameters belonging to the small block $(x_{it1}, \dots, x_{itm})$ and conversely for Coverage-large. Cluster loss-large and Cluster loss-small are defined similarly. Classification and coverage are worse for the block of smaller dimension $m$ when $m / (p-m)$ is very small, but quickly equalize as $m$ gets larger. Results are shown in Table (ref).

design[Covariate Dimension] In this simulation, we take $B = p$ (one latent variable for each covariate) and study the effect of increasing $p$. We let $k_{\ell} = 2$ for $1 \leq \ell \leq p$ and clusters $\theta_{\ell a} = \pm 1$ ($a \in \{1, 2\}$). Performance only slightly deteriorates as $p$ increases. Results are shown in Table (ref)

Model Selection

design[Model Selection - Number of Clusters] We implement the $C_p$ criterion and study its performance on the DGP in design (ref) above. We use penalty sequence $f(N, T) = \widehat{\sigma}^2 \frac{\log T}{T}$, as in section (ref). Model loss is calculated using average (over $k_1, \dots k_{B}$) distance from the truth $\|\widehat{k} - k^0\|_1 / B$. We set $k_{max} = (6, 6)$ and use $200$ independent samples. For $k^0=(2, 3)$, we estimate $\mathbb{E}\frac{\|\widehat{k}-k^0\|}{B} = 0.03$. For $k^0 = (4, 4)$, we find $\mathbb{E}\frac{\|\widehat{k}-k^0\|}{B} = 0.73$, with all estimates $\widehat{k} = (3, 3), (3, 4)$, or $(4, 3)$. We view this performance as reasonable given that the clusters are quite close in this design, though the results suggest we may be slightly over-penalizing.

Conclusion

Clustering methods have recently become popular as a way of modeling limited heterogeneity in panel data. This paper motivates a family of panel structures, nesting the standard regression clustering model, that have significant cross-sectional homogeneity but are nevertheless ill-suited to estimation by the clustering methods currently considered in the literature We propose a modified procedure that simultaneously clusters on multiple discrete latent types, significantly expanding the set of panel structures that can be accommodated by these methods. We employ Lloyd's algorithm to compute the estimator, and give consistency and asymptotic normality results for the resulting estimates.