# Clustered Covariate Regression††thanks: This paper benefited from the
valuable support of Brantly Callaway, Weige Huang, Chengye Jia, Feiyu Jiang,
Sonia Karami, Mark Robinson, Oleg Rytchkov, Vira Semenova, Rahul Singh, and
Dai Zusai. The paper is based on the second chapter of author 2’s
dissertation.

Abdul-Nasah Soale111email<EMAIL_ADDRESS>Department of Applied and
Computational Mathematics and Statistics, University of Notre Dame Emmanuel
Selorm Tsyawo222Corresponding author, email<EMAIL_ADDRESS>AIRESS &
FGSES, Université Mohammed VI Polytechnique, Morocco

(August 28, 2024)

###### Abstract

High covariate dimensionality is an increasingly occurrent phenomenon in model
estimation. A common approach to handling high-dimensionality is
regularisation, which requires sparsity of model parameters. However, sparsity
may not always be supported by economic theory or easily verified in some
empirical contexts; severe bias and misleading inference can occur. This paper
introduces a grouped parameter estimator (GPE) that circumvents this problem
by using a parameter clustering technique. The large sample properties of the
GPE hold under fairly standard conditions including a compact parameter
support that can be bounded away from zero. Monte Carlo simulations
demonstrate the excellent performance of the GPE relative to competing
estimators in terms of bias and size control. Lastly, an empirical application
of the GPE to the estimation of price and income elasticities of demand for
gasoline illustrates the practical utility of the GPE.

Keywords: clustering, high-dimension, approximation error
JEL classification: C01, C55

##  1 Introduction

High covariate dimensionality, where the covariate dimension $p$ is large
relative to the sample size $n$, is increasingly occurrent in model
estimation. High-dimensionality occurs for a variety of reasons: non-
parametric functions [3], unobserved interaction structures [24], fixed
effects in panel data [8], unit-specific model parameters in panel data [13],
reduced residual variation via more controls [5], or sub-sample analyses. High
dimensional settings are challenging for estimation and inference using
standard methods such as Ordinary Least Squares (OLS). When $p$ is large
relative to $n$, OLS, for example, is highly imprecise or infeasible when
$p>n$.

Several methods designed to tackle high dimensionality in model estimation are
regularisation-based, e.g., the (post-)LASSO (pLASSO hereafter) of [3, 6, 17,
5], and the Dantzig selector of [9, 14]. Regularisation-based estimators rely
on a(n) (approximate) sparsity condition where a single group parameter,
typically zero, is shared by a large subset of the covariates while the
remaining individual parameters are “freely” determined from the data. As an
alternative to regularisation-based estimators, clustering-based estimators,
e.g., [8, 13, 1, 15, 23, 22, 26], rely on discrete heterogeneity of the
parameters, well-separated parameter groups, and panel data.333[7], a notable
exception, allows for continuous heterogeneity and controls for approximation
error. The underlying assumptions of these methods may not always be
sufficiently borne out by economic theory or be easily verified in empirical
contexts since the data-generating process (DGP hereafter) is unobserved. For
instance, the support of a parameter vector may be bounded away from zero or
its heterogeneity may not be discrete. This paper introduces the Grouped
Parameter Estimator (GPE), which obviates sparsity and discrete heterogeneity
conditions.

The GPE incorporates parameter clustering in model estimation, where $k$
parameter groups are “freely” determined from the data. Alternatively, the GPE
can be viewed as a constrained least squares estimator with $p-k$ parameter
equality constraints that are optimally allocated from the data. The GPE
presents a number of practical advantages relative to existing estimators.
First, the GPE is a parsimonious estimator that only requires at least one
group for reliable model estimation and inference. Second, the GPE is robust
to different configurations of the unobservable parameter vector, namely, (1)
sparsity, approximate sparsity, and non-sparsity and (2) continuous, mixed,
and discrete heterogeneity. Third, the large sample properties of the GPE are
established under both fixed and increasing $p$. Lastly, this paper proposes a
simple yet consistent data-driven selection rule for the number of groups that
controls approximation error.

The remainder of the paper is laid out as follows. Section 2 presents the GPE
and the estimation algorithm. The large sample properties of the grouped
parameter estimator are established under standard regularity conditions in
Section 3. Section 4 examines the finite sample performance of the GPE
relative to competing estimators using simulation experiments. In Section 5,
we apply the GPE in analyses of price and income elasticities of demand for
gasoline. Section 6 concludes. The proofs of all theoretical results are
collected in the appendix.

##  2 The GPE

###  2.1 A population-level treatment

Let ${\bm{x}}$ be a $p$-dimensional covariate vector and $y$ denote a scalar
response. Consider the linear model

(2.1) $\displaystyle y={\bm{x}}\beta_{o}+\varepsilon,$

where $\beta_{o}\in\mathbb{R}^{p}$ is a vector of parameters, and the random
noise $\varepsilon$ satisfies $\mathbb{E}[\varepsilon]=0$,
$\mathrm{var}[\varepsilon]<\infty$, and
$\mathbb{E}[{\bm{x}}^{\prime}\varepsilon]=0$. One of the popular estimators of
$\beta_{o}$ in (2.1) is OLS. Some elements of $\beta_{o}$, however, may not be
significantly different from each other (in the statistical sense given finite
samples). For instance, if
$\beta_{oj}\approx\beta_{oj^{\prime}}\approx\beta_{oj^{\prime\prime}}$
sufficiently, one can assign a common parameter to the corresponding
$(j,j^{\prime},j^{\prime\prime})$’th covariates in ${\bm{x}}$. Moreover, $p$
may be large relative to or even larger than $n$. The foregoing present a
strong case for estimating fewer than $p$ parameters for $\beta_{o}$. More
formally, for $p\geq 2$, we consider the case where $\beta_{o}$ can be
partitioned into $k$ groups, for $k\leq p$. By (2.1), the covariates in each
group can be aggregated and assigned their respective group parameters
$\delta_{ol},\ l\in\\{1,\ldots,k\\}$. Thus, we fit (2.1) by estimating only
$k$ parameters. When $k=p$, each group becomes a singleton and OLS applies.

The proposed GPE relies on group assignments which are not known _a priori_
but can be determined from data. Let ${\bm{m}}\in\mathcal{M}_{k}$ be a binary
group assignment matrix which represents the group assignments of $\beta$ to
$k$ groups, and let $\mathcal{M}_{k}\subset\mathbb{R}^{p\times k}$ represent
the space of $p\times k$ binary matrices – details follow in Condition 2.1.
The $(j,l)$’th element of ${\bm{m}}$ equals $1$ if $\beta_{j}$ belongs to
group $l$ and zero otherwise. As each row of ${\bm{m}}$ contains a single one
and $k-1$ zeros, ${\bm{m}}^{\prime}{\bm{m}}\in\mathbb{R}^{k\times k}$ is a
diagonal matrix with non-zero diagonal entries comprising group sizes and
$\mathrm{tr}({\bm{m}}^{\prime}{\bm{m}})=p$, where $\mathrm{tr}(A)$ denotes the
trace of matrix $A$. Moreover, because each element of $\beta$ belongs to only
one group and empty groups are not allowed, the inner product of any two
columns of ${\bm{m}}$ is zero. $\beta_{o}$ can then be expressed as
$\beta_{o}={\bm{m}}_{o}\delta_{o}+b_{k}$, where
$b_{k}:=\beta_{o}-{\bm{m}}_{o}\delta_{o}$ is a non-stochastic _approximation
error_ – Lemma 2.1 below states the bound on $b_{k}$.

For $k$ fixed, denote group assignment indicators by
$g_{j}\in\\{1,\ldots,k\\}$. The GPE estimand, viz. ${\bm{m}}_{o}\delta_{o}$,
solves the following clustering problem:

$\min_{{\bm{m}}\in\mathcal{M}_{k},\
\delta\in\Delta}||{\bm{m}}\delta-\beta_{o}||^{2}=\min_{\\{\tilde{\delta},g_{1},\ldots,g_{p}\\}}\sum_{j=1}^{p}(\tilde{\delta}(g_{j})-\beta_{oj})^{2}=\sum_{j=1}^{p}(\delta_{og_{j}^{*}}-\beta_{oj})^{2}$

where
$\tilde{\delta}:=(\tilde{\delta}(1),\ldots,\tilde{\delta}(k))^{\prime}\in\Delta\subset\mathbb{R}^{k}$
is the vector of group parameters. To characterise the bound on the
approximation error $b_{k}$, we introduce the following standard assumption:

######  Assumption 2.1.

The support of the $p\times 1$ parameter vector $\beta$ is compact.

Assumption 2.1 is a standard assumption in the literature – cf. [8, Assumption
1(a)], [13, Assumption R(ii) ], and [7, Lemma 2 ]. It is useful in deriving a
bound on the approximation error $b_{k}$.

The GPE estimand ${\bm{m}}_{o}\delta_{o}$ is invariant to the relabelling of
groups. Following the literature, e.g., [8, 13], we define
${\bm{m}}\in\mathcal{M}_{k}$ up to permutations that leave
${\bm{m}}\delta_{o}$ unchanged. To avoid overly dominant groups that may
induce ill-conditioning in the matrix
${\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}$ for all
${\bm{m}}\in\mathcal{M}_{k}$, we impose the following condition.

######  Condition 2.1.

The eigenvalues of ${\bm{m}}^{\prime}{\bm{m}}$ for all
${\bm{m}}\in\mathcal{M}_{k}$ are bounded above and away from zero uniformly in
$k$.

The eigenvalues of ${\bm{m}}^{\prime}{\bm{m}}$ are group sizes induced by the
group assignment matrix ${\bm{m}}\in\mathcal{M}_{k}$. Condition 2.1 can always
be guaranteed by splitting up overly dominant groups that may cause the non-
invertibility of
$A_{m}:={\bm{m}}^{\prime}\mathbb{E}[{\bm{x}},{\bm{x}}]{\bm{m}}$. Condition 2.1
is without loss of generality as increasing the number of groups by simply
splitting an overly dominant group into, say two, does not increase the
approximation error. As Condition 2.1 is without loss of generality, it is
implicitly assumed throughout the paper.

#### Notation:

For any two sequences of non-negative numbers $\\{a_{n}:n\geq 1\\}$ and
$\\{b_{n}:n\geq 1\\}$, $a_{n}\lesssim b_{n}$ means $a_{n}\leq cb_{n}$ for some
finite $c>0$; $a_{n}\asymp b_{n}$ means $a_{n}\lesssim b_{n}$ and
$b_{n}\lesssim a_{n}$; $\rho_{\mathrm{\max}}(A)$ denotes the largest
eigenvalue of a square matrix $A$; $||\cdot||$ is the usual Euclidean norm
when applied to vectors and the spectral norm, i.e.,
$||A||:=\sqrt{\rho_{\mathrm{\max}}(A^{\prime}A)}$ when applied to matrices;
and $||\Upsilon||_{W}^{2}:=\Upsilon^{\prime}W\Upsilon$ denotes the weighted
Euclidean norm of a vector $\Upsilon$ for some conformable positive definite
matrix $W$. Define
$\mathcal{S}_{p}:=\\{\uptau\in\mathbb{R}^{p}:||\uptau||=1\\}$ as the space of
vectors with unit Euclidean norm, and let $|\mathrm{supp}(\beta_{o})|$ denote
the number of support points of $\beta_{o}$.

The following result derives the bound of the approximation error $b_{k}$ in
terms of $p$ and $k$.

######  Lemma 2.1.

Let Assumption 2.1 hold, then $||b_{k}||^{2}\asymp p/k^{2}$ if
$k<|\mathrm{supp}(\beta_{o})|$ and zero otherwise.

Lemma 2.1 shows that $||b_{k}||$ is decreasing in $k$. Lemma 2.1 is a non-
stochastic analogue of the stochastic bound on the $k$-means clustering
approximation error in [19, Theorem 5.3 ] – cf. [7, Lemma 2 ].

######  Remark 2.1.

Lemma 2.1 is a worst-case bound under continuous heterogeneity and increasing
$p$ asymptotics. Under discrete or mixed heterogeneity in $\beta_{o}$ with
$|\mathrm{supp}(\beta_{o})|<p$, zero approximation error $b_{k}$ is attainable
with $k\geq|\mathrm{supp}(\beta_{o})|$.

Next, we motivate the GPE from a regression point of view. Define the
following least squares criterion:

(2.2)
$Q_{o}({\bm{m}},\delta):=\mathbb{E}[(y-{\bm{x}}{\bm{m}}\delta)^{2}]/\sqrt{p}.$

Given ${\bm{m}}\in\mathcal{M}_{k}$,
$\delta_{o}({\bm{m}}):=({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}y]$
is a deterministic function of ${\bm{m}}$ that minimises
$Q_{o}({\bm{m}},\delta)$. The following standard assumptions are imposed in
order to characterise the identification of the GPE estimand
${\bm{m}}_{o}\delta_{o}$.

######  Assumption 2.2.

$\mathbb{E}[{\bm{x}}^{\prime}\varepsilon]=0$.

######  Assumption 2.3.

There exist constants $B>0$ and $d>2$ such that for all
$\uptau_{p}\in\mathcal{S}_{p}$, $\mathbb{E}[|{\bm{x}}\uptau_{p}|^{d}]\leq B$,
$B^{-1}\leq\mathbb{E}[|{\bm{x}}\uptau_{p}|^{2}]$, and
$\mathbb{E}[|\varepsilon|^{d}|{\bm{x}}]\leq B\ a.s.$ uniformly in $n$.

Assumption 2.2 is a weak identification assumption that admits
heteroskedasticity of arbitrary form. Assumption 2.3 requires that the
skedastic function be bounded, i.e., $\mathbb{E}[\varepsilon^{2}|{\bm{x}}]\leq
B\ a.s.$ – cf. [4, Theorem 4.1 ]. The conditional moment restriction imposed
on $\varepsilon$ in Assumption 2.3 is standard; it is slightly stronger than
Condition SM(i) of [3]. The following outlines an important implication of
Assumption 2.3 that is used throughout the text.

######  Remark 2.2.

Assumption 2.3 implies that the eigenvalues of
$\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]$ are uniformly bounded above and away
from zero as for any $\uptau_{p}\in\mathcal{S}_{p}$,
$\mathbb{E}[|{\bm{x}}\uptau_{p}|^{2}]=\uptau_{p}^{\prime}\mathbb{E}[|{\bm{x}}^{\prime}{\bm{x}}]\uptau_{p}$,

$B^{-1}\leq\underline{\uptau}_{p}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]\underline{\uptau}_{p}=\rho_{\mathrm{\min}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])\leq\uptau_{p}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]\uptau_{p}\leq\bar{\uptau}_{p}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]\bar{\uptau}_{p}=\rho_{\mathrm{\max}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])\leq
B$

where $\underline{\uptau}_{p}$ and $\bar{\uptau}_{p}$, respectively, denote
the eigenvectors associated with the smallest and largest eigenvalues of
$\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]$.

In view of Remark 2.2, Assumption 2.3 rules out perfect or very high
collinearity in ${\bm{x}}$ – cf. [3, Condition SM(i) ]. It also ensures that
${\bm{m}}\mathbb{E}[{\bm{x}},{\bm{x}}]{\bm{m}}$ is always positive-definite
for all ${\bm{m}}\in\mathcal{M}_{k}$.

Notice that $Q_{o}({\bm{m}},\delta)$ in (2.2) consists of terms involving
$b_{k}$. For the terms not depending on $b_{k}$, define the auxiliary
criterion
$\tilde{Q}_{o}({\bm{m}},\delta):=({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})/\sqrt{p}+\sigma^{2}/\sqrt{p}$
which is uniquely minimised at
$\\{{\bm{m}},\delta\\}=\\{{\bm{m}}_{o},\delta_{o}\\}$. The following result
characterises the identification of $\\{{\bm{m}}_{o},\delta_{o}\\}$ up to a
term involving the approximation error $b_{k}$.

######  Theorem 2.1.

Under Assumptions 2.2, 2.3 and 2.1,
$Q_{o}({\bm{m}},\delta)=\tilde{Q}_{o}({\bm{m}},\delta)+O_{p}(||b_{k}||)$ and
$\tilde{Q}_{o}({\bm{m}},\delta)$ is uniquely minimised at
$\\{{\bm{m}},\delta\\}=\\{{\bm{m}}_{o},\delta_{o}\\}$.

Theorem 2.1 defines the GPE estimand ${\bm{m}}_{o}\delta_{o}$ up to a bias
term $||b_{k}||$. Identification thus holds exactly when the approximation
error $b_{k}=0$ or approximately when $||b_{k}||=o(1)$.

######  Corollary 2.1.

Suppose Assumptions 2.2, 2.3 and 2.1 hold, then (a) the minimiser of
$Q_{o}({\bm{m}},\delta)$ solves the parameter clustering problem:
$\displaystyle\min_{m\in\mathcal{M}_{k},\
\delta\in\Delta}||{\bm{m}}\delta-\beta_{o}||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}$
and (b) $\displaystyle||b_{k}||^{2}\asymp\min_{m\in\mathcal{M}_{k},\
\delta\in\Delta}||{\bm{m}}\delta-\beta_{o}||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}$.

Corollary 2.1 has two important takeaways: (1) the GPE estimand solves a
(weighted) univariate $k$-means parameter clustering problem and (2) the
approximation errors from the weighted and unweighted clustering problems are
of the same order. In practice, the latter is particularly helpful as it is
exploited in the computational algorithm to lower the computational cost –
details follow in Remark 2.3.

###  2.2 Computation

The GPE estimand ${\bm{m}}_{o}\delta_{o}$ is well-defined thanks to Theorem
2.1 and Corollary 2.1 and $k$ fixed. Based on these results, we provide a
computation algorithm for the GPE. Let $\\{({\bm{x}}_{i},y_{i}),1\leq i\leq
n\\}$ be a random sample of $({\bm{x}},y)$ and $\mathbb{E}_{n}[\cdot]$ denote
the sample average. Define
$Q_{n}({\bm{m}},\delta):=\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}{\bm{m}}\delta)^{2}]/\sqrt{p}$
as the sample version of (2.2). The GPE is given by
$\hat{\beta}_{n}:=\hat{{\bm{m}}}_{n}\hat{\delta}_{n}$, where

$\\{\hat{{\bm{m}}}_{n},\hat{\delta}_{n}\\}=\operatorname*{arg\,min}_{{\bm{m}}\in\mathcal{M}_{k},\
\delta\in\Delta}Q_{n}({\bm{m}},\delta).$

In the population, computing $\\{{\bm{m}}_{o},\delta_{o}\\}$ can be
conceptualised in two steps: (1) given $\beta_{o}$, solve the clustering
problem in Corollary 2.1 for ${\bm{m}}_{o}$ and (2) solve the least squares
problem for $\delta_{o}$ given ${\bm{m}}_{o}$. In finite samples, however,
this two-step approach cannot be followed for a few reasons: (1) $\beta_{o}$
is unknown; (2) $\hat{{\bm{m}}}_{n}$ and $\hat{\delta}_{n}$ are interdependent
and cannot be computed separately; and (3) the computational cost of an
exhaustive search is prohibitive, especially for $p$ large. In light of the
foregoing, we follow the literature in applying an iterative algorithm – e.g.,
[22, 8, 13, 26]. The following algorithm is adapted from [8].

######  Algorithm 1.

  1. (a)

Fix $k\in[1,\min\\{p,n\\}-1]$, initialise $\hat{\beta}_{n}^{(0)},\
\hat{{\bm{m}}}^{(0)}$, $\hat{\delta}_{n}^{(0)}$, and counter $s=1$.

  2. (b)

For each $j\in\\{1,\ldots,p\\}$, compute $\hat{\beta}_{n,j}^{(s)}$ and update
the $j$’th row of $\hat{{\bm{m}}}^{(s)}$.

  3. (c)

Compute
$\displaystyle\hat{\delta}_{n}^{(s)}=\operatorname*{arg\,min}_{\delta\in\Delta}Q_{n}(\hat{{\bm{m}}}^{(s)},\delta)$.

  4. (d)

Check for numerical convergence, else set $s\leftarrow s+1$ and return to Step
(b).

$\hat{\beta}_{n,j}^{(s)}$ in Step (b) is computed as the simple OLS regression
of $y-{\bm{x}}_{-j}\hat{\beta}_{n,-j}^{(s-1)}$ on $x_{j}$ where
${\bm{x}}_{-j}^{\prime}$ and $\hat{\beta}_{n,-j}^{(s-1)}$ are $(p-1)\times 1$
vectors formed by removing the $j$’th elements from ${\bm{x}}^{\prime}$ and
$\hat{\beta}_{n}^{(s-1)}$, respectively. Updating the $j$’th row of
$\hat{{\bm{m}}}^{(s)}$ involves assigning $\hat{\beta}_{n,j}^{(s)}$ to the
nearest (in absolute value) of the $k$ elements in $\hat{\delta}_{n}^{(s-1)}$,
and updating the group mean. This is the step that differs from [8, Algorithm
1 ]; it is a faster yet (asymptotically) equivalent step thanks to Corollary
2.1(b). Assigning $\hat{\beta}_{n,j}^{(s)}$ to the nearest group in Step (b)
achieves the same goal up to the same order of approximation error as
verifying the group assignment of the $j$’th covariate for which the objective
function is minimised. That latter can be computationally prohibitive.

######  Remark 2.3.

Corollary 2.1 shows that the objective function is a weighted within-group sum
of squares deviations up to an approximation bias term. Adopting the
unweighted group assignment saves much computational cost as univariate
clustering algorithms are readily available. Another important difference
between our Algorithm 1 and that of [8] is that ${\bm{x}}$ can be non-binary,
thus in the case of the GPE, panel data are not _sine qua non_.

As is well known in the literature, the solution to iterative algorithms like
Algorithm 1 depend on the starting values; it is hence essential to try
different starting schemes.444See Section S.3.1 or Section S.3 in the
supplement of [7] for examples. Following the literature, e.g., [7, 630 ] or
[13, Proof of Lemma 3.1 ], our asymptotic theory focusses on the global
minimum (indexed by $k$) while abstracting away from optimisation error.

######  Remark 2.4.

A researcher may sometimes be interested in the partial effects of a finite
set of covariates and may not want to group them. Such a restriction is easily
incorporated in ${\bm{m}}$ by adding columns corresponding to singleton groups
of covariates of interest whose parameters are updated in Step (c) while group
assignments are held fixed in Step (b) across iterations. We handle the
intercept in this way.

##  3 Large Sample Properties

###  3.1 Convergence of $\hat{{\bm{m}}}_{n}$

We consider a sequence of models indexed by $n$ with observed data
$\\{D_{i},1\leq i\leq n\\}$, where $D_{i}:=({\bm{x}}_{i},y_{i})$,
$y_{i}=y_{i,n}$, and ${\bm{x}}_{i}={\bm{x}}_{i,n}$. The large sample
properties of the GPE are developed under the following additional
assumptions.

######  Assumption 3.1.

Observed data $\\{D_{i},1\leq i\leq n\\}$ are independent and identically
distributed $(iid)$ for each $n$.

######  Assumption 3.2.

$p\log(p)=o(n)$.

The dependence of the sampling process on $n$ in Assumption 3.1 allows $p$ to
grow with $n$ while allowing for fixed-$p$ asymptotics as well. Assumption 3.2
characterises the rate at which $p$ is allowed to grow with $n$.

The following provides a crucial result on the rate of convergence of
$||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||$.

######  Theorem 3.1.

Suppose Assumptions 2.1, 2.2, 3.1, 2.3 and 3.2 hold, then
$p||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||=O_{p}(||b_{k}||^{2})$.

From Theorem 3.1 above, one observes that the consistency of
$\hat{{\bm{m}}}_{n}$ depends on the magnitude of the approximation error
$b_{k}$.

######  Remark 3.1.

The rate condition in Assumption 3.2 required under the GPE is more
restrictive relative to the $\log(p)=o(n^{1/3})$ rate of, e.g., [3, 17, 5],
which allows up to an exponential growth of $p$ in $n$.

By Remark 3.1, the slower growth rate of $p$ in $n$ allowed by the GPE is the
price paid by not explicitly exploiting sparsity and model selection as in
(post)-LASSO estimators.

###  3.2 Consistency and Asymptotic Normality

For a given ${\bm{m}}\in\mathcal{M}_{k}$, the corresponding estimator is

$\hat{\delta}_{n}({\bm{m}})=({\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}].$

Define
$\hat{A}_{m}:={\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}$
and $A_{m}:={\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}$. A
consistency condition is imposed on $\hat{A}_{m_{o}}$ in the following
assumption:

######  Assumption 3.3.

$||\hat{A}_{m_{o}}-A_{m_{o}}||=o_{p}(1)$.

Assumption 3.3 is a high-level assumption; we verify it in Lemma S.2.6 under
an additional sufficient condition – see [29], [4, Lemma 6.2 ], and [12,
Corrolary 4.1 ] for applicable matrix weak law of large numbers results.

Define the terms
$\mathcal{B}_{k}=\hat{\beta}_{n}-{\bm{m}}_{o}\hat{\delta}_{n}({\bm{m}}_{o})+(\mathrm{I}_{p}+{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}])b_{k}$;
$\mathcal{A}_{k}:={\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$;
the parameter $\theta:=\uptau_{p}^{\prime}\beta_{o}$ (where dependence on
$\uptau_{p}$ is suppressed for notational ease); and

$\sigma_{\theta}:=(\uptau_{p}^{\prime}{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\varepsilon^{2}]{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\uptau_{p})^{1/2}.$

The following result provides the convergence rate of the GPE
$\hat{\beta}_{n}$ to $\beta_{o}$ and characterises its asymptotic normality
for any $\uptau_{p}\in\mathcal{S}_{p}$. It is based on the decomposition:
$\hat{\beta}_{n}-\beta_{o}=\mathcal{A}_{k}+\mathcal{B}_{k}$.

######  Theorem 3.2 (Consistency and Asymptotic Normality).

Suppose Assumptions 2.3, 2.3, 2.1, 3.3, 3.1, 2.2 and 3.2 hold, then for all
$\uptau_{p}\in\mathcal{S}_{p}$, (a)
$\uptau_{p}^{\prime}\mathcal{B}_{k}=O_{p}(n^{-1/2})+O_{p}(||b_{k}||)$; (b)
$\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o}-\mathcal{B}_{k})=O_{p}(n^{-1/2})$;
and (c)
$\sigma_{\theta}^{-1}\sqrt{n}\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o}-\mathcal{B}_{k})\xrightarrow{d}\mathcal{N}(0,1).$

Theorem 3.2 allows $p>n$ in finite samples and $(p,k)$ may be changing with
$n$.

###  3.3 Selection rule

The underlying principle of the GPE is that (2.1) be estimable with fewer than
$p$ parameters and conserve degrees of freedom. Although Lemma 2.1 could be
used to fix a choice rule on $k$ that ensures that the approximation error is
controlled, such a choice rule would be agnostic of the DGP, lack flexibility,
and poorly adapt to the DGP. For instance, under discrete, mixed, or even
continuous heterogeneity on a short diameter of $\beta_{o}$, a small $k$
should suffice to control approximation error. Thus, instead of driving
approximation error approximately to zero at the cost of a larger estimation
variance as $k$ is large and degrees of freedom are more quickly depleted, we
rather, in view of Theorem 3.2, opt for the smallest $k$ in a data-driven way
such that the approximation error is controlled – see, e.g., [7, 630 ] for a
selection rule along similar lines.

We select the number of groups per the rule

$\hat{k}:=\min_{k\geq
1}\\{k:n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})\leq\gamma\\}$

where

$\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1}):=\frac{\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k})^{2}]-\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k+1})^{2}]}{\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k+1})^{2}]},$

$\hat{\beta}_{n}^{k}$ is the GPE with $k$ groups (for notational emphasis),
and $\gamma>0$ is a user-defined constant.555We recommend $\gamma=2.7$ – see
Section S.2.7 for details. As
$\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k})^{2}]$ is decreasing in
the number of groups $k$,
$\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})$ is non-increasing in
$k$. The proposed selection rule is similar in spirit to model selection
criteria such as the Bayesian Information Criterion (BIC) without the penalty
term for overfitting – selecting the minimum for which the inequality holds
plays this role in our case – see, e.g., [26, 279 ]. The following result
shows that this simple selection rule is consistent, i.e., it selects $k$ such
that the random (approximation) error term $\mathcal{B}_{k}$ in Theorem 3.2
does not dominate the estimation error.

######  Theorem 3.3.

Fix $k$ and let $\\{a_{k}:k\geq 1\\}$ be an increasing sequence of positive
numbers such that $||\mathcal{B}_{k}||=O_{p}(n^{-1/a_{k}})$. Suppose
Assumptions 2.3, 2.3, 2.2, 2.1, 3.1 and 3.2 hold, then (a)
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})=O_{p}(1)$ if and only
if $a_{k}\geq 2$; and (b) $\displaystyle
n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})\rightarrow\infty$ as
$n\rightarrow\infty$ otherwise.

An important implication of Theorem 3.3 and Theorem 3.2(a) is that
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})=O_{p}(1)$ implies
$||\mathcal{B}_{k}||=O_{p}(n^{-1/2})$ which in turn implies
$||b_{k}||=O(n^{-1/2})$. Thus, fixing $\gamma$ to a suitable positive constant
that does not depend on $n$ ensures the approximation error does not dominate
estimation error, i.e., $||\mathcal{B}_{k}||=O_{p}(n^{-1/2})$ at $k=\hat{k}$.

It is instructive to draw parallels between $k$ and its approximate analogue –
the number $s:=|\mathrm{supp}(\beta_{o})|-1$ of non-zero parameters in the
context of the pLASSO provided zero is on the support of $\beta_{o}$.

######  Remark 3.2.

In addition to requiring that zero be on the support of $\beta_{o}$, the
pLASSO also imposes the rate condition $s^{2}(\log(p\vee n))^{2}/n=o(1)$. As
the GPE remains consistent even when $\mathrm{supp}(\beta_{o})$ is bounded
away from zero and $\hat{k}<|\mathrm{supp}(\beta_{o})|$ is allowed as long as
approximation error is controlled, the pLASSO is more restrictive in this
sense.

The choice of $k$ is user-determined in a data-driven way whereas the pLASSO
rate condition on $s$ is a restriction imposed on the support of the
unobserved $\beta_{o}$. Remarks 3.1 and 3.2 highlight the trade-off between
the pLASSO and the GPE, and the complementary roles both estimators play in
handling high-dimensionality in model estimation.

##  4 Simulation Experiment

We investigate the finite sample performance of the GPE using simulated data.
We compare the bias and rejection rates of competing estimators under (1)
under discrete, mixed, and continuous heterogeneity; (2) sparse, approximately
sparse, and non-sparse $\beta_{o}$; and (3) increasing and fixed $p$ designs.
The following estimators are considered: the GPE, pLASSO.BCCH of [3, 6] using
the recommended tuning parameter, the LASSO version (LASSO.BCCH), the LASSO
with tuning parameter chosen via 10-fold cross-validation, see, e.g., [17,
Sect. 6.1 ] (LASSO.CV), the post-LASSO version of LASSO.CV (pLASSO.CV), and
OLS.

The scalar-valued parameter of interest in this simulation exercise is
$\theta:=p^{-1/2}\sum_{j=1}^{p}\beta_{j}$. The comparison of the
aforementioned estimators is based on the following metrics: (1) mean bias
(MnB) which is $||\hat{\beta}-\beta_{o}||/\sqrt{p}$ averaged over all
simulated samples; (2) median absolute deviation (MAD) of $\hat{\theta}_{n}$;
(3) root mean squared error (RMSE) of $\hat{\theta}_{n}$; (4) the rejection
rate (Rej.) for a 5%-level test of $\mathbb{H}_{o}:\theta=\theta_{o}$; and (5)
the most frequently selected number of estimated parameters across simulated
samples ($\bar{k}$). In evaluating the rejection rate for the pLASSO.BCCH and
pLASSO.CV, parameters shrunk to zero in the first step are treated as
constants.

Define $\Phi_{.5}^{-1}(\tau_{j})$ as the $\tau_{j}$’th quantile of the normal
distribution with mean $0$ and standard deviation $0.5$ where
$\tau_{j}=0.9(j-1)/(p-1)+0.05$. The following DGPs are considered:

  1. DGP I: $y={\bm{x}}\beta+\varepsilon,\quad\beta_{oj}=\Phi_{.5}^{-1}(\tau_{j})$;

  2. DGP II: $y={\bm{x}}\beta+\varepsilon\sqrt{(1+x_{1}^{2})/2};\quad\beta_{oj}=0.7^{j-1}$;

  3. DGP III: $y={\bm{x}}\beta+\varepsilon\sqrt{(1+x_{1}^{2})/2};\quad\beta_{oj}=\mathrm{I}(j\leq 5)$; and

  4. DGP IV: $y={\bm{x}}\beta+\varepsilon\sqrt{(1+x_{1}^{2})/2};\quad\beta_{oj}=\mathrm{I}(j\leq s)\Phi_{.5}^{-1}(\tau_{j}),\quad p=100^{0.1}n^{0.9}$.

Choices of $p$ and $s$ which are functions of $n$ are rounded to the nearest
whole number. For each $(n,p)$ combination, the vector of covariates is
generated as ${\bm{x}}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$, with
$\bm{\Sigma}_{jj^{\prime}}=0.5^{|j-j^{\prime}|}$ which allows for
multicollinearity in ${\bm{x}}$. $\varepsilon\sim(\chi_{1}^{2}-1)/\sqrt{2}$ in
DGP I while $\varepsilon\sim\mathcal{N}(0,1)$ in DGPs II through IV. DGPs II
and III are adapted from [3] in order to compare the GPE to the pLASSO.BCCH
under approximate sparsity and sparsity, respectively. DGP I is intended to
examine the finite sample performance of estimators under continuous
heterogeneity, non-sparsity, homoskedasticity, and an asymmetrically
distributed $\varepsilon$. Heteroskedasticity is imposed in all other DGPs.
DGP II allows approximate sparsity with continuous heterogeneity, DGP III
allows exact sparsity on a two-point support. DGP IV allows sparsity with
limited or increasing heterogeneity. Pairwise combinations
$(n,p)\in\\{100,400,800\\}\times\\{75,150\\}$ are considered for DGPs I
through III while $p$ is allowed to grow with $n$ in DGP IV. All results are
based on 1000 simulated samples.

Table 1: DGP I

| $p=75$
---|---
$n=100$ | MnB | MAD | RMSE | Rej. | $\bar{k}$
GPE | 0.283 | 0.110 | 0.177 | 0.052 | 4
pLASSO.BCCH | 0.348 | 0.137 | 0.192 | 0.264 | 24
LASSO.BCCH | 0.402 | 0.196 | 0.285 | – | 24
pLASSO.CV | 0.207 | 0.064 | 0.104 | 0.008 | 67
LASSO.CV | 0.183 | 0.062 | 0.094 | – | 67
OLS | 0.262 | 0.075 | 0.122 | 0.001 | 75
$n=400$ |  |  |  |  |
GPE | 0.189 | 0.042 | 0.076 | 0.046 | 4
pLASSO.BCCH | 0.128 | 0.038 | 0.057 | 0.220 | 48
LASSO.BCCH | 0.176 | 0.055 | 0.080 | – | 48
pLASSO.CV | 0.070 | 0.022 | 0.033 | 0.033 | 72
LASSO.CV | 0.068 | 0.022 | 0.033 | – | 72
OLS | 0.071 | 0.022 | 0.033 | 0.027 | 75
$n=800$ |  |  |  |  |
GPE | 0.145 | 0.025 | 0.041 | 0.043 | 4
pLASSO.BCCH | 0.070 | 0.021 | 0.031 | 0.173 | 57
LASSO.BCCH | 0.122 | 0.032 | 0.050 | – | 57
pLASSO.CV | 0.047 | 0.015 | 0.022 | 0.033 | 74
LASSO.CV | 0.047 | 0.015 | 0.022 | – | 74
OLS | 0.048 | 0.015 | 0.022 | 0.032 | 75

$p=150$
---
MnB | MAD | RMSE | Rej. | $\bar{k}$
0.393 | 0.207 | 0.333 | 0.058 | 3
0.412 | 0.173 | 0.271 | 0.234 | 35
0.559 | 0.280 | 0.459 | – | 35
0.392 | 0.139 | 0.271 | 0.017 | 92
0.299 | 0.123 | 0.183 | – | 92
– | – | – | – | –
|  |  |  |
0.295 | 0.085 | 0.155 | 0.062 | 4
0.184 | 0.053 | 0.081 | 0.210 | 86
0.254 | 0.090 | 0.136 | – | 86
0.080 | 0.026 | 0.037 | 0.013 | 145
0.077 | 0.026 | 0.037 | – | 145
0.082 | 0.025 | 0.037 | 0.008 | 150
|  |  |  |
0.238 | 0.053 | 0.101 | 0.056 | 4
0.092 | 0.027 | 0.041 | 0.205 | 108
0.137 | 0.045 | 0.066 | – | 108
0.050 | 0.016 | 0.023 | 0.042 | 147
0.049 | 0.016 | 0.023 | – | 147
0.051 | 0.016 | 0.023 | 0.039 | 150

Table 2: DGP II

| $p=75$
---|---
$n=100$ | MnB | MAD | RMSE | Rej. | $\bar{k}$
GPE | 0.058 | 0.044 | 0.066 | 0.066 | 3
pLASSO.BCCH | 0.070 | 0.085 | 0.093 | 0.906 | 3
LASSO.BCCH | 0.091 | 0.227 | 0.228 | – | 3
pLASSO.CV | 0.076 | 0.042 | 0.066 | 0.253 | 12
LASSO.CV | 0.052 | 0.068 | 0.076 | – | 12
OLS | 0.265 | 0.081 | 0.123 | 0.000 | 75
$n=400$ |  |  |  |  |
GPE | 0.041 | 0.021 | 0.032 | 0.058 | 4
pLASSO.BCCH | 0.035 | 0.043 | 0.045 | 0.891 | 5
LASSO.BCCH | 0.047 | 0.116 | 0.117 | – | 5
pLASSO.CV | 0.040 | 0.021 | 0.032 | 0.216 | 15
LASSO.CV | 0.027 | 0.036 | 0.040 | – | 15
OLS | 0.071 | 0.022 | 0.034 | 0.042 | 75
$n=800$ |  |  |  |  |
GPE | 0.037 | 0.016 | 0.022 | 0.053 | 4
pLASSO.BCCH | 0.025 | 0.030 | 0.032 | 0.881 | 6
LASSO.BCCH | 0.038 | 0.075 | 0.078 | – | 6
pLASSO.CV | 0.028 | 0.016 | 0.024 | 0.220 | 17
LASSO.CV | 0.020 | 0.027 | 0.030 | – | 17
OLS | 0.048 | 0.016 | 0.022 | 0.036 | 75

$p=150$
---
MnB | MAD | RMSE | Rej. | $\bar{k}$
0.041 | 0.048 | 0.066 | 0.060 | 3
0.051 | 0.064 | 0.069 | 0.936 | 3
0.066 | 0.168 | 0.168 | – | 3
0.058 | 0.036 | 0.058 | 0.311 | 11
0.039 | 0.058 | 0.064 | – | 11
– | – | – | – | –
|  |  |  |
0.029 | 0.021 | 0.031 | 0.053 | 4
0.026 | 0.032 | 0.034 | 0.921 | 5
0.034 | 0.087 | 0.087 | – | 5
0.032 | 0.019 | 0.029 | 0.272 | 15
0.020 | 0.031 | 0.034 | – | 15
0.082 | 0.024 | 0.036 | 0.011 | 150
|  |  |  |
0.026 | 0.015 | 0.021 | 0.049 | 4
0.018 | 0.023 | 0.024 | 0.900 | 6
0.027 | 0.057 | 0.059 | – | 6
0.023 | 0.014 | 0.021 | 0.298 | 20
0.015 | 0.023 | 0.025 | – | 20
0.051 | 0.015 | 0.022 | 0.030 | 150

Table 3: DGP III

| $p=75$
---|---
$n=100$ | MnB | MAD | RMSE | Rej. | $\bar{k}$
GPE | 0.020 | 0.042 | 0.063 | 0.060 | 2
pLASSO.BCCH | 0.035 | 0.014 | 0.025 | 0.076 | 5
LASSO.BCCH | 0.116 | 0.219 | 0.224 | – | 5
pLASSO.CV | 0.068 | 0.038 | 0.062 | 0.249 | 8
LASSO.CV | 0.048 | 0.046 | 0.056 | – | 8
OLS | 0.265 | 0.081 | 0.123 | 0.000 | 75
$n=400$ |  |  |  |  |
GPE | 0.006 | 0.020 | 0.030 | 0.066 | 2
pLASSO.BCCH | 0.017 | 0.007 | 0.011 | 0.052 | 5
LASSO.BCCH | 0.050 | 0.073 | 0.081 | – | 5
pLASSO.CV | 0.034 | 0.020 | 0.030 | 0.242 | 8
LASSO.CV | 0.023 | 0.022 | 0.026 | – | 8
OLS | 0.071 | 0.022 | 0.034 | 0.042 | 75
$n=800$ |  |  |  |  |
GPE | 0.004 | 0.015 | 0.021 | 0.045 | 2
pLASSO.BCCH | 0.011 | 0.005 | 0.007 | 0.041 | 5
LASSO.BCCH | 0.038 | 0.033 | 0.044 | – | 5
pLASSO.CV | 0.024 | 0.015 | 0.022 | 0.272 | 8
LASSO.CV | 0.016 | 0.016 | 0.018 | – | 8
OLS | 0.048 | 0.016 | 0.022 | 0.036 | 75

$p=150$
---
MnB | MAD | RMSE | Rej. | $\bar{k}$
0.015 | 0.043 | 0.062 | 0.058 | 2
0.025 | 0.010 | 0.018 | 0.072 | 5
0.085 | 0.163 | 0.166 | – | 5
0.053 | 0.031 | 0.053 | 0.299 | 8
0.036 | 0.037 | 0.045 | – | 8
– | – | – | – | –
|  |  |  |
0.005 | 0.020 | 0.030 | 0.057 | 2
0.011 | 0.005 | 0.008 | 0.066 | 5
0.036 | 0.056 | 0.062 | – | 5
0.027 | 0.017 | 0.027 | 0.324 | 8
0.017 | 0.018 | 0.021 | – | 8
0.082 | 0.024 | 0.036 | 0.011 | 150
|  |  |  |
0.003 | 0.013 | 0.020 | 0.045 | 2
0.008 | 0.004 | 0.006 | 0.056 | 5
0.027 | 0.026 | 0.033 | – | 5
0.019 | 0.012 | 0.019 | 0.316 | 8
0.012 | 0.013 | 0.015 | – | 8
0.051 | 0.015 | 0.022 | 0.030 | 150

Figure 1: DGP IV (a): $s=100^{0.1}\sqrt{n}/(2\log(n\vee p))$

Figure 2: DGP IV (b): $s=5$

Figure 3: Power Curve: DGP II

Figure 4: Power Curve: DGP III

Notes: Both figures plot power curves on $\theta$ derived from DGPs II and III
at $(n,p)=(100,75)$. The deviation from the null is generated as
$\theta-\theta_{o}:=-h\sum_{j=1}^{p}\beta_{oj}/\sqrt{p},\ h\in[0,1)$.

Tables 1, 2 and 3 present results corresponding to DGPs I, II, and III.
Figures 1 and 2 plot results of two variants of DGP IV with
$s=100^{0.1}\sqrt{n}/(2\log(n\vee p))$ and $s=5$, respectively. Figures 4 and
4 present power curves on $\theta$ corresponding to DGPs II and III,
respectively. As OLS is infeasible under $p>n$, OLS results under such
configurations are not provided. Also, the rejection rates for the LASSO.BCCH
and LASSO.CV are not provided as inference on the frequentist LASSO is not
fully resolved – see [11] for a discussion. Fixing $p$, one observes a general
decrease in bias, namely, MnB, MAD, and RMSE for all estimators as $n$
increases in Tables 1, 2 and 3. Unsurprisingly, one observes a slight
deterioration in performance holding $n$ fixed with $p$ doubled. The rejection
rates are useful in showing which estimators provide reliable inference – this
is crucial in model estimation. In DGP I where $\beta_{o}$ has continuous
heterogeneity and is non-sparse, only the GPE delivers reliable inference over
all $(n,p)$ configurations.

Even under approximate sparsity (Table 2) the pLASSO.BCCH does not deliver
reliable inference while the GPE continues to perform reasonably. It is under
exact sparsity (Table 3) that inference based on the pLASSO.BCCH becomes
reliable. This appears to suggest that the performance of the LASSO-based
estimators is very sensitive to even approximately small deviations from
sparsity.666A similar observation emerges in Section S.3.2 where sparsity is
perturbed with a $O(n^{-1/2})$ term. As sparsity is not easily verifiable in
typical empirical contexts, this constitutes a major drawback of LASSO-based
estimators for model estimation. It ought to be borne in mind, however, that
the pLASSO.BCCH, among other shrinkage-based methods, is reliable in
estimating first-stage relationships in instrumental variable procedures –
see, e.g., [3, 17, 5, 21, 10]. Interestingly, the GPE chooses $k=2$ most often
in DGP III where $\beta_{o}$ has exactly $2$ support points. In DGPs I and II
where $\beta_{o}$ has continuous heterogeneity, one observes that $k$
increases weakly in $n$ but not $p$. pLASSO.CV does not deliver acceptable
rejection rates under any of the settings. This is not surprising as cross-
validation is not suitable for choosing tuning parameters for model estimation
and inference even under exact sparsity – see, e.g., [6].

Two variants of DGP IV, namely, DGP IV(a) with an increasing $s$ and DGP IV(b)
with fixed $s$ present scenarios that are favourable to the pLASSO.BCCH. Save
the GPE, pLASSO.BCCH, and OLS, no other estimator appears viable under DGP IV.
At $n=100$, one observes poor empirical size control of the pLASSO.BCCH as
evidenced by the rejection rates in Figures 1 and 2; it improves as the sample
size increases. Although feasible for sample sizes greater than $100$, OLS
performs poorly as $p$ is large relative to $n$. Figures 4 and 4 show that the
good size control of the GPE is not at the expense of power. The pLASSO.BCCH
in Figure 4, where it controls size, has more power than the GPE. This is not
surprising as the pLASSO.BCCH is based on fewer covariates that are correctly
selected in the first step. Notice in Figure 4 that although, the pLASSO.BCCH
and pLASSO.CV have more power than the GPE, they fail to control size.

In sum, the simulation exercise shows very robust performance in terms of
small bias and good size control of the GPE with its data-driven selection
rule under different empirically relevant designs viz. (1) sparse,
approximately sparse, or non-sparse $\beta_{o}$; (2) discrete, mixed, or
continuous heterogeneity; (3) homoskedastic or heteroskedastic $\varepsilon$;
and (4) skewed or symmetrically distributed $\varepsilon$. Robustness to all
the aforementioned features is an important property of an estimator as the
DGP is typically unbeknownst to the practitioner. The exercise also shows a
very high sensitivity of the pLASSO.BCCH to deviations from sparsity while the
pLASSO.CV proves unreliable for model estimation and inference.

##  5 Empirical Application

In this section, we estimate price and income elasticities of demand for
gasoline using sample selection and variable construction as in [30, 14, 27].
We use data from the 1994-1996 Canadian National Private Vehicle Use
Survey.777The data set is available on Professor Yatchew’s website:
https://economics.utoronto.ca/yatchew/. Table 4 reports summary statistics of
gas consumption, price, and income.

Table 4: Summary Statistics

Age-group | # Obs. | Gas Consumption | Price | Income
---|---|---|---|---
All | 5001 | 291.898 (273.327) | 0.567 (0.050) | 49.140 (24.878)
20-33 | 713 | 259.976 (186.457) | 0.568 (0.048) | 47.433 (22.329)
34-37 | 552 | 291.555 (238.942) | 0.569 (0.050) | 51.857 (23.802)
38-40 | 455 | 312.438 (243.983) | 0.569 (0.051) | 57.505 (25.250)
41-46 | 806 | 340.280 (291.522) | 0.566 (0.049) | 56.526 (24.670)
47-52 | 683 | 356.302 (329.025) | 0.573 (0.052) | 57.152 (25.778)
53-59 | 604 | 317.312 (307.942) | 0.565 (0.048) | 51.325 (25.997)
60-65 | 437 | 262.965 (315.203) | 0.567 (0.050) | 39.863 (20.657)
$\geq$66 | 751 | 195.908 (213.674) | 0.563 (0.049) | 32.124 (17.029)

Notes: Income is reported in thousands of the previous year’s combined annual
household pre-tax income. The table reports means and standard deviations in
parentheses.

The approach taken in constructing the covariate matrix follows [27, 14],
namely, we use $\log$ price, its square, $\log$ income, its square, distance
driven, its square, and 27 time, geographical, and household composition
dummies. The aforementioned covariates are augmented with an interaction of
the dummies with $\log$ price and its square in the first set of results
(Figure 6), and $\log$ income and its square in the second set of results
(Figure 6), respectively, resulting in $87$ covariates in each model. In this
empirical exercise, we focus on price and income elasticities by age groups,
namely, (20-33), (34-37), (38-40), (41-46), (47-52), (53-59), (60-65), and
above 65.

Figure 5: Price Elasticity of demand

Figure 6: Income Elasticity of demand

Notes: The thick black line is the estimate, 95% point-wise confidence bands
are dashed and the thick outermost lines are the 95% simultaneous confidence
bands using the [25] approach. The estimated average price elasticity of
demand on the full sample is $-0.022$ (standard error $0.227\times 10^{3}$)
while the estimated average income elasticity of demand on the full sample is
$0.063$ (standard error $0.007$).

The following conclusions can be drawn from Figure 6. First, the elasticity of
demand is strictly negative across all age groups; it is statistically
different from zero at the 5% level in all age groups. Second, the elasticity
of demand is non-monotone in age; consumers are fairly less sensitive to price
in the young (20-33) and old (above 65) age groups, and least sensitive in the
(53-59) group.

As expected, the income elasticity of demand is non-negative. Save for the
second age group where the income elasticity of demand is statistically
insignificant (point-wise) at the 5% level, income elasticity of demand is
point-wise significant in all other age groups. Although income elasticity of
demand is non-monotone in age, as in Figure 6 above, the sensitivity of demand
to income is generally higher in older age groups than in younger ones.

##  6 Conclusion

There is a growing need for robust and flexibly-adapted estimators in the face
of increasingly high covariate dimension in empirical work. This paper makes a
useful addition to the practitioner’s toolkit of high-dimensional model
estimation methods. Quite importantly, the GPE in this paper obviates sparsity
which is crucial for (post-)LASSO estimators. The GPE also flexibly adapts to
the data thanks to the simple data-driven selection rule for the number of
groups.

This paper develops the GPE of the linear model with exogenous covariates in
an $iid$ setting. It will be interesting to consider important extensions that
can enhance the GPE’s utility to the practitioner, e.g., endogeneity with
possibly many (weak) instruments, clustered or weakly dependent data, and non-
linear models with both smooth and non-smooth objective functions.

## References

  * [1] Tomohiro Ando and Jushan Bai  “Panel data models with grouped factor structure under unknown group membership”  In _Journal of Applied Econometrics_ 31.1 Wiley Online Library, 2016, pp. 163–191
  * [2] Maurice S Bartlett  “An inverse matrix adjustment arising in discriminant analysis”  In _The Annals of Mathematical Statistics_ 22.1 JSTOR, 1951, pp. 107–111
  * [3] Alexandre Belloni, Daniel Chen, Victor Chernozhukov and Christian Hansen  “Sparse models and methods for optimal instruments with an application to eminent domain”  In _Econometrica_ 80.6 Wiley Online Library, 2012, pp. 2369–2429
  * [4] Alexandre Belloni, Victor Chernozhukov, Denis Chetverikov and Kengo Kato  “Some new asymptotic theory for least squares series: Pointwise and uniform results”  In _Journal of Econometrics_ 186.2 Elsevier, 2015, pp. 345–366
  * [5] Alexandre Belloni, Victor Chernozhukov, Ivan Fernández-Val and Christian Hansen  “Program evaluation and causal inference with high-dimensional data”  In _Econometrica_ 85.1 Wiley Online Library, 2017, pp. 233–298
  * [6] Alexandre Belloni, Victor Chernozhukov and Christian Hansen  “High-dimensional methods and inference on structural and treatment effects”  In _Journal of Economic Perspectives_ 28.2, 2014, pp. 29–50
  * [7] Stéphane Bonhomme, Thibaut Lamadon and Elena Manresa  “Discretizing unobserved heterogeneity”  In _Econometrica_ 90.2 Wiley Online Library, 2022, pp. 625–643
  * [8] Stéphane Bonhomme and Elena Manresa  “Grouped patterns of heterogeneity in panel data”  In _Econometrica_ 83.3 Wiley Online Library, 2015, pp. 1147–1184
  * [9] Emmanuel Candes and Terence Tao  “The Dantzig selector: Statistical estimation when p is much larger than n”  In _The annals of Statistics_ 35.6 Institute of Mathematical Statistics, 2007, pp. 2313–2351
  * [10] Marine Carrasco and Guy Tchuente  “Regularized LIML for many instruments”  In _Journal of Econometrics_ 186.2 Elsevier, 2015, pp. 427–442
  * [11] George Casella, Malay Ghosh, Jeff Gill and Minjung Kyung  “Penalized regression, standard errors, and Bayesian lassos”  In _Bayesian Analysis_ 5.2 International Society for Bayesian Analysis, 2010, pp. 369–411  DOI: 10.1214/10-BA607
  * [12] Xiaohong Chen and Timothy M Christensen  “Optimal uniform convergence rates and asymptotic normality for series estimators under weak dependence and weak conditions”  In _Journal of Econometrics_ 188.2 Elsevier, 2015, pp. 447–465
  * [13] Xu Cheng, Frank Schorfheide and Peng Shao  “Clustering for multi-dimensional heterogeneity”  Working paper, 2021
  * [14] Victor Chernozhukov, Whitney Newey and Rahul Singh  “De-biased machine learning of global and local parameters using regularized Riesz representers”  In _arXiv preprint arXiv:1802.08667_ , 2018
  * [15] Denis Chetverikov and Elena Manresa  “Spectral and post-spectral estimators for grouped panel data models”  In _arXiv preprint arXiv:2212.13324_ , 2022
  * [16] James Davidson  “Stochastic limit theory: An introduction for econometricians”  Oxford University Press, 2021
  * [17] Max H Farrell  “Robust inference on average treatment effects with possibly more covariates than observations”  In _Journal of Econometrics_ 189.1 Elsevier, 2015, pp. 1–23
  * [18] Walter D Fisher  “On grouping for maximum homogeneity”  In _Journal of the American statistical Association_ 53.284 Taylor & Francis, 1958, pp. 789–798
  * [19] Siegfried Graf and Harald Luschgy  “Rates of convergence for the empirical quantization error”  In _Annals of probability_ JSTOR, 2002, pp. 874–897
  * [20] Bruce Hansen  “Econometrics”  Princeton University Press, 2022
  * [21] Christian Hansen and Damian Kozbur  “Instrumental variables estimation with many weak instruments using regularized JIVE”  In _Journal of Econometrics_ 182.2 Elsevier, 2014, pp. 290–308
  * [22] Chang-Ching Lin and Serena Ng  “Estimation of panel data models with parameter heterogeneity when group membership is unknown”  In _Journal of Econometric Methods_ 1.1 De Gruyter, 2012, pp. 42–55
  * [23] Ruiqi Liu, Zuofeng Shang, Yonghui Zhang and Qiankun Zhou  “Identification and estimation in panel models with overspecified number of groups”  In _Journal of econometrics_ 215.2 Elsevier, 2020, pp. 574–590
  * [24] Elena Manresa  “Estimating the structure of social interactions using panel data”  In _Unpublished Manuscript. CEMFI, Madrid_ , 2013
  * [25] José Luis Montiel Olea and Mikkel Plagborg-Møller  “Simultaneous confidence bands: Theory, implementation, and an application to SVARs”  In _Journal of Applied Econometrics_ 34.1 Wiley Online Library, 2019, pp. 1–17
  * [26] Vasilis Sarafidis and Neville Weber  “A partially heterogeneous framework for analyzing panel data”  In _Oxford Bulletin of Economics and Statistics_ 77.2 Wiley Online Library, 2015, pp. 274–296
  * [27] Vira Semenova and Victor Chernozhukov  “Debiased machine learning of conditional average treatment effects and other causal functions”  In _The Econometrics Journal_ 24.2 Oxford University Press, 2021, pp. 264–289
  * [28] Jack Sherman and W.J. Morrison  “Adjustment of an inverse matrix corresponding to changes in the elements of a given column or a given row of the original matrix”  In _Annals of Mathematical Statistics_ 20.4, 1949, pp. 621
  * [29] Joel A Tropp  “User-friendly tail bounds for sums of random matrices”  In _Foundations of computational mathematics_ 12 Springer, 2012, pp. 389–434
  * [30] Adonis Yatchew and Joungyeo Angela No  “Household gasoline demand in Canada”  In _Econometrica_ 69.6 JSTOR, 2001, pp. 1697–1709

Appendix

Abdul-Nasah Soale & Emmanuel Selorm Tsyawo

The supplementary material contains the proofs of all technical results and
useful details. The proofs of the results in the main text are first presented
in Appendix S.1, followed by supporting lemmata in Appendix S.2. Appendix S.3
presents other useful details related to the implementation of the GPE and
simulation experiments.

##  Appendix S.1 Proof of Main Theoretical Results

###  S.1.1 Proof of Lemma 2.1

Let $\tilde{\beta}_{o}(l),\ l\in\\{1,\ldots,k\\}$ denote the $l$’th group
average of $\beta_{o}$, then for $k<|\mathrm{supp}(\beta_{o})|$,

$\displaystyle||{\bm{m}}_{o}\delta_{o}-\beta_{o}||^{2}$
$\displaystyle=\sum_{j=1}^{p}\sum_{l=1}^{k}(\delta_{ol}-\beta_{oj})^{2}\cdot\mathrm{I}(g_{j}^{*}=l)$
$\displaystyle=\sum_{j=1}^{p}(\tilde{\beta}_{o}(g_{j}^{*})-\beta_{oj})^{2}\leq
p\times\max_{j}(\tilde{\beta}_{o}(g_{j}^{*})-\beta_{oj})^{2}\lesssim p/k^{2}$

by Assumption 2.1 and Lemma S.2.1.

First, since $0<||{\bm{m}}_{o}\delta_{o}-\beta_{o}||^{2}$ if
$k<|\mathrm{supp}(\beta_{o})|$,
$||{\bm{m}}_{o}\delta_{o}-\beta_{o}||^{2}\asymp p/k^{2}$ from above. Second,
$||{\bm{m}}_{o}\delta_{o}-\beta_{o}||^{2}=0$ if
$k\geq|\mathrm{supp}(\beta_{o})|$, i.e., the approximation error is exactly
zero if $k$ is at least equal to the number of support points of $\beta_{o}$.
This concludes the proof. ∎

###  S.1.2 Proof of Theorem 2.1

First, let us obtain the following useful expression of
$Q_{o}({\bm{m}},\delta)$ in terms of the approximation error
$b_{k}:=\beta_{o}-{\bm{m}}_{o}\delta_{o}$ and
$\sigma^{2}:=\mathbb{E}[\varepsilon^{2}]$. Define the auxiliary criterion
$\tilde{Q}_{o}({\bm{m}},\delta):=({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})/\sqrt{p}+\sigma^{2}/\sqrt{p}$.

$\displaystyle
Q_{o}({\bm{m}},\delta)=\mathbb{E}[(\varepsilon-{\bm{x}}({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o}-b_{k}))^{2}]/\sqrt{p}$
$\displaystyle=\sigma^{2}/\sqrt{p}+({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o}-b_{k})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o}-b_{k})/\sqrt{p}$
$\displaystyle=({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})/\sqrt{p}+\sigma^{2}/\sqrt{p}$
$\displaystyle-2({\bm{m}}\delta-{\bm{m}}_{o}\delta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]b_{k}/\sqrt{p}+b_{k}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]b_{k}/\sqrt{p}$
(S.11)
$\displaystyle=\tilde{Q}_{o}({\bm{m}},\delta)+O(||b_{k}||)+O(||b_{k}||^{2}/\sqrt{p})=\tilde{Q}_{o}({\bm{m}},\delta)+O(||b_{k}||)$

The second equality follows from Assumption 2.2 and the penultimate equality
follows from Assumptions 2.3 and 2.1. The last follows from Lemma 2.1 as
$||b_{k}||\asymp\sqrt{p}/k$ whence
$||b_{k}||^{2}/\sqrt{p}\asymp\sqrt{p}/k^{2}\leq\sqrt{p}/k$ for $k\geq 1$,
i.e., the $O(||b_{k}||)$ term dominates the $O(||b_{k}||^{2}/\sqrt{p})$ term.
Notice that by the invertibility of $\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]$
which holds under Assumption 2.3 and holding ${\bm{m}}\in\mathcal{M}_{k}$
fixed, $Q_{o}({\bm{m}},\delta)$ is minimised at

$\displaystyle\delta_{o}({\bm{m}})$
$\displaystyle=({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}y]$
$\displaystyle=({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}_{o}\delta_{o}-({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]b_{k}$
$\displaystyle=({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}_{o}\delta_{o}+O(||b_{k}||)$
$\displaystyle=:\tilde{\delta}_{o}({\bm{m}})+O(||b_{k}||)$

thanks to Assumption 2.2.
$({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]b_{k}=O(||b_{k}||)$
under Assumption 2.3 because

(S.12)
$\displaystyle||({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]^{1/2}||$
$\displaystyle=\sqrt{||({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}||^{2}}$
$\displaystyle=\sqrt{||({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}||^{2}}\lesssim
1$

for all ${\bm{m}}\in\mathcal{M}_{k}$ as there are no empty groups and
$||b_{k}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]^{1/2}||\leq||b_{k}||\cdot\sqrt{\rho_{\mathrm{\max}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])}\lesssim||b_{k}||$
by Assumption 2.3.

$\tilde{Q}_{o}({\bm{m}},\tilde{\delta}_{o}({\bm{m}}))=\delta_{o}^{\prime}\psi({\bm{m}})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]\psi({\bm{m}})\delta_{o}/\sqrt{p}+\sigma^{2}/\sqrt{p}=||\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]^{1/2}\psi({\bm{m}})\delta_{o}||^{2}/\sqrt{p}+\sigma^{2}/\sqrt{p}$

where
$\psi({\bm{m}}):={\bm{m}}({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}}_{o}-{\bm{m}}_{o}$.
Observe that ${\bm{x}}\psi({\bm{m}})\delta_{o}$ is the complement of the
linear projection of ${\bm{x}}{\bm{m}}$ on ${\bm{x}}{\bm{m}}_{o}\delta_{o}$.
As both ${\bm{x}}{\bm{m}}$ and ${\bm{x}}{\bm{m}}_{o}$ are full column rank by
construction (Assumption 2.3), ${\bm{x}}\psi({\bm{m}})\delta_{o}\neq 0$
($a.s.$) for all ${\bm{m}}\neq{\bm{m}}_{o}$ up to permutations such that
${\bm{m}}\delta_{o}={\bm{m}}_{o}\delta_{o}$. It follows that
$\displaystyle\min_{{\bm{m}}\in\mathcal{M}_{k}\setminus{\bm{m}}_{o}}\tilde{Q}_{o}({\bm{m}},\tilde{\delta}_{o}({\bm{m}}))>\sigma^{2}/\sqrt{p}$.
∎

###  S.1.3 Proof of Corollary 2.1

Part (a): Under Assumption 2.2, the expected value of the objective function
(2.2) is

$\displaystyle Q_{o}({\bm{m}},\delta)$
$\displaystyle=\mathbb{E}[(y-{\bm{x}}{\bm{m}}\delta)^{2}]/\sqrt{p}=\mathbb{E}[(\varepsilon-{\bm{x}}({\bm{m}}\delta-\beta_{o}))^{2}]/\sqrt{p}$
$\displaystyle=\mathbb{E}[\varepsilon^{2}]/\sqrt{p}+({\bm{m}}\delta-\beta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-\beta_{o})/\sqrt{p}$
$\displaystyle=\sigma^{2}/\sqrt{p}+||{\bm{m}}\delta-\beta_{o}||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}/\sqrt{p}$

which is a (weighted) clustering criterion, i.e., assigning elements of
$\beta_{o}$ into $k$ groups which minimise the squared weighted norm
$||\cdot||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}$.

Part (b): By Assumption 2.3,
$B^{-1}\leq\rho_{\mathrm{\min}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])\leq\rho_{\mathrm{\max}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])\leq
B$ whence

$\displaystyle||{\bm{m}}\delta-\beta_{o}||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}$
$\displaystyle=||{\bm{m}}\delta-\beta_{o}||^{2}\frac{({\bm{m}}\delta-\beta_{o})^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]({\bm{m}}\delta-\beta_{o})}{||{\bm{m}}\delta-\beta_{o}||^{2}}=\tilde{\uptau}_{p}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]\tilde{\uptau}_{p}\cdot||{\bm{m}}\delta-\beta_{o}||^{2}$
$\displaystyle\asymp||{\bm{m}}\delta-\beta_{o}||^{2}$

for some $\tilde{\uptau}_{p}\in\mathcal{S}_{p}$ whenever
$||{\bm{m}}\delta-\beta_{o}||\neq 0$ and
$0=||{\bm{m}}\delta-\beta_{o}||^{2}=||{\bm{m}}\delta-\beta_{o}||_{\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]}^{2}$
otherwise. The conclusion follows from Lemma 2.1. ∎

###  S.1.4 Proof of Theorem 3.1

From the DGP (2.1),

(S.13) $\displaystyle
y={\bm{x}}{\bm{m}}_{o}\delta_{o}+\varepsilon+{\bm{x}}b_{k}.$

For ${\bm{m}}\in\mathcal{M}_{k}$, define

$\hat{\delta}_{n}({\bm{m}})=\operatorname*{arg\,min}_{\delta\in\Delta}Q_{n}({\bm{m}},\delta)=({\bm{m}}^{\prime}{\bf
X}^{\prime}{\bf X}{\bm{m}})^{-1}({\bm{m}}^{\prime}{\bf X}^{\prime}Y).$

Denote the $n\times 1$ vector of residuals by
$\mathcal{E}_{n}({\bm{m}}):=Y-{\bf
X}{\bm{m}}\hat{\delta}_{n}({\bm{m}})=(\mathrm{I}_{n}-{\bf
X}{\bm{m}}({\bm{m}}^{\prime}{\bf X}^{\prime}{\bf
X}{\bm{m}})^{-1}{\bm{m}}^{\prime}{\bf X}^{\prime})Y=:\mu_{n}({\bm{m}})Y$ where
$\mathrm{I}_{n}$ is the $n\times n$ identity matrix. Substitute in the DGP
(2.1) and (S.13), respectively, to obtain
$\mathcal{E}_{n}({\bm{m}})=\mu_{n}({\bm{m}}){\bf
X}\beta_{o}+\mu_{n}({\bm{m}})\bm{\varepsilon}$ and
$\mathcal{E}_{n}({\bm{m}}_{o})=\mu_{n}({\bm{m}}_{o})(\bm{\varepsilon}+{\bf
X}b_{k})$ observing that

$\mu_{n}({\bm{m}}_{o}){\bf X}{\bm{m}}_{o}\delta_{o}=(\mathrm{I}_{n}-{\bf
X}{\bm{m}}_{o}({\bm{m}}_{o}^{\prime}{\bf X}^{\prime}{\bf
X}{\bm{m}}_{o})^{-1}{\bm{m}}_{o}^{\prime}{\bf X}^{\prime}){\bf
X}{\bm{m}}_{o}\delta_{o}={\bf X}{\bm{m}}_{o}\delta_{o}-{\bf
X}{\bm{m}}_{o}\delta_{o}=0.$

Define $\hat{\bm{\varepsilon}}_{m}:=\mu_{n}({\bm{m}})\bm{\varepsilon}$ which
is the $n\times 1$ vector of residuals from regressing ${\bf X}{\bm{m}}$ on
$\bm{\varepsilon}$. For any $\eta_{m}>0$, there exists $\eta_{\varepsilon}>0$
such that

$\displaystyle\mathbb{P}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\geq\eta_{m})\leq\mathbb{P}(Q_{n}({\bm{m}}_{o},\hat{\delta}_{n}({\bm{m}}_{o}))-Q_{n}(\hat{{\bm{m}}}_{n},\hat{\delta}_{n}(\hat{{\bm{m}}}_{n}))\geq\eta_{\varepsilon}\sqrt{p}\big{)}$
$\displaystyle=\mathbb{P}(\mathcal{E}_{n}({\bm{m}}_{o})^{\prime}\mathcal{E}_{n}({\bm{m}}_{o})-\mathcal{E}_{n}(\hat{{\bm{m}}}_{n})^{\prime}\mathcal{E}_{n}(\hat{{\bm{m}}}_{n})\geq\eta_{\varepsilon}np\big{)}$
$\displaystyle=\mathbb{P}\big{(}\hat{\bm{\varepsilon}}_{m_{o}}^{\prime}\hat{\bm{\varepsilon}}_{m_{o}}-\hat{\bm{\varepsilon}}_{\hat{m}_{n}}^{\prime}\hat{\bm{\varepsilon}}_{\hat{m}_{n}}+2b_{k}^{\prime}{\bf
X}^{\prime}\hat{\bm{\varepsilon}}_{m_{o}}+b_{k}^{\prime}{\bf
X}^{\prime}\mu_{n}({\bm{m}}_{o}){\bf X}b_{k}-2\beta_{o}^{\prime}{\bf
X}^{\prime}\hat{\bm{\varepsilon}}_{\hat{m}_{n}}\geq
np\tau_{n}(\hat{{\bm{m}}}_{n})\big{)}$
$\displaystyle\leq\mathbb{P}\big{(}O_{p}(np)^{-1}+2b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{m_{o},i}]/p+b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}/p-2\beta_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{\hat{m}_{n},i}]/p\geq\tau_{n}(\hat{{\bm{m}}}_{n})\big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{4}{\tau_{n}(\hat{{\bm{m}}}_{n})}\max\Big{\\{}O_{p}(np)^{-1},\
2b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{m_{o},i}]/p,\
b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}/p,\
-2\beta_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{\hat{m}_{n},i}]/p\Big{\\}}\geq
1\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{4}{\tau_{n}(\hat{{\bm{m}}}_{n})}|O_{p}((np)^{-1})|\geq
1\Big{)}+\mathbb{P}\Big{(}\frac{8}{\tau_{n}(\hat{{\bm{m}}}_{n})p}|b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{m_{o},i}]|\geq
1\Big{)}$
$\displaystyle+\mathbb{P}\Big{(}\frac{4}{\tau_{n}(\hat{{\bm{m}}}_{n})p}b_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}\geq
1\Big{)}+\mathbb{P}\Big{(}\frac{8}{\tau_{n}(\hat{{\bm{m}}}_{n})p}|\beta_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{\hat{m}_{n},i}]|\geq
1\Big{)}$
$\displaystyle=O((np)^{-d})+O(||b_{k}||^{2}/(p^{2}n))+O(||b_{k}||^{2}/p)+O((np)^{-d/2})$
$\displaystyle=O_{p}(||b_{k}||^{2}/p)+O_{p}((np)^{-d/2})$

where
$0<\tau_{n}(\hat{{\bm{m}}}_{n}):=\eta_{\varepsilon}+\beta_{o}^{\prime}{\bf
X}^{\prime}\mu_{n}(\hat{{\bm{m}}}_{n}){\bf X}\beta_{o}/(np)=O_{p}(1)$ under
Assumptions 2.3 and 2.1. The second inequality follows from Lemma S.2.3(a),
and the penultimate equality follows from Lemma S.2.3(b), Assumptions 2.1,
2.3, and the Markov inequality.

As $\eta_{m}>0$ is arbitrary, it can be chosen sufficiently small such that
$\sqrt{2}\geq\eta_{m}>0$, then

$\displaystyle||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||$
$\displaystyle=||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\cdot\mathrm{I}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\geq\eta_{m})+||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\cdot\mathrm{I}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||<\eta_{m})$
$\displaystyle=||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\cdot\mathrm{I}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\geq\eta_{m})$
$\displaystyle\lesssim
2\cdot\mathrm{I}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\geq\eta_{m})\ a.s.$

The second equality follows because ${\bm{m}}\in\mathcal{M}_{k}$ is a binary
matrix with $\\{0,1\\}$ entries whence
$\mathrm{I}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||<\sqrt{2})=0$ by construction.
The last line follows from the triangle inequality and the definition of
${\bm{m}}\in\mathcal{M}_{k}$ (Condition 2.1). It follows from the above and
the Markov inequality that

$\mathbb{E}[||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||]\lesssim
2\cdot\mathbb{P}(||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||\geq\eta_{m})=O_{p}(||b_{k}||^{2}/p)+O_{p}((np)^{-d/2}).$

To complete the proof, note that

$\frac{1}{(np)^{d/2}}<\frac{1}{np}\lesssim\frac{1}{p^{2}}=\Big{(}\frac{p}{p^{2}}\Big{)}/p\leq\Big{(}\frac{p}{k^{2}}\Big{)}/p\asymp||b_{k}||^{2}/p.$

The first relation follows since $d>2$ (Assumption 2.3), the second holds by
Assumption 3.2 as it implies $o(1)=\frac{p}{n}=\frac{p^{2}}{np}\lesssim 1$,
the last inequality holds by the natural condition $k\leq p$, and the last
equality follows from Lemma 2.1 under Assumption 2.1. Thus, the
$O_{p}(||b_{k}||^{2}/p)$ term dominates the $O_{p}((np)^{-d/2})$ term.

∎

###  S.1.5 Proof of Theorem 3.2

The following decomposition is useful for studying the convergence of the GPE:

(S.14)
$\displaystyle\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o})=\uptau_{p}^{\prime}\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})+\uptau_{p}^{\prime}\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))+\uptau_{p}^{\prime}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})\delta_{o}+\uptau_{p}^{\prime}b_{k}.$

From the proof of Lemma S.2.4,
$\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}=\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}+\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$.

A further decomposition of (S.1.5) gives:

$\displaystyle\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o})$
$\displaystyle=\uptau_{p}^{\prime}{\bm{m}}_{o}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})$
$\displaystyle+\uptau_{p}^{\prime}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})+\uptau_{p}^{\prime}\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))+\uptau_{p}^{\prime}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})\delta_{o}+\uptau_{p}^{\prime}b_{k}$
$\displaystyle=:\uptau_{p}^{\prime}{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+\uptau_{p}^{\prime}\mathcal{B}_{k}.$

$\mathcal{B}_{k}={\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}+(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})+\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))+(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})\delta_{o}+b_{k}$

simplifies to
$\mathcal{B}_{k}=\hat{\beta}_{n}-{\bm{m}}_{o}\hat{\delta}_{n}({\bm{m}}_{o})+(\mathrm{I}_{p}+{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}])b_{k}$
given in the main text.

(a) By Lemma S.2.4,
$\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}=\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}+\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]=O_{p}(||b_{k}||)+O_{p}(n^{-1/2})$.
As any ${\bm{m}}\in\mathcal{M}_{k}$ is a binary matrix, it simply assigns
elements of the $k\times 1$ vector $\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}$
(with replacement) to a $p$-dimensional vector, thus each element of
${\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}$
is of order at most $O_{p}(||b_{k}||)$, i.e.,
$\uptau_{p}^{\prime}{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}=O_{p}(||b_{k}||)$.

By a similar argument as above:

(1)
$\uptau_{p}^{\prime}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})=\uptau_{p}^{\prime}\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})-\uptau_{p}^{\prime}{\bm{m}}_{o}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o})=O_{p}(n^{-1/2})+O_{p}(||b_{k}||)$
by Lemma S.2.4, and the triangle inequality;

(2)
$\uptau_{p}^{\prime}\hat{{\bm{m}}}_{n}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))=O_{p}(||b_{k}||^{2}/\sqrt{p})$
by Lemma S.2.5; and

(3)
$\uptau_{p}^{\prime}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})\delta_{o}=O_{p}(||b_{k}||^{2}/p)$
by Theorem 3.1 and Assumption 2.1 since elements of $\delta_{o}$ are averages
of sub-vectors of $\beta_{o}$. Putting terms together,
$\uptau_{p}^{\prime}\mathcal{B}_{k}=O_{p}(n^{-1/2})+O_{p}(||b_{k}||)+O_{p}(||b_{k}||^{2}/\sqrt{p})=O_{p}(n^{-1/2})+O_{p}(||b_{k}||)$
by (S.1.2).

(b) This part follows from part (a) above and Lemma S.2.4 as

$\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o}-\mathcal{B}_{k})=\uptau_{p}^{\prime}{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]=O_{p}(n^{-1/2}).$

(c) From part (b) above, Assumption 3.3, the continuity of the inverse at a
non-singular matrix, and the continuous mapping theorem,

$\sqrt{n}\uptau_{p}^{\prime}(\hat{\beta}_{n}-\beta_{o}-\mathcal{B}_{k})=\uptau_{p}^{\prime}{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}(\sqrt{n}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}])+o_{p}(1).$

Under Assumptions 3.1 and 2.2,

$\mathrm{var}[\uptau_{p}^{\prime}{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}(\sqrt{n}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}])]=\uptau_{p}^{\prime}{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\varepsilon^{2}]{\bm{m}}_{o}A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\uptau_{p}=\sigma_{\theta}^{2}.$

The conclusion follows from Lemma S.2.4. ∎

###  S.1.6 Proof of Theorem 3.3

Recall the decomposition:
$\hat{\beta}_{n}-\beta_{o}=\mathcal{A}_{k}+\mathcal{B}_{k}$ used in Theorem
3.2, and observe that under the DGP (2.1), the average squared GPE residuals
can be expressed as

(S.15) $\displaystyle\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n})^{2}]$
$\displaystyle=\mathbb{E}_{n}[\varepsilon_{i}^{2}]-2(\hat{\beta}_{n}-\beta_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+(\hat{\beta}_{n}-\beta_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}](\hat{\beta}_{n}-\beta_{o})$
$\displaystyle=\mathbb{E}_{n}[\varepsilon_{i}^{2}]-2\mathcal{A}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+\mathcal{A}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{A}_{k}+2\mathcal{B}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{A}_{k}-2\mathcal{B}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle+\mathcal{B}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{B}_{k}$
$\displaystyle=\mathbb{E}_{n}[\varepsilon_{i}^{2}]-\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle-2\mathcal{B}_{k}^{\prime}(\mathrm{I}_{p}-\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime})\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+\mathcal{B}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{B}_{k}.$

It follows from (S.15) that

(S.16)
$\displaystyle\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k})^{2}]-\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k+1})^{2}]$
$\displaystyle=\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{(k+1)}\hat{A}_{m_{o}^{k+1}}^{-1}{\bm{m}}_{o}^{(k+1)\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]-\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{k}\hat{A}_{m_{o}^{k}}^{-1}{\bm{m}}_{o}^{k\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle-2\Big{(}\mathcal{B}_{k}^{\prime}\big{(}\mathrm{I}_{p}-\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}^{k}\hat{A}_{m_{o}^{k}}^{-1}{\bm{m}}_{o}^{k\prime}\big{)}-\mathcal{B}_{k+1}^{\prime}\big{(}\mathrm{I}_{p}-\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}^{(k+1)}\hat{A}_{m_{o}^{k+1}}^{-1}{\bm{m}}_{o}^{(k+1)\prime}\big{)}\Big{)}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle+\mathcal{B}_{k}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{B}_{k}-\mathcal{B}_{k+1}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\mathcal{B}_{k+1}$
$\displaystyle=O_{p}(n^{-1})+O_{p}(n^{-1/2}\mathbb{E}[||\mathcal{B}_{k}||])+O_{p}(\mathbb{E}[||\mathcal{B}_{k}||^{2}])$

and
$\displaystyle\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\hat{\beta}_{n}^{k+1})^{2}]\leq\sup_{\beta}\mathbb{E}_{n}[(y_{i}-{\bm{x}}_{i}\beta)^{2}]=O_{p}(1)$
under Assumptions 2.1, 2.3, 2.2 and 3.1. Under Assumptions 2.3 and 2.2 the
first term is always $O_{p}(n^{-1})$ by Lemma S.2.2. Putting together the
foregoing implies
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})=O_{p}(1)$ if and only
if $||\mathcal{B}_{k}||=O_{p}(n^{-1/a_{k}})=O_{p}(n^{-1/2})$ for $a_{k}\geq
2$. The second term of (S.16) is dominant if $a_{k}>2$, but the third term is
dominant when the $a_{k}<2$.888Both are of the same order at $a_{k}=2$. This
implies
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})\rightarrow\infty$ as
$n\rightarrow\infty$ when $a_{k}<2$ since the third term is positive $a.s.$ as
$n\rightarrow\infty$ under Assumption 3.2. ∎

##  Appendix S.2 Useful Lemmata

Lemmata contained in this section are used in the proofs of technical results
in the main text.

###  S.2.1 Convergence rate: non-stochastic k-means clustering

Consider the kmeans-clustering problem for a non-stochastic sequence
$\mathscr{B}_{p}=\\{b_{j}:1\leq j\leq p\\}$ in $\mathbb{R}$. The problem
involves partitioning $\mathscr{B}_{p}$ into $k$ groups
$\bm{S}=\\{S_{1},\ldots,S_{k}\\}$, $k\geq 1$, in order to minimise the within-
group sum of squared deviations:

$\\{g_{1}^{*},\ldots,g_{p}^{*}\\}=\displaystyle\operatorname*{arg\,min}_{\\{g_{1},\ldots,g_{p}\\}}\sum_{j=1}^{p}(b_{j}-\tilde{b}(g_{j}))^{2}$

where $\tilde{b}(l)$ is the mean of elements $b_{j}$ in group
$S_{l}:=\\{b_{j}\in\mathscr{B}_{p}:g_{j}=l\\}$, and
$g_{j}\in\\{1,\ldots,k\\}$. Let $\\{S_{l}^{*},1\leq l\leq k\\}$ denote the
sets of optimal group assignments and $|S_{l}|$ denote the cardinality of
$S_{l}$. The following lemma establishes a worst-case bound on the maximum
within-group deviation and the average squared deviation in terms of $k$.

######  Lemma S.2.1.

Suppose $\mathscr{B}_{p}=\\{b_{j}:j=1,\ldots,p\\}$ has a compact support in
$\mathbb{R}$,
$\displaystyle\\{g_{1}^{*},\ldots,g_{p}^{*}\\}:=\operatorname*{arg\,min}_{\\{g_{1},\ldots,g_{p}\\}}\sum_{j=1}^{p}(b_{j}-\tilde{b}(g_{j}))^{2}$
are the optimal group assignments, group sizes are bounded above and away from
zero uniformly in $k$ where $k<|\mathrm{supp}(\mathscr{B}_{p})|$, i.e., $k$ is
less than the number of support points in $\mathscr{B}_{p}$, and
$g_{j}\in\\{1,\ldots,k\\}$, then (a) $\displaystyle k\max_{1\leq j\leq
p}|b_{j}-\tilde{b}(g_{j}^{*})|\lesssim 1$ and (b)
$\displaystyle\min_{\\{g_{1},\ldots,g_{p}\\}}\frac{1}{p}\sum_{j=1}^{p}(b_{j}-\tilde{b}(g_{j}))^{2}\lesssim
1/k^{2}$.

#### Proof:

$\displaystyle\max_{1\leq j\leq p}|b_{j}|<\infty$ by the boundedness of the
support of $\mathscr{B}_{p}$, hence let $C_{b}$ define an upper-bound on the
length of the support of $\mathscr{B}_{p}$, i.e., $\displaystyle C_{b}\geq
2\max_{1\leq j\leq p}|b_{j}|$ uniformly in $p$.

$\displaystyle\begin{split}\min_{\\{g_{1},\ldots,g_{p}\\}}\frac{1}{p}\sum_{j=1}^{p}(b_{j}-\tilde{b}(g_{j}))^{2}&=\min_{\bm{S}}\frac{1}{p}\sum_{l=1}^{k}\sum_{b\in
S_{l}}(b-\tilde{b}(l))^{2}\\\ &=\frac{1}{p}\sum_{l=1}^{k}\sum_{b\in
S_{l}^{*}}(b-\tilde{b}(l))^{2}\leq\sum_{l=1}^{k}\frac{|S_{l}^{*}|}{p}\max_{b\in
S_{l}^{*}}|b-\tilde{b}(l)|^{2}\\\
&=\sum_{l=1}^{k}\tilde{p}_{l}\tilde{r}_{l}^{2}\end{split}$

where
$\displaystyle\\{S_{1}^{*},\ldots,S_{k}^{*}\\}=\operatorname*{arg\,min}_{\bm{S}}\frac{1}{p}\sum_{l=1}^{k}\sum_{b\in
S_{l}}(b-\tilde{b}(l))^{2}$, $\ \tilde{p}_{l}:=\frac{|S_{l}^{*}|}{p}$, and
$\displaystyle\tilde{r}_{l}:=\max_{b\in S_{l}^{*}}|b-\tilde{b}(l)|$ is the
maximum within-group deviation from the group mean.
$\sum_{l=1}^{k}\tilde{p}_{l}=1$ such that $\tilde{p}_{l}>0$ for each
$l=1,\ldots,k$ as group sizes are bounded away from zero.

Let the diameter of group $S_{l}^{*}$ be
$\displaystyle\mathrm{diam}(S_{l}^{*})=\max_{b_{1},b_{2}\in
S_{l}^{*}}|b_{1}-b_{2}|$. At the minimised sum of within-group sum of squared
deviations, the sum of group diameters cannot exceed the length of the
support, i.e., $\displaystyle\sum_{l=1}^{k}\mathrm{diam}(S_{l}^{*})\leq
C_{b}$, because groups do not overlap – see, e.g., [18] on the contiguity of
least squares partitions. Since maximum within-group deviation from the group
mean cannot exceed the diameter of the group, i.e., $\displaystyle\max_{1\leq
l\leq k}\big{(}\tilde{r}_{l}-\mathrm{diam}(S_{l}^{*})\big{)}\leq 0$, it
follows that $\sum_{l=1}^{k}\tilde{r}_{l}\leq C_{b}$.

Define a parameter vector
$\bm{r}=(r_{1},\ldots,r_{k})^{\prime}\in\mathbb{R}_{+}^{k}$, then

$\sum_{l=1}^{k}\tilde{p}_{l}\tilde{r}_{l}^{2}\leq\max_{\bm{r}}\sum_{l=1}^{k}\tilde{p}_{l}r_{l}^{2}\text{
subject to }\sum_{l=1}^{k}r_{l}\leq C_{b}.$

The Lagrangean corresponding to the constrained maximisation problem above is
given by

$\mathcal{L}(\bm{r},\lambda)=\sum_{l=1}^{k}\tilde{p}_{l}r_{l}^{2}-\lambda(\sum_{l=1}^{k}r_{l}-C_{b})$

where $\lambda\geq 0$ is the Lagrange multiplier. Taking partial derivatives
gives

(S.21) $\displaystyle\frac{\partial\mathcal{L}(\bm{r})}{\partial r_{l}}$
$\displaystyle:=2\tilde{p}_{l}r_{l}-\lambda,\ l=1,\ldots,k\text{ and}$ (S.22)
$\displaystyle\frac{\partial\mathcal{L}(\bm{r})}{\partial\lambda}$
$\displaystyle:=-(\sum_{l=1}^{k}r_{l}-C_{b}).$

The case of $\lambda=0$ implies the maximiser of $\mathcal{L}(\bm{r})$ is
$r_{l}^{*}=0$ for each $l=1,\ldots,k$ from (S.21) but this rather minimises
the objective function. This case is hence ruled out. For $\lambda>0$, (S.22)
gives $\sum_{l=1}^{k}r_{l}-C_{b}=0$ which is equivalently expressed as
$r_{k}=C_{b}-\sum_{l=1}^{k-1}r_{l}$. Setting the set of equations (S.21) to
zero results in $2\tilde{p}_{l}r_{l}=\lambda$ for each $l=1,\ldots,k$ and
$\tilde{p}_{l}r_{l}=\tilde{p}_{k}r_{k}$ for each $l=1,\ldots,k-1$ with
$r_{k}=C_{b}-\sum_{l^{\prime}=1}^{k-1}r_{l^{\prime}}$. Substituting
$r_{k}=C_{b}-\sum_{l^{\prime}=1}^{k-1}r_{l^{\prime}}$ into the $k-1$ equations
gives
$\tilde{p}_{l}r_{l}=\tilde{p}_{k}(C_{b}-\sum_{l^{\prime}=1}^{k-1}r_{l^{\prime}}),\
l=1,\ldots,k-1$. Dividing through by $\tilde{p}_{k}$ and rearranging terms,
$\bm{A}\bm{r}_{-k}=\bm{\nu}C_{b}$ where

$\bm{A}=\left(\begin{array}[]{cccc}1+\frac{\tilde{p}_{1}}{\tilde{p}_{k}}&1&\ldots&1\\\
\vdots&\ddots&\ldots&\vdots\\\
1&1&\ldots&1+\frac{\tilde{p}_{k-1}}{\tilde{p}_{k}}\\\ \end{array}\right)$

, $\bm{r}_{-k}=(r_{1},\ldots,r_{k-1})^{\prime}$, and $\bm{\nu}$ is a
$(k-1)\times 1$ vector of ones. $\bm{A}\bm{r}_{-k}=\bm{\nu}C_{b}$ is solved by
$\bm{r}_{-k}^{*}:=\bm{A}^{-1}\bm{\nu}C_{b}$. The next step is to obtain the
explicit expression for $\bm{A}^{-1}$.

Note that $\bm{A}=\bm{D}+\bm{\nu}\bm{\nu}^{\prime}$ where $\bm{D}$ is a
diagonal matrix with diagonal entries
$\\{D_{ll}:=\frac{\tilde{p}_{l}}{\tilde{p}_{k}}$, $1\leq l\leq k-1\\}$.
Applying the formula of [28], see e.g., [2, eqn. 2],

$(\bm{A}^{-1})_{ll^{\prime}}=\begin{cases}\frac{\tilde{p}_{k}}{\tilde{p}_{l}}-\frac{\tilde{p}_{k}^{2}}{q\tilde{p}_{l}^{2}}&\text{
if }l=l^{\prime}\\\
-\frac{\tilde{p}_{k}^{2}}{q\tilde{p}_{l}\tilde{p}_{l^{\prime}}}&\text{ if
}l\not=l^{\prime}\end{cases}$

where $q=1+\sum_{\iota=1}^{k-1}\frac{\tilde{p}_{k}}{\tilde{p}_{\iota}}$. It
hence follows that for each $l\in\\{1,\ldots,k-1\\}$,

$\displaystyle r_{l}^{*}$
$\displaystyle=C_{b}\sum_{l^{\prime}=1}^{k-1}(\bm{A}^{-1})_{ll^{\prime}}=C_{b}\Big{(}\frac{\tilde{p}_{k}}{\tilde{p}_{l}}-\frac{\tilde{p}_{k}^{2}}{q\tilde{p}_{l}^{2}}-\sum_{l^{\prime}\not=l}^{k-1}\frac{\tilde{p}_{k}^{2}}{q\tilde{p}_{l}\tilde{p}_{l^{\prime}}}\Big{)}=C_{b}\frac{\tilde{p}_{k}}{q\tilde{p}_{l}}\Big{(}q-\sum_{l^{\prime}=1}^{k-1}\frac{\tilde{p}_{k}}{\tilde{p}_{l^{\prime}}}\Big{)}$
$\displaystyle=C_{b}\frac{\tilde{p}_{k}}{(1+\sum_{\iota=1}^{k-1}\frac{\tilde{p}_{k}}{\tilde{p}_{\iota}})\tilde{p}_{l}}=C_{b}\frac{1}{\sum_{\iota=1}^{k}\frac{\tilde{p}_{l}}{\tilde{p}_{\iota}}}=C_{b}\frac{1}{\sum_{\iota=1}^{k}\frac{|S_{l}^{*}|}{|S_{\iota}^{*}|}}\lesssim
1/k$

under conditions $C_{b}<\infty$ and the uniform boundedness of group sizes
$\\{|S_{l}^{*}|,1\leq l\leq k\\}$ above and away from zero. Also, observe that
$r_{k}^{*}=C_{b}-\sum_{l=1}^{k-1}r_{l}^{*}=C_{b}\frac{1}{\sum_{\iota=1}^{k}\frac{\tilde{p}_{k}}{\tilde{p}_{\iota}}}=C_{b}\frac{1}{\sum_{\iota=1}^{k}\frac{|S_{k}^{*}|}{|S_{\iota}^{*}|}}\lesssim
1/k$. Thus, $r_{l}^{*}\lesssim 1/k$ for all $l\in\\{1,\ldots,k\\}$, and that
concludes the proof of the first part.

(b) From part (a), note that
$\displaystyle\min_{\\{g_{1},\ldots,g_{p}\\}}\frac{1}{p}\sum_{j=1}^{p}(b_{j}-\tilde{b}(g_{j}))^{2}\leq\sum_{l=1}^{k}\tilde{p}_{l}\tilde{r}_{l}^{2}\leq\max_{1\leq
l\leq k}\tilde{r}_{l}^{2}\sum_{l=1}^{k}\tilde{p}_{l}=\max_{1\leq l\leq
k}\tilde{r}_{l}^{2}\lesssim 1/k^{2}$ noting that
$\sum_{l=1}^{k}\tilde{p}_{l}=1$.

∎

###  S.2.2 Triangular array CLT

Define the triangular array $\\{S_{in},1\leq i\leq n,n\in\mathbb{N}\\}$ where
$S_{in}=\tilde{S}_{in}/s_{n}$,
$\tilde{S}_{in}=\uptau_{p}^{\prime}{\bm{x}}_{i}^{\prime}\varepsilon_{i}$,
$s_{n}^{2}=\sum_{i=1}^{n}\uptau_{p}^{\prime}\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}\varepsilon_{i}^{2}]\uptau_{p}$,
and $\uptau_{p}\in\mathcal{S}_{p}$.

######  Lemma S.2.2.

Under Assumptions 3.1, 2.2 and 2.3,
$\sum_{i=1}^{n}S_{in}\xrightarrow{d}\mathcal{N}(0,1)$.

#### Proof:

First,
$\mathbb{E}[\tilde{S}_{in}]=\uptau_{p}^{\prime}\mathbb{E}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]=0$
by Assumption 2.2 for each $i=1,\ldots,n$. Second, under independent sampling
(Assumption 3.1),
$\mathbb{E}[(\sum_{i=1}^{n}S_{in})^{2}]=\sum_{i=1}^{n}\mathbb{E}[S_{in}^{2}]=s_{n}^{2}/s_{n}^{2}=1$.
Third, it suffices to verify the Lyapunov condition
$\displaystyle\lim_{n\rightarrow\infty}\sum_{i=1}^{n}\mathbb{E}[|S_{in}|^{d}]=0$
for $d>2$ in order to complete the proof.

By the law of iterated expectations (LIE) and Assumption 2.3

(S.23)
$\mathbb{E}[\tilde{S}_{in}^{2}]=\uptau_{p}^{\prime}\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}\varepsilon_{i}^{2}]\uptau_{p}=\uptau_{p}^{\prime}\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}\mathbb{E}[\varepsilon_{i}^{2}|{\bm{x}}_{i}]]\uptau_{p}\leq
B\rho_{\mathrm{\max}}(\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}])\lesssim 1.$

This implies $s_{n}=O(n^{1/2})$. Since
$\mathbb{E}[|\uptau_{p}^{\prime}{\bm{x}}_{i}^{\prime}|^{d}]\leq B$ by
Assumption 2.3, it follows from the LIE that
$\mathbb{E}[|\tilde{S}_{in}|^{d}]=\mathbb{E}[|\uptau_{p}^{\prime}{\bm{x}}_{i}^{\prime}|^{d}\mathbb{E}[|\varepsilon_{i}|^{d}|{\bm{x}}_{i}]]\leq
B^{2}<\infty$. Thus, $\sum_{i=1}^{n}\mathbb{E}[|\tilde{S}_{in}|^{d}]=O(n)$.
Combining the foregoing,
$\sum_{i=1}^{n}\mathbb{E}[|S_{in}|^{d}]=s_{n}^{-d}\sum_{i=1}^{n}\mathbb{E}[|\tilde{S}_{in}|^{d}]=O(n^{-(d/2-1)})=o(1)$
since $d>2$ (Assumption 2.3). The conclusion follows from [16, Theorems 24.6
and 24.11 ].

∎

###  S.2.3 Auxiliary Result I

The following lemma is used in the proof of Theorem 3.1. Let
$\hat{\varepsilon}_{m_{o},i}$ and $\hat{\varepsilon}_{m,i}$, respectively,
denote residuals from regressing ${\bf X}{\bm{m}}$ and ${\bf X}{\bm{m}}_{o}$
on $\bm{\varepsilon}$

######  Lemma S.2.3.

Let $\hat{\varepsilon}_{i}$ be the $i$’th residual from regressing ${\bf
X}{\bm{m}}$ on $\bm{\varepsilon}$, then under Assumptions 2.2, 2.3 and 3.1,
(a)
$\sum_{i=1}^{n}(\hat{\varepsilon}_{m_{o},i}^{2}-\hat{\varepsilon}_{m,i}^{2})\leq
O_{p}(1)$; and
(b)$\sqrt{n}\uptau_{p}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{m,i}]=\sqrt{n}\uptau_{p}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+O_{p}(1)$
for any $\uptau_{p}\in\mathcal{S}_{p}$ and ${\bm{m}}\in\mathcal{M}_{k}$.

#### Proof:

(a) Define
$\hat{\vartheta}_{m,n}:=({\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}})^{-1}{\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$.
Under the orthogonality condition (Assumption 2.2), the resulting parameter
estimates, say $\hat{\vartheta}_{m_{o},n}$ and $\hat{\vartheta}_{m,n}$, from
both regressions converge to zero at the $O_{p}(n^{-1/2})$ rate, i.e.,
$||\hat{\vartheta}_{m_{o},n}||=O_{p}(n^{-1/2})$ and
$||\hat{\vartheta}_{m,n}||=O_{p}(n^{-1/2})$. To see this, observe that under
Assumptions 2.3 and 3.1

$||n\mathrm{var}[\hat{\vartheta}_{m,n}]||\rightarrow||({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\mathbb{E}[\varepsilon^{2}|{\bm{x}}]]{\bm{m}})({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}||\leq
B||({\bm{m}}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]{\bm{m}})^{-1}||$

which is finite under Assumption 2.3 and that there are no empty groups.
Further,

$\displaystyle\sum_{i=1}^{n}(\hat{\varepsilon}_{m_{o},i}^{2}-\hat{\varepsilon}_{m,i}^{2})$
$\displaystyle=n\hat{\vartheta}_{m_{o},n}^{\prime}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}\hat{\vartheta}_{m_{o},n}-n\hat{\vartheta}_{m,n}^{\prime}{\bm{m}}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}\hat{\vartheta}_{m,n}-2n(\hat{\vartheta}_{m_{o},n}^{\prime}{\bm{m}}_{o}-\hat{\vartheta}_{m,n}{\bm{m}}^{\prime})\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle\leq
n\hat{\vartheta}_{m_{o},n}^{\prime}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}\hat{\vartheta}_{m_{o},n}-2\sqrt{n}(\hat{\vartheta}_{m_{o},n}^{\prime}{\bm{m}}_{o}-\hat{\vartheta}_{m,n}{\bm{m}}^{\prime})(\sqrt{n}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}])$
$\displaystyle=O_{p}(1)$

under Assumptions 3.1, 2.2 and 2.3. The inequality holds because
$\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]$ is positive semi-definite.

(b) Under Assumptions 2.3 and 2.2, it follows from part (a) above that
$\uptau_{p}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\hat{\varepsilon}_{m,i}]=\uptau_{p}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]-\uptau_{p}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}\hat{\vartheta}_{m,n}=O_{p}(n^{-1/2})$.

∎

###  S.2.4 Auxiliary Result II

The following is essential, in addition to Theorem 3.1, in establishing the
rate of convergence of the GPE
$\hat{\beta}_{n}:=\hat{{\bm{m}}}_{n}\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})$ and
characterising its asymptotic normality. Define the bias term
$\hat{B}_{k}:=\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}$
and the covariance matrix
$V_{\delta}:=A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\varepsilon^{2}]{\bm{m}}_{o}A_{m_{o}}^{-1}$.

######  Lemma S.2.4.

Let Assumptions 3.1, 2.2, 2.3 and 3.3 hold, then,

$\displaystyle\text{(a)
}||\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}||=O_{p}(n^{-1/2})+O_{p}(||b_{k}||)\text{
and (b)
}\sqrt{n}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}-\hat{B}_{k})\xrightarrow{d}\mathcal{N}(0,V_{\delta}).$

#### Proof:

(a) Recall
$y_{i}={\bm{x}}_{i}{\bm{m}}_{o}\delta_{o}+{\bm{x}}_{i}b_{k}+\varepsilon_{i}$.
Thus,

$\displaystyle\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}$
$\displaystyle=\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}+\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle=:\hat{B}_{k}+\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}].$

The conclusion of part (a) of the lemma follows by Assumption 2.3,
$\hat{A}_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]b_{k}=O_{p}(||b_{k}||)$
– see (S.12) – and that
$n^{1/2}{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]=O_{p}(1)$
under the conditions of Lemma S.2.2.

(b) The second part of the proof follows from the first part, the Cramér-Wold
device, the continuity of the inverse at a non-singular matrix, the continuous
mapping theorem, Assumptions 3.1 and 2.2, and Lemma S.2.2, i.e.,

$\sqrt{n}(\hat{\delta}_{n}({\bm{m}}_{o})-\delta_{o}-\hat{B}_{k})=A_{m_{o}}^{-1}{\bm{m}}_{o}^{\prime}\sqrt{n}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+o_{p}(1)\xrightarrow{p}\mathcal{N}(0,V_{\delta}).$

∎

###  S.2.5 Auxiliary Result III

######  Lemma S.2.5.

Under Assumptions 2.1, 2.2, 3.1, 2.3 and 3.2,
$\uptau_{k}^{\prime}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))=O_{p}(||b_{k}||^{2}/\sqrt{p})$
for all $\uptau_{k}\in\mathcal{S}_{k}$.

#### Proof:

Fixing $k$ for both ${\bm{m}}_{o}$ and $\hat{{\bm{m}}}_{n}$, it is easy to see
that $\hat{{\bm{m}}}_{n}$ is a row permutation of the rows of ${\bm{m}}_{o}$.
Therefore,
$\hat{{\bm{m}}}_{n}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\hat{{\bm{m}}}_{n}={\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}+\bm{\zeta}_{m}$,
where $\bm{\zeta}_{m}$ is a $k\times k$ perturbation matrix of
${\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o}$,
namely,

(S.24) $\displaystyle\bm{\zeta}_{m}$
$\displaystyle=(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\hat{{\bm{m}}}_{n}+{\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}](\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})$

Notice that
$({\bm{m}}_{o}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]{\bm{m}}_{o})^{-1}-(\hat{{\bm{m}}}_{n}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\hat{{\bm{m}}}_{n})^{-1}=\hat{A}_{m_{o}}^{-1}-\hat{A}_{\hat{m}_{n}}^{-1}=\hat{A}_{\hat{m}_{n}}^{-1}\bm{\zeta}_{m}\hat{A}_{m_{o}}^{-1}$.
One then obtains

(S.25)
$\displaystyle\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o})$
$\displaystyle=-\hat{A}_{\hat{m}_{n}}^{-1}\bm{\zeta}_{m}\hat{A}_{m_{o}}^{-1}\hat{{\bm{m}}}_{n}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}]+\hat{A}_{m_{o}}^{-1}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}]$
$\displaystyle=-\hat{A}_{\hat{m}_{n}}^{-1}\bm{\zeta}_{m}\hat{A}_{m_{o}}^{-1}\hat{{\bm{m}}}_{n}^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}]+\big{(}\mathrm{I}_{k}-\hat{A}_{\hat{m}_{n}}^{-1}\bm{\zeta}_{m}\big{)}\hat{A}_{m_{o}}^{-1}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}].$

Note that
$\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}]=\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]\beta_{o}+\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]=O_{p}(\sqrt{p})+O_{p}(n^{-1/2})=O_{p}(\sqrt{p})$
under Assumptions 2.3, 2.2, 2.1 and 3.1.

First, since
$||\hat{A}_{m_{o}}^{-1}||=O_{p}(1)=||\hat{A}_{\hat{m}_{n}}^{-1}||$ under
Assumption 2.3 and that no group is empty, it follows from (S.24) that
$||\hat{A}_{\hat{m}_{n}}^{-1}\bm{\zeta}_{m}\hat{A}_{m_{o}}^{-1}||=O_{p}(\mathbb{E}[||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||])$.
Therefore, the first summand of
$\uptau_{k}^{\prime}(\hat{\delta}_{n}(\hat{{\bm{m}}}_{n})-\hat{\delta}_{n}({\bm{m}}_{o}))$
is of the stochastic order
$O_{p}(\sqrt{p}\mathbb{E}[||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||])$. Second, by a
similar token,
$\hat{A}_{m_{o}}^{-1}(\hat{{\bm{m}}}_{n}-{\bm{m}}_{o})^{\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}y_{i}]=O_{p}(\sqrt{p}\mathbb{E}[||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||])$.
The conclusion follows from Theorem 3.1 as
$O_{p}(\sqrt{p}\mathbb{E}[||\hat{{\bm{m}}}_{n}-{\bm{m}}_{o}||])=O_{p}(||b_{k}||^{2}/\sqrt{p})$.

∎

###  S.2.6 Sufficient condition for Assumption 3.3

Let $\\{r_{n}:n\geq 1\\}$ denote a non-stochastic sequence that satisfies
$\displaystyle r_{n}\geq\max_{1\leq i\leq
n}||{\bm{m}}_{o}^{\prime}({\bm{x}}_{i}^{\prime}{\bm{x}}_{i}-\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]){\bm{m}}_{o}||/\sqrt{p}\
a.s.$ We impose the following sufficient condition for verifying Assumption
3.3.

######  Assumption S.2.1.

$r_{n}=o\big{(}\sqrt{\frac{n}{p\log p}}\big{)}$.

Assumption S.2.1 controls the rate at which $\displaystyle\max_{1\leq i\leq
n}||{\bm{m}}_{o}^{\prime}({\bm{x}}_{i}^{\prime}{\bm{x}}_{i}-\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]){\bm{m}}_{o}||/\sqrt{p}$
is allowed to diverge. Put together, Assumptions S.2.1 and 3.2 allow
${\bm{x}}$ to have an unbounded support. Assumptions S.2.1 and 3.2 are
important sufficient conditions for the matrix weak law of large numbers of
[29] – see also [12, Corollary 4.1 ]. The following result states and verifies
sufficient conditions for Assumption 3.3.

######  Lemma S.2.6.

Assumption 3.3 holds under Assumptions 3.2, 2.3, 3.1 and S.2.1.

#### Proof:

Under Assumptions 2.3 and 3.1,
$\varsigma_{n}^{2}:=||\sum_{i=1}^{n}\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]||=n||\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}]||\lesssim
n$. Let $R_{n}:=r_{n}\sqrt{p}$, then by Assumptions S.2.1 and 3.2,

$\displaystyle\frac{R_{n}\sqrt{\log k}}{\varsigma_{n}}$
$\displaystyle\leq\frac{R_{n}\sqrt{\log p}}{\varsigma_{n}}\lesssim
R_{n}\sqrt{\frac{\log p}{n}}=r_{n}\sqrt{\frac{p\log p}{n}}=o(1).$

The foregoing and Corollary 4.1 of [12] imply

$||\hat{A}_{m_{o}}-A_{m_{o}}||=||{\bm{m}}_{o}^{\prime}(\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]-\mathbb{E}[{\bm{x}}_{i}^{\prime}{\bm{x}}_{i}]){\bm{m}}_{o}||=O_{p}\Big{(}\sqrt{\frac{\log
p}{n}}\Big{)}.$

The conclusion then follows from Assumption 3.2.

∎

###  S.2.7 Calibration of $\gamma$

The crux of this exercise is to characterise the limiting distribution of the
first term in (S.16) under plausible conditions from which
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})$ is asymptotically
pivotal when $||\mathcal{B}_{k}||=o_{p}(n^{-1/2})$. $\gamma$ can then be set
to a suitable quantile from the limiting distribution.

Define the matrix
$\mathcal{C}_{n}^{k}:=\big{(}{\bm{m}}_{o}^{(k+1)}({\bm{m}}_{o}^{(k+1)\prime}{\bf
X}^{\prime}{\bf
X}{\bm{m}}_{o}^{(k+1)})^{-1}{\bm{m}}_{o}^{(k+1)\prime}-{\bm{m}}_{o}^{k}({\bm{m}}_{o}^{k\prime}{\bf
X}^{\prime}{\bf X}{\bm{m}}_{o}^{k})^{-1}{\bm{m}}_{o}^{k\prime}\big{)}$ and
observe that ${\bf X}\mathcal{C}_{n}^{k}{\bf X}^{\prime}$ is the difference
between two projection matrices. Recall the rank of a projection matrix equals
its trace and the eigenvalues are either zero or one – e.g., [20, Theorem 3.3
]. The following two conditions are imposed.

######  Condition S.2.1 (Homoskedasticity).

$\mathbb{E}[\varepsilon^{2}|{\bm{x}}]=\sigma^{2}$.

######  Condition S.2.2 (Rank $\mathcal{C}_{n}^{k}$).

$\mathrm{rank}[\mathcal{C}_{n}^{k}]=\mathrm{tr}[{\bf X}\mathcal{C}_{n}^{k}{\bf
X}^{\prime}]$ uniformly in $n,k$.

######  Lemma S.2.7.

Let $a_{k}<2$, i.e., $||\mathcal{B}_{k}||=o_{p}(n^{-1/2})$, Assumptions 3.1,
2.2, 2.3 and 3.3 hold. (a) Under Condition S.2.1,
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})\xrightarrow{d}\chi_{k+1}^{2}-\chi_{k}^{2}$
where $\chi_{k+1}^{2}$ and $\chi_{k}^{2}$ are mutually dependent chi-squared
distributed random variables with degrees of freedom $k+1$ and $k$,
respectively. (b) Further, suppose Condition S.2.2 holds, then
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})\xrightarrow{d}F(1,\infty)$.

The $F(1,\infty)$ limiting distribution of
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})$ under
$||\mathcal{B}_{k}||=o_{p}(n^{-1/2})$ does not depend on $k$. We set $\gamma$
to the $(1-\alpha)$’th quantile of the $F(1,\infty)$ distribution. For a
significance level $\alpha=0.1$, $\gamma=2.7$ (rounded to the first decimal
place).

#### Proof:

First,
$\mathrm{var}[\sqrt{n}{\bm{m}}_{o}^{k\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]]=\sigma^{2}A_{m_{o}^{k}}$
under Condition S.2.1. Thus, under $a_{k}<2$, i.e.,
$||\mathcal{B}_{k}||=o_{p}(n^{-1/2})$, Assumptions 3.1, 2.2, 2.3 and 3.3, the
Weak Law of Large Numbers, Lemma S.2.4, and Condition S.2.1,

$\displaystyle n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})$
$\displaystyle=n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{(k+1)}\hat{A}_{m_{o}^{k+1}}^{-1}{\bm{m}}_{o}^{(k+1)\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]/\mathbb{E}_{n}[\varepsilon_{i}^{2}]$
$\displaystyle-n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{k}\hat{A}_{m_{o}^{k}}^{-1}{\bm{m}}_{o}^{k\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]/\mathbb{E}_{n}[\varepsilon_{i}^{2}]+o_{p}(1)$
$\displaystyle=n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{(k+1)}(\sigma^{2}A_{m_{o}^{k+1}})^{-1}{\bm{m}}_{o}^{(k+1)\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]$
$\displaystyle-n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]{\bm{m}}_{o}^{k}(\sigma^{2}A_{m_{o}^{k}})^{-1}{\bm{m}}_{o}^{k\prime}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+o_{p}(1)$
$\displaystyle\xrightarrow{d}\chi_{k+1}^{2}-\chi_{k}^{2}.$

Second, by the non-negativity (in matrix notation) of
$n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1}):=\frac{1}{\mathbb{E}_{n}[\varepsilon_{i}^{2}]}\bm{\varepsilon}^{\prime}{\bf
X}\mathcal{C}_{n}^{k}{\bf X}^{\prime}\bm{\varepsilon}+o_{p}(1)$, the symmetric
$p\times p$ matrix $\mathcal{C}_{n}^{k}$ is positive semi-definite (up to some
$o_{p}(1)$ term). As ${\bf X}{\bm{m}}_{o}^{k}({\bm{m}}_{o}^{k\prime}{\bf
X}^{\prime}{\bf X}{\bm{m}}_{o}^{k})^{-1}{\bm{m}}_{o}^{k\prime}{\bf
X}^{\prime}$ is a projection matrix,

$\mathrm{tr}[{\bf X}\mathcal{C}_{n}^{k}{\bf
X}^{\prime}]=\mathrm{tr}\big{[}{\bf
X}{\bm{m}}_{o}^{(k+1)}({\bm{m}}_{o}^{(k+1)\prime}{\bf X}^{\prime}{\bf
X}{\bm{m}}_{o}^{(k+1)})^{-1}{\bm{m}}_{o}^{(k+1)\prime}{\bf
X}^{\prime}\big{]}-\mathrm{tr}\big{[}{\bf
X}{\bm{m}}_{o}^{k}({\bm{m}}_{o}^{k\prime}{\bf X}^{\prime}{\bf
X}{\bm{m}}_{o}^{k})^{-1}{\bm{m}}_{o}^{k\prime}{\bf X}^{\prime}\big{]}=1.$

It therefore follows from Condition S.2.2 that
$\mathrm{rank}[\mathcal{C}^{k}]=1$ where
$||\mathcal{C}^{k}-\mathcal{C}_{n}^{k}||=o_{p}(1)$. Combining the foregoing
with Assumption 2.3 implies
$\mathbb{E}[\sigma^{2}{\bm{x}}^{\prime}{\bm{x}}]^{1/2}\mathcal{C}^{k}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\sigma^{2}]^{1/2}$
has only one non-zero eigenvalue which equals $1$ thus,

$\displaystyle n\varphi_{n}(\hat{\beta}_{n}^{k},\hat{\beta}_{n}^{k+1})$
$\displaystyle=n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]\mathcal{C}^{k}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+o_{p}(1)$
$\displaystyle=n\mathbb{E}_{n}[\varepsilon_{i}{\bm{x}}_{i}]\mathbb{E}[\sigma^{2}{\bm{x}}^{\prime}{\bm{x}}]^{-1/2}\mathbb{E}[\sigma^{2}{\bm{x}}^{\prime}{\bm{x}}]^{1/2}\mathcal{C}^{k}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\sigma^{2}]^{1/2}\mathbb{E}[{\bm{x}}^{\prime}{\bm{x}}\sigma^{2}]^{-1/2}\mathbb{E}_{n}[{\bm{x}}_{i}^{\prime}\varepsilon_{i}]+o_{p}(1)$
$\displaystyle\xrightarrow{d}\chi_{\mathrm{rank}(\mathcal{C}^{k})}^{2}=\chi_{1}^{2}\overset{d}{=}F(1,\infty)$

where $\overset{d}{=}$ denotes equality in distribution.

∎

##  Appendix S.3 Additional Results and Details

###  S.3.1 Starting schemes of Algorithm 1

Instead of employing random starting schemes, our experiments show a faster
and better convergence based on preliminary estimates. Consider a block size,
say $\lfloor\kappa p\rfloor<n$ where $\kappa\in[1/p,1)$. Randomly assign $p$
covariates to groups with at most $\lfloor\kappa p\rfloor$ covariates. Regress
$y$ on each group of covariates separately and collect together to form
$\hat{\beta}_{n}^{(0)}$ in Algorithm 1. The special case of the above
procedure which we implement throughout in this paper is to set $\kappa=1/p$.
In this case, there is only one starting value possible as there is only one
way to assign $p$ covariates to $p$ groups. On the other hand, whenever
$\kappa>1/p$, one is able to generate several starting values and then select
the GPE that yields the least sum of squared residuals.

###  S.3.2 $O(n^{-1/2})$ perturbation of sparsity and discrete heterogeneity

We introduce a perturbation term in $\beta_{o}$ in DGP III in order to
evaluate the performance of estimators, including the GPE, and the sensitivity
of the GPE selection rule to “local” violations of discrete heterogeneity.
This exercise also serves to check the sensitivity of post-LASSO estimators to
“local” violations of sparsity. To this end, DGP III is modified with
$\beta_{oj}=\mathrm{I}(j\leq 5)+F_{\chi_{1}^{2}}^{-1}(\tau_{j})/\sqrt{2n}$
where $F_{\chi_{1}^{2}}^{-1}(\tau)$ denotes the $\tau$’th quantile of the
$\chi_{1}^{2}$ distribution.

Table S.5: DGP III - $O(n^{-1/2})$ perturbation.

|  $5\%$ Rejection Rate
---|---
$n$ | GPE | pLASSO.BCCH | pLASSO.CV | OLS
100 | 0.070 | 1 | 0.164 | 0.000
150 | 0.057 | 1 | 0.114 | 0.007
200 | 0.057 | 1 | 0.111 | 0.022
250 | 0.052 | 1 | 0.120 | 0.020
300 | 0.051 | 1 | 0.129 | 0.023
350 | 0.048 | 1 | 0.113 | 0.033
400 | 0.061 | 1 | 0.115 | 0.042

$\bar{k}$
---
GPE | pLASSO.BCCH | pLASSO.CV | OLS
2 | 5 | 34 | 75
2 | 5 | 32 | 75
2 | 5 | 35 | 75
2 | 5 | 38 | 75
2 | 5 | 36 | 75
2 | 5 | 37 | 75
2 | 5 | 39 | 75

Table S.5 presents simulation results from this exercise. Only the GPE
controls size meaningfully across all sample sizes; the OLS begins to show
good performance at larger sample sizes $350$ and $400$. Quite importantly,
although the GPE rule does not select more than $2$ groups most often, in
which case $||b_{k}||=O(n^{-1/2})$, the approximation error – being of the
same order as estimation error – does not affect inference as shown by the
good size control of the GPE. The rejection rate of the pLASSO.BCCH, which
equals $1$ throughout all sample sizes considered, shows a very high
sensitivity to “local” $O(n^{-1/2})$ violations of sparsity. It is, however,
worth emphasising that the choice of tuning parameter used for the pLASSO.BCCH
is robust to the $O(n^{-1/2})$ perturbation as it selects the correct number
of relevant variables. In sum, the above exercise suggests the robustness of
the GPE to local $O(n^{-1/2})$ perturbations to sparsity, discrete
heterogeneity, or both.