Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,269 characters · 11 sections · 55 citation commands
Clustered Covariate Regression
Let $y$ be a scalar outcome and ${\bm{x}}\in\mathbb{R}^p$ be a vector of covariates. Consider the standard linear model
where the random noise $\varepsilon$ satisfies $\mathbb{E}[\varepsilon|{\bm{x}}] = 0$ almost surely ($a.s.$). To estimate $\bm{\beta}_o$ with a finite sample of size $n$, one typically uses Ordinary Least Squares (OLS) when $p < n$. However, when $p$ is large relative to $n$, OLS gives highly imprecise estimates when $p<n$ or is infeasible when $p > n$. Unfortunately, the problem of high dimensionality in model estimation is increasingly becoming common in empirical practice.
A popular approach to tackling the problem of high dimensionality is to assume sparsity where only a few covariates are relevant and determine them using a regularisation technique, e.g., LASSO tibshirani1996regression.\footnote{Sparsity-at-zero on the absolute values of parameters is the notion of sparsity considered in this paper for simplicity. While it is conceivable to view clustered structures as some form of sparsity on pairwise differences between parameters, e.g., \`a la ke2015homogeneity, this risks belabouring the terminology.} Several sparsity-dependent estimators have been proposed in the literature, e.g., the post-LASSO estimators of belloni2014inference,belloni2017program. The Dantzig selector candes2007dantzig,bickel2009simultaneous,chernozhukov2022biased constitutes another important category of sparsity-dependent estimators. The rationale behind the aforementioned methods is that the penalty promotes sparsity in the estimate. Thus, a single parameter value, typically zero, can be assigned to a large subset of the covariates while parameters of the remaining covariates are “freely" estimated. As li-muller-2021linear points out, (approximate) sparsity in social science applications is not a priori obvious. Moreover, sparsity is not invariant to linear reparametrisations of covariates li-muller-2021linear,giannone2021economic.\footnote{ For example, sparsity in the rotated space of principal components, e.g, the first principal component only, corresponds to density in the un-transformed space of covariates giannone2021economic.}
Clustering is another useful technique used to deal with high dimensionality. Discrete heterogeneity, i.e., finite number of support points in $\bm{\beta}_o$, and well-separated groups which are commonly assumed in clustering-based estimators, improve upon the sparsity assumption, see, e.g., bonhomme2015grouped,ke2015homogeneity,sarafidis2012cross,su2016identifying,cheng2021clustering. While the aforementioned methods improve upon sparsity-dependent methods by achieving group-wise heterogeneity, they typically assume discrete heterogeneity of $\bm{\beta}_o$ (characterised by a finite number of latent-types) in panel data settings with a large time dimension.\footnote{To have a sense of the improvement, note that a clustering-based estimator is consistent under the discrete-non-sparse configuration with $\bm{\beta}_o=(1,\ldots,1)'$ using at least one group of covariates while sparsity-dependent estimators fail.} Notable exceptions include bonhomme2022discretizing which allows continuous heterogeneity and ke2015homogeneity which considers a cross-sectional setting but imposes discrete heterogeneity. Besides some game-theoretic settings, discrete heterogeneity does not appear to enjoy much support from, e.g., economic theory hahn2010panel.
The grouped parameter estimator (GPE) is fundamentally different from the aforementioned clustering-based estimators. (1) In this paper, each element of $\bm{\beta}_o$ corresponds to a different covariate. This contrasts with the literature on the above clustering-based estimators, which typically imposes clustering on a single parameter---such as an intercept or coefficients of particular covariates---that would otherwise be treated as homogeneous across units. (2) Our objective is to flexibly estimate a potentially high-dimensional regression function, in line with non-parametric approaches. In contrast to the literature on clustering-based estimators, which focuses on identifying discrete groupings and conducting inference on cluster-specific parameters, our goal is not to uncover latent classifications of units. (3) Clustering is a pure approximation tool in this paper as discrete heterogeneity is not assumed, groups need not be separated, and the resulting approximation error is controlled via a consistent choice of the number of groups. Thus, this paper is correctly viewed as complementing clustering-based estimators by not assuming discrete heterogeneity, not requiring well-separated groups, characterising the resulting approximation error, and choosing the number of groups such that the estimator is asymptotically normal in the presence of approximation error.
This paper bears some similarity with the recent paper of chernozhukov2023inference as both do not impose sparsity and treat the estimand as a product of two matrix-valued parameters. However, the focus of chernozhukov2023inference differs substantially from ours. (1) Their parameter is high dimensional $(n\times p)$ matrix-valued while $\bm{\beta}_o$ is $p\times 1$. Thus, we do not need to impose a rank condition on $\bm{\beta}_o$. The product representation in this paper comes from assigning $\bm{\beta}_o$ (up to an approximation bias) to $k\leq p$ groups. (2) It is not obvious how a standard textbook linear model such as (ref) with cross-sectional data can be cast in their proposed framework since their outcome, covariate matrix, parameter, and random noise is each $n\times p$ matrix-valued -- see chernozhukov2023inference. (3) They necessarily require $p \rightarrow \infty$ while this paper does not. (4) Unlike this paper, Assumption 4.3 of chernozhukov2023inference requires a bounded support of ${\bm{x}}$ with substantial restrictions on its variation -- see chernozhukov2023inference.
Since $\bm{\beta}_o$ is never observed in practice, wrongly assumed restrictions can result in misleading inference. For this reason, there is an emerging need in the literature for methods that are agnostic of and robust to different degrees of (non)-sparsity and heterogeneity of $\bm{\beta}_o$. The Grouped Parameter Estimator (GPE) introduced in this paper is developed using clustering as an approximation tool under the natural condition that the support of $\bm{\beta}_o$ be bounded. Relative to existing estimators, GPE is robust to pairwise configurations of (1) sparsity, approximate sparsity, or non-sparsity on the one hand and (2) continuous, mixed, or discrete heterogeneity of $\bm{\beta}_o$ on the other hand. Therefore, GPE eliminates the burden of assuming specific configurations of the unobservable $\bm{\beta}_o$, e.g., sparsity, and the attendant consequences for inference. A simple data-driven selection rule is proposed in this paper to determine the number of groups in order to control approximation error. Under mild conditions, the large sample properties of GPE are established under both fixed and increasing $p$ asymptotics.
The remainder of the paper is laid out as follows. (ref) presents the GPE estimand and estimator. The large sample properties of the grouped parameter estimator are established in (ref). (ref) examines the finite sample performance of GPE relative to competing estimators using simulation experiments. In (ref), we apply GPE in the analyses of price and income elasticities of demand for gasoline and conclude in (ref). The proofs of all theoretical results are collected in the appendix. An online appendix contains extra theoretical and simulation results.
The Grouped Parameter Estimator (GPE) exploits proximity between elements of $\bm{\beta}_o$ in order to decrease the size of the approximating model from $p$ to $k$ where $k\in \{\underline{k}_p,\ldots,\bar{k}_p\} $ and $1\leq \underline{k}_p\leq \bar{k}_p \leq p$. For instance, if $\beta_{oj}\approx \beta_{oj'}\approx \beta_{oj''}$ “sufficiently", one can assign a common parameter to the corresponding $(j,j',j'')$'th covariates in ${\bm{x}}$.\footnote{Technically, “sufficiently" as used here requires that the resulting approximation error not dominate the estimation error.} This implies that for $p\geq 2$, one can consider partitioning $\bm{\beta}_o$ into $k$ groups. Owing to the linear index ${\bm{x}}\bm{\beta}_o$ in (ref), grouping parameters translates into grouping covariates by design. Covariates in each group are aggregated and assigned their respective group parameters $\{\delta_{ol}, 1\leq l\leq k\} $. In essence, the idea behind GPE is to fit (ref) using $k$ instead of $p$ parameters. GPE simplifies to the OLS estimator when $k=p$.
Let ${\bm{m}}_o \in \mathcal{M}_k$ be a binary matrix that represents the optimal assignment (in the $k$-means clustering sense) of elements in $\bm{\beta}_o$ to $k$ groups, where $ \mathcal{M}_k \subset \mathbb{R}^{p\times k} $ denotes a space of $ p\times k $ binary matrices. The $(j,l)$'th entry of ${\bm{m}}_o$ equals $1$ if $\beta_{oj}$ belongs to group $l$ and zero otherwise. Since each row of ${\bm{m}}_o$ contains a single one and $k-1$ zeros, ${\bm{m}}_o'{\bm{m}}_o \in \mathbb{R}^{k\times k}$ is a diagonal matrix with group sizes on the diagonal and $\mathrm{tr}({\bm{m}}_o'{\bm{m}}_o) = p$, where $\mathrm{tr}(A)$ denotes the trace of matrix $A$. Also, because each element of $\bm{\beta}_o$ belongs to only one group and empty groups are not allowed, the inner product of any two columns of ${\bm{m}}_o$ is zero.
For a fixed $k \in \{\underline{k}_p,\bar{k}_p\} $, $\bm{\beta}_o$ is decomposed as
where $\bm{b}_k: = \bm{\beta}_o - {\bm{m}}_o\bm{\delta}_o $ is a non-stochastic approximation bias. The GPE estimand, viz. $ {\bm{m}}_o\bm{\delta}_o $ solves the following clustering problem:
where $\lVert \cdot \rVert$ denotes the Euclidean norm when applied to a vector and $\Delta \subset \mathbb{R}^k $. We introduce the following standard assumption in order to characterise the bound on the approximation bias $\bm{b}_k$.
(ref) is a standard assumption in the literature, which is used to bound the approximation bias $\bm{b}_k$ -- cf. su2016identifying, bonhomme2015grouped, cheng2021clustering, and bonhomme2022discretizing.
The GPE estimand ${\bm{m}}_o\bm{\delta}_o$ is unique since it is invariant to the relabelling/permutation of groups. Moreover, neither ${\bm{m}}_o$ nor $\bm{\delta}_o$ is of independent interest in this paper. We, nonetheless, follow the literature, e.g., bonhomme2015grouped,cheng2021clustering, in defining ${\bm{m}}\in\mathcal{M}_k$ up to permutations that leave ${\bm{m}}\bm{\delta}_o$ unchanged. To avoid overly dominant or empty groups that can induce ill-conditioning in the matrix $A_m:={\bm{m}}'\mathbb{E}[{\bm{x}}'{\bm{x}}]{\bm{m}}$ for all $ {\bm{m}} \in \mathcal{M}_k $, we maintain the following condition throughout the paper.
Group sizes induced by the group assignment matrix ${\bm{m}}\in \mathcal{M}_k$ constitute the diagonal elements of ${\bm{m}}'{\bm{m}}$. Since the columns of ${\bm{m}}$ are orthogonal, group sizes are also the eigenvalues of ${\bm{m}}'{\bm{m}}$. Thus, (ref) implies that $||{\bm{m}}|| $ is bounded for all ${\bm{m}}\in \mathcal{M}_k$ where $||\cdot||$ denotes the spectral norm when applied to matrices. In practice, (ref) can always be guaranteed by splitting up overly dominant groups and ensuring each group has at least one element. (ref) hence requires that the number of groups increases if $p$ is increasing. This also means that the lower and upper bounds on $k$ ($\underline{k}_p$ and $\bar{k}_p$), which are $p$-indexed, ought to increase under increasing-$p$ asymptotics. To alleviate concerns of scaling as covariates are aggregated into groups, we implicitly assume throughout the paper that $\mathbb{E}[{\bm{x}}]=0$.
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 $ and $a_n \asymp b_n$ means $a_n\lesssim b_n$ and $b_n\lesssim a_n$. Also, let $|\mathrm{supp}(\bm{\beta}_o)|$ denote the number of support points of $\bm{\beta}_o$. The following result derives the bound on $ ||\bm{b}_k|| $ in terms of $p$ and $k$.
(ref) shows $ ||\bm{b}_k||$ is decreasing in $k$ but increasing in $p$. Thus, the more granular the partitions of $\bm{\beta}_o$, the smaller the approximation error. (ref) is a univariate non-asymptotic and non-stochastic analogue of the stochastic bound on the $k$-means clustering approximation bias in graf2002rates -- cf. bonhomme2022discretizing.
Next, we motivate GPE from a regression point of view. Define the following least squares criterion:
Given ${\bm{m}}\in\mathcal{M}_k$, $\bm{\delta}_o({\bm{m}}):=({\bm{m}}'\mathbb{E}[{\bm{x}}'{\bm{x}}]{\bm{m}})^{-1}{\bm{m}}'\mathbb{E}[{\bm{x}}'y]$ is a deterministic function of ${\bm{m}}$ that minimises $ Q_o({\bm{m}}\bm{\delta}) $ with respect to $\bm{\delta}\in\Delta$ holding ${\bm{m}}\in\mathcal{M}_k$ fixed. Define $ \mathcal{S}_p := \{ \uptau \in \mathbb{R}^p: ||\uptau||=1 \} $ as the space of vectors with unit length. The following standard assumptions are imposed in order to characterise the identification of the GPE estimand ${\bm{m}}_o\bm{\delta}_o$.
(ref) is a standard exogeneity condition that precludes neither non-Gaussianity nor arbitrary heteroskedasticity of the model error $\varepsilon$. (ref) requires that the conditional variance be bounded, i.e., $ \mathbb{E}[\varepsilon^2|{\bm{x}}] \leq B \ a.s. $ The conditional moment restriction, $ \mathbb{E}[|\varepsilon|^d|{\bm{x}}] \leq B \ a.s. $, is fairly weak, cf. sarafidis2012cross. It is weaker than the sub-Gaussian tail condition imposed on $\varepsilon$ in ke2015homogeneity but slightly stronger than Condition SM(i) of belloni2012sparse. Let $\rho_{\mathrm{\max}}(A)$ denote the largest eigenvalue of the square matrix $A$.
In view of (ref), (ref) rules out perfect or very high collinearity in $ {\bm{x}} $ -- cf. belloni2015some and belloni2012sparse. It also ensures that $A_m:={\bm{m}}'\mathbb{E}[{\bm{x}}'{\bm{x}}]{\bm{m}}$ is positive-definite for all ${\bm{m}}\in\mathcal{M}_k$ and $k\in \{\underline{k}_p,\bar{k}_p\} $.
By the decomposition (ref), $Q_o({\bm{m}}\bm{\delta})$ in (ref) consists of terms involving $\bm{b}_k$ and terms not dependent on $\bm{b}_k$. Define the auxiliary criterion
where $ \sigma^2:= \mathbb{E}[\varepsilon^2] $. $\widetilde{Q}_o({\bm{m}}\bm{\delta})$ is uniquely minimised at $ {\bm{m}}\bm{\delta} = {\bm{m}}_o\bm{\delta}_o $. Thus, we can characterise the identification of the GPE estimand ${\bm{m}}_o\bm{\delta}_o $ up to the bias term $||\bm{b}_k||$.
By (ref), the identification of $ {\bm{m}}_o\bm{\delta}_o $ holds exactly if $||\bm{b}_k||=0$ or approximately if $||\bm{b}_k||/\sqrt{p}=o(1)$.
(ref) has two important takeaways: (1) the GPE estimand solves a (weighted) univariate $k$-means parameter clustering problem, and (2) the approximation biases from the weighted and unweighted clustering problems are of the same order. The latter is exploited in this paper to lower the computational cost.
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 mean. Define $ Q_n({\bm{m}}\bm{\delta}) :=\mathbb{E}_n[( y_i - {\bm{x}}_i{\bm{m}}\bm{\delta})^2]/p$ as the sample version of (ref). GPE is given by $ \widehat{\bm{\beta}}_n := \hat{{\bm{m}}}_n\hat{\bm{\delta}}_n $, where
Unlike in the population, a two-step estimation of $\hat{{\bm{m}}}_n$ and $\hat{\bm{\delta}}_n$ cannot be followed in finite samples because (1) $ \bm{\beta}_o $ is unknown; (2) $\hat{{\bm{m}}}_n$ and $\hat{\bm{\delta}}_n$ are interdependent and cannot be computed separately; and (3) the computational cost of an exhaustive search is prohibitive, especially for large $p$. In light of the foregoing, we propose the following iterative algorithm for GPE.
$k\leq (n-2)$ ensures there is at least one degree of freedom for estimation. $\hat{\beta}_{n,j}^{(s)} $ in Step (ref) is the OLS slope estimate of $ y-{\bm{x}}_{-j}\widehat{\bm{\beta}}_{n,-j}^{(s-1)} $ regressed on the $j$'th covariate $ x_j $ where $\widehat{\bm{\beta}}_{n,-j}^{(s-1)}:=(\hat{\beta}_{n,1}^{(s)},\ldots,\hat{\beta}_{n,j-1}^{(s)},\hat{\beta}_{n,j+1}^{(s-1)},\ldots,\hat{\beta}_{n,p}^{(s-1)})'$ and ${\bm{x}}_{-j}'$ is a $(p-1)\times 1$ vector formed from ${\bm{x}}$ by removing the $j$'th covariate. Updating the $ j $'th row of $ \hat{{\bm{m}}}^{(s)} $ involves assigning $ \hat{\beta}_{n,j}^{(s)} $ to the nearest element (in absolute value) in $ \hat{\bm{\delta}}_n^{(s-1)} $, and updating the group mean. This is the step that differs from, e.g., bonhomme2015grouped; it is a faster yet (asymptotically) equivalent step thanks to (ref)(b). Assigning $ \hat{\beta}_{n,j}^{(s)} $ to the nearest group in Step (ref) achieves the same goal up to the same order of approximation bias as obtaining the group assignment of $x_j$ which minimises the objective function. The latter involves computing the objective function $p-1$ times while the former only requires computing $k-1$ scalar absolute differences.
It is well known in the literature that the solution to iterative algorithms like (ref) is sensitive to starting values. Hence, it is essential to have reliable starting schemes. Our starting values are obtained from subgroup analyses based on the regression model as proposed by ma2017concave, except that in our case, we apply the concave penalty function to the pairwise differences of the model parameters instead of the intercepts. We proceed to solve the constrained optimisation problem using a modified alternating direction method of multipliers (ADMM) algorithm of boyd2011distributed -- see (ref) for details. Following the literature, e.g., bonhomme2022discretizing, our asymptotic theory focuses on the global minimum (indexed by $k$) while abstracting away from optimisation error.
We consider a sequence of models indexed by $ n $ with observed data $ \{({\bm{x}}_i,y_i), 1\leq i\leq n \} $, $ y_i = y_{i,n} $, $ {\bm{x}}_i = {\bm{x}}_{i,n} $, and $p=p_n$. The large sample properties of GPE are developed under the following additional assumptions.
The dependence of the sampling process on $ n $ in (ref) allows $ p $ to grow with $ n $ while allowing for fixed-$p$ asymptotics as well. (ref) characterises the rate at which $p$ is allowed to grow with $n$. As, (ref) is an asymptotic rate condition, it does not rule out $p>n$ in finite samples.
By (ref), the slower growth rate of $p$ in $n$ allowed by GPE, unlike (post)-LASSO estimators, is the price to pay for not explicitly exploiting sparsity (if it holds) which eliminates irrelevant covariates in a first step.
The following provides a characterisation of the convergence rate of $ ||\hat{{\bm{m}}}_n - {\bm{m}}_o|| $.
From (ref) above, one observes that the convergence of $\hat{{\bm{m}}}_n$ to ${\bm{m}}_o$ depends on $||\bm{b}_k||$, $n$, and $p$.
For a given $ {\bm{m}} \in \mathcal{M}_k $, the corresponding estimator is
The closed-form expression of $\hat{\bm{\delta}}_n({\bm{m}})$ is useful in concentrating out $\bm{\delta}$ from the objective function $Q_n({\bm{m}}\bm{\delta})$. Define $ \hat{A}_m:= {\bm{m}}'\mathbb{E}_n[{\bm{x}}_i'{\bm{x}}_i]{\bm{m}} $ as the sample analogue of $ A_m$. A consistency condition is imposed on $ \hat{A}_{m_o} $ in the following assumption.
(ref) is a high-level assumption -- cf. tropp2012user, chen2015optimal, and belloni2015some. (ref) in the Online Appendix verifies (ref) under sufficient conditions that allow an unbounded support of ${\bm{x}}$ and possibly increasing $p$. Also, observe that for a given $k$, (ref) is imposed on ${\bm{m}}_o$ and not on the entirety of $\mathcal{M}_k$.
Using the decomposition of $\widehat{\bm{\beta}}_n - \bm{\beta}_o$ into the estimation error $ \widehat{\mathcal{A}}_k $ and the approximation error $\widehat{\mathcal{B}}_k$, namely $\widehat{\bm{\beta}}_n - \bm{\beta}_o = \widehat{\mathcal{A}}_k + \widehat{\mathcal{B}}_k $ where
we characterise the convergence rate of GPE $\widehat{\bm{\beta}}_n$ and its asymptotic normality for any $ \uptau_p \in \mathcal{S}_p $. Define $\sigma_\theta: = (\uptau_p'{\bm{m}}_oA_{m_o}^{-1}{\bm{m}}_o'\mathbb{E}[{\bm{x}}'{\bm{x}}\varepsilon^2]{\bm{m}}_oA_{m_o}^{-1}{\bm{m}}_o'\uptau_p)^{1/2}$ where the dependence of $\sigma_\theta$ on $\uptau_p\in\mathcal{S}_p$ is suppressed for notational ease.
The underlying principle of GPE is that (ref) be estimable with fewer than $p$ parameters and thus conserve degrees of freedom. Although (ref) could be used to specify a selection rule of $k$ that ensures that the approximation error is controlled, such a choice rule would lack flexibility and be poorly adapted to the configuration of $\bm{\beta}_o$. For instance, under discrete or mixed heterogeneity of $\bm{\beta}_o$, a small $k$ should suffice to control approximation error. Further, simulations under a continuous-non-sparse configuration of $\bm{\beta}_o$ in (ref) suggest that a small number of groups $k$ suffices to control approximation error. Therefore, instead of using a large value of $k$ to drive approximation bias almost to zero at the cost of imprecise estimates and poor inference, we choose the smallest $k$ in a data-driven way to control approximation error -- see bonhomme2022discretizing for a selection rule based on a similar concept.
We select the number of groups per the rule $$ \hat{k}_n:= \min_{k\geq 1} \{ k: n\widehat{\varphi}_n(k)\leq C \} $$ where $$ \widehat{\varphi}_n(k):= \frac{\mathbb{E}_n[(y_i - {\bm{x}}_i\widehat{\bm{\beta}}_n^k)^2] - \mathbb{E}_n[(y_i - {\bm{x}}_i\widehat{\bm{\beta}}_n^{k+1})^2]}{\mathbb{E}_n[(y_i - {\bm{x}}_i\widehat{\bm{\beta}}_n^{k+1})^2]}, $$ $\widehat{\bm{\beta}}_n^k$ is GPE with $k$ groups (for notational emphasis), and $ C > 0 $ is a user-defined constant.\footnote{We recommend using $C=2.7$ -- see (ref) for details.} As $\mathbb{E}_n[(y_i - {\bm{x}}_i\widehat{\bm{\beta}}_n^k)^2]$ is non-increasing in $k$, $\widehat{\varphi}_n(k)$ is non-negative.\footnote{To see why, it suffices to create a $(k+1)$'th singleton group using the covariate whose parameter estimate generates the largest within-group absolute deviation.} Also, $\hat{k}_n$ is scale-invariant since $\widehat{\varphi}_n(k)$ is a ratio. The above selection rule bears an interesting similarity to model selection criteria such as the Bayesian Information Criterion (BIC) -- see, e.g., sarafidis2015partially. Unlike the BIC which has a penalty term to avoid over-fitting, setting $C$ to a constant and choosing $\hat{k}_n$ as the minimum admissible $k$ guards against over-fitting in our case. The following result shows that this simple selection rule is consistent, i.e., it selects $k$ such that the approximation error $\widehat{\mathcal{B}}_k$ in (ref) does not dominate the estimation error $\widehat{\mathcal{A}}_k$. This is similar in spirit to the choice of bandwidth in non-parametric methods. In view of (ref) and (ref)(a), let $\{a_k: k\geq 1\}$ be a decreasing sequence of positive numbers such that $||\widehat{\mathcal{B}}_k|| = O_p(n^{-1/a_k}) + o_p(n^{-1/2})$.
An important implication of (ref) is that $ n\widehat{\varphi}_n(k) = O_p(1) $ implies $||\widehat{\mathcal{B}}_k|| = O_p(n^{-1/2})$. $k=\hat{k}_n$ is chosen to ensure inference is not contaminated by the approximation error. Thus, $C > 0$ is set to a suitable constant under $a_k<2$, so that the approximation error is dominated by the estimation error, i.e., $||\widehat{\mathcal{B}}_k|| = o_p(n^{-1/2})$.\footnote{$C=2.7$ is the $90$’th percentile of the $\chi_1^2$ limiting distribution of $n\widehat{\varphi}_n(k)$ under two technical conditions and $||\widehat{\mathcal{B}}_k|| = o_p(n^{-1/2})$ -- see (ref).}
It is instructive to draw parallels between $k$ and its approximate analogue in the post-LASSO estimator of belloni2014inference (pLASSO hereafter), namely the sparsity index $s$, which is the bound on the number of covariates in a sparse model whose corresponding elements in $\bm{\beta}_o$ are non-zero.
The choice of $k$ is user-determined in a data-driven way whereas the pLASSO rate condition on $s$ is an assumption imposed on the unobservable $\bm{\beta}_o$. (ref) highlight the complementarity and trade-off between sparsity-dependent estimators and GPE in handling high covariate dimensionality. The consistent choice $k_n$ controls the approximation error and avoids unverifiable assumptions on the approximation error. In contrast, panel-data clustering-based methods, e.g., bonhomme2015grouped,ke2015homogeneity,sarafidis2015partially assume zero approximation error with well-separated groups subject to correctly choosing $k$, and regularisation-based methods such as belloni2014inference,chernozhukov2023inference impose order conditions on the approximation error.
This section focuses on the finite sample performance of GPE using simulated data. We compare the bias and rejection rates of competing estimators under Continuous-non-Sparse (CnS), Continuous-approximately-Sparse (CaS), and Discrete-Sparse (D-S) configurations of $ \bm{\beta}_o $.\footnote{See (ref) of the Online Appendix for details on the possible support configurations of $ \bm{\beta}_o $.} Estimators considered include (1) the proposed GPE; (2) the post-LASSO estimator of belloni2014inference (pLASSO); (3) post-Generalised Dantzig Selector (GDS) estimator of chernozhukov2022biased (pGDS); (4) OLS; (5) an “oracle" estimator which is OLS on sample size $3n$ (Orac.OLS); and an infeasible GPE that uses $\bm{\beta}_o$ as the starting value in (ref) (Orac.GPE).\footnote{pLASSO and pGDS estimators are based on tuning parameters recommended in belloni2014inference,chernozhukov2022biased, respectively.}$^,$\footnote{Simulations in an earlier draft suggest that pLASSO and pGDS with the tuning parameter selected via 10-fold cross-validation do not provide meaningful inference in any setting considered.} The infeasible Orac.OLS provides a useful benchmark for examining the gains in bias reduction and inference using a feasible estimator relative to OLS if the sample size were thrice as large.\footnote{In the context of experimental data, a practical consideration is whether a researcher would use, e.g., GPE or pLASSO to estimate a high-dimensional model or incur the added cost of collecting twice as much data in order to apply OLS.} The inclusion of Orac.GPE is helpful in gauging the performance of (ref) and the proposed selection rule of $\hat{k}_n$. In comparing it to GPE, one is also able to examine the reliability of the starting scheme.
$ \theta_o:= p^{-1/2}\sum_{j=1}^{p}\beta_{oj} $ is the scalar-valued parameter of interest in this simulation exercise. The comparison of competing estimators is based on the following metrics: (1) mean bias (MnB) namely $ ||\widehat{\bm{\beta}}_n - \bm{\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.) of a 5%-level $t$-test $ \mathbb{H}_o: \theta - \theta_o = 0 $; and (5) the median size of the approximating model across simulated samples ($ \mathrm{med}(\hat{k}_n) $). In computing the standard errors of $\hat{\theta}_n$ and rejection rates of pLASSO and pGDS, parameters shrunk to zero in the first step are treated as constants.
Define $ \Phi^{-1}(\tau_j) $ as the $ \tau_j $'th quantile of the standard normal distribution where $ \tau_j = 0.9(j-1)/(p-1) + 0.05 $. (ref) is the data-generating process with $ \varepsilon=U\sqrt{(1+x_1^2)/2}$ and the following specifications of $\bm{\beta}_o$:
The vector of covariates is generated as ${\bm{x}} \sim \mathcal{N}(\bm 0, \bm \Sigma)$, with $\bm \Sigma_{jj'} = 0.5^{|j-j'|}$. $ U \sim (\chi_1^2-1)/\sqrt{2} $ in DGPs CnS and $\mathrm{CaS_1}$ while $ U \sim \mathcal{N}(0,1) $ in DGPs $\mathrm{CaS_2}$ and D-S$_1$. The configuration in DGP CnS is dense; heterogeneity is continuous and the support is bounded away from zero. DGP $\mathrm{CaS_1}$ introduces approximate sparsity by allowing zero to fall within the bounds of $\mathrm{supp}(\bm{\beta}_o)$. DGPs $\mathrm{CaS_2}$ and D-S$_1$ are adapted from belloni2012sparse in order to compare GPE to pLASSO under approximately sparse and sparse configurations, respectively. Save DGP D-S$_1$ where heterogeneity is discrete, all other DGPs have continuous heterogeneity.\footnote{See (ref) for simulation results on other discrete and mixed heterogeneity configurations.} Heteroskedasticity is imposed in all DGPs at pairs $ (n,p) \in\{100,400\} \times \{75,150\} $. As OLS is not feasible when $ p>n $, OLS results are not reported for $(n,p)=(100,150)$. All results are based on 1000 simulated samples.
(ref) and (ref) present the findings of the simulation exercise. Fixing $ p $, one observes a general decrease in MnB, MAD, and RMSE for all estimators as $ n $ increases in (ref). Unsurprisingly, one observes a slight deterioration in performance as $ p $ is doubled and $ n $ is held fixed. Rejection rates are useful in showing which estimators provide reliable inference -- this is crucial in model estimation and hypothesis tests. Out of the four feasible estimators, only GPE delivers meaningful inference over all $ (n,p) $ configurations in DGPs Cns, CaS$_1$, and CaS$_2$. Orac.OLS, in spite of being run on samples of size $3n$, does not have good rejection rates across all $(n,p)$ configurations -- the Orac.OLS severely under-rejects in (ref) at $(n,p)=(100,150)$.
Even under approximate sparsity ((ref)), neither of the sparsity-dependent methods (pLASSO and pGDS), delivers reliable inference. GPE, in contrast, continues to perform reasonably. It is under exact sparsity ((ref)) that inference based on pLASSO and pGDS becomes reliable. This appears to suggest that the performance of sparsity-dependent estimators can be very sensitive to almost negligible deviations from sparsity.\footnote{A similar observation emerges in (ref) of the Online Appendix where sparsity is perturbed with a $O(n^{-1/2})$ term.} As sparsity is not easily verifiable in typical empirical applications, this constitutes a major drawback of pLASSO and pGDS for model estimation and inference. It ought to be borne in mind, however, that pLASSO and pGDS, among other sparsity-dependent methods, provide good approximations to the optimal instruments in the first stage of instrumental variable methods -- see belloni2012sparse,farrell2015robust,belloni2017program,hansen2014instrumental,carrasco2015regularized or the effect of a few covariates of interest in approximately sparse high-dimensional settings, e.g., belloni2014inference,chernozhukov2022biased.
Interestingly, $ \mathrm{med}(\hat{k}_n) = 2$ for GPE in DGP D-S$_1$ where $\bm{\beta}_o$ has exactly $2$ support points -- see (ref) and the bottom-right box-plot in (ref). In the other DGPs where $\bm{\beta}_o$ has continuous heterogeneity, one observes that $\mathrm{med}(\hat{k}_n)$ increases weakly in $n$ but not in $p$. Turning to the box plots in (ref), one observes that GPE and OLS' $\hat{\theta}_n-\theta_o$ appear well centred around zero across all four DGPs while that of pLASSO and pGDS are well-centred only in DGP D-S$_1$. This property of OLS does not translate into good size control; OLS severely under-rejects at $(n,p)=(100,75)$ and $(n,p)=(400,150)$ in all DGPs. Although tripling sample size and using OLS reduces bias relative to GPE, rejection rates of Orac.OLS can sometimes be very low, e.g., (ref) at $(n,p)=(100,150)$. The second rows of box plots in (ref) suggest GPE, relative to OLS, better estimates $\bm{\beta}_o$ (judging by the $||\widehat{\bm{\beta}}_n - \bm{\beta}_o||/\sqrt{p}$) as the configuration gets more sparse. With respect to the sparsity-dependent pLASSO and pGDS, a clear pattern does not emerge. The fairly small number of groups, see (ref), used by GPE while having a comparable or sometimes better performance than competing estimators suggests that its parsimony is not achieved at the expense of substantial bias or unreliable inference.
(ref) demonstrate that the good size control of GPE in (ref) is not at the expense of power. GPE is the only feasible estimator considered in this section that controls size meaningfully and has non-trivial power across all four DGPs. In (ref) for example, OLS has rejection rates equal to zero, and it has power under the alternative. pLASSO and pGDS control size meaningfully in (ref) and possess non-trivial power under the alternative. It is interesting to note that GPE in (ref) remains competitive in terms of size control and power under the alternative although the setting is favourable to pLASSO and pGDS.
In sum, this simulation exercise shows very robust performance in terms of bias reduction and good size control of GPE with its data-driven selection rule under (1) different empirically relevant but unknowable configurations of $\bm{\beta}_o$; (2) heteroskedastic $\varepsilon$; (3) Gaussian and non-Gaussian $\varepsilon$; and (4) skewed and symmetrically distributed $\varepsilon$. Robustness to all the aforementioned features is an important property of an estimator especially in the high dimensional setting as the configuration of $\bm{\beta}_o$ is typically unbeknownst to the practitioner or difficult to ascertain via economic theory. The exercise also shows a very high sensitivity of the sparsity-dependent pLASSO and pGDS to almost negligible deviations from sparsity.
In this section, we estimate price and income elasticities of demand for gasoline following yatchew2001household,chernozhukov2022biased,semenova2021debiased. We use data from the 1994-1996 Canadian National Private Vehicle Use Survey.\footnote{The data set is available on Professor Yatchew's website: \url{https://economics.utoronto.ca/yatchew/}.} The outcome variable is $\log$ consumption and the covariates of interest are $\log$ price and $\log$ income. The sparsity assumption for the GPE is not required, unlike chernozhukov2022biased,semenova2021debiased, which rely on it in this empirical setting.
The approach taken in constructing the covariate vector ${\bm{x}}$ follows semenova2021debiased,chernozhukov2022biased. The covariates include $\log $ price, $\log $ price squared, $\log $ income, $\log $ income squared, distance driven, distance driven squared, 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 ((ref)), and $\log$ income and its square in the second set of results ((ref)), respectively, resulting in $87$ covariates for each model. In this empirical exercise, we focus on price and income elasticities within sub-samples defined by age groups in (ref). This approach is useful in learning the variation in consumers' price and income elasticities of demand for gasoline by age group.
(ref), respectively, present the price and income elasticities of demand for gasoline by age group. One observes from (ref) that the price elasticity of demand is strictly negative and statistically different from zero at the 5% level for all age groups. Demand is, however, price inelastic for all age groups. The price elasticity of demand is non-monotone in age; consumers are most sensitive to price in the last two age groups where the price elasticity of demand does not vary (at the 5% level), unlike younger age groups that exhibit significant variation by age group in the price elasticity of demand.
As expected, the income elasticity of demand is non-negative. Save for the first age group where the income elasticity of demand is not statistically significant (point-wise) at the 5% level, the income elasticity of demand is point-wise significant for all other age groups. One cannot reject the hypothesis that income elasticity of demand equals a constant $c \in [0.029, 0.041] $ for all age groups. Moreover, the income elasticity of demand is not significantly different from zero for all age groups. One, however, rejects the hypothesis that consumers' demand for gasoline is not sensitive to income at all age groups. The $95\%$ simultaneous confidence bands suggest that gasoline is a necessity for all age groups.
Lastly, (ref) offers a comparison of price and income elasticities of demand using GPE and pLASSO. Although both GPE and pLASSO estimates are qualitatively similar, GPE estimates are more precise. It is also worth emphasising that $\log$ price and $\log$ income are not selected by the first stages of pLASSO in the full sample or any sub-sample; they are added in the second stage as the amelioration set -- see belloni2014inference.
As the specific configuration of a high-dimensional $\bm{\beta}_o$ cannot be known for sure in several empirical settings, there is a growing need for reliable, robust, and flexibly-adapted estimators. Although sparsity-dependent estimators such as pLASSO constitute the workhorse of econometric and applied work in high-dimensional settings and their theoretical properties are now well understood, their reliance on sparsity remains a major drawback. Existing clustering-based estimators are largely limited to fixed effects models with large $T$ panel data and do typically impose discrete heterogeneity. The aforementioned assumptions are restrictive and hard to verify or justify theoretically as they are imposed on the unobservable $\bm{\beta}_o$.
This paper makes a useful addition to the practitioner's toolkit of high-dimensional model estimation methods. The proposed GPE eliminates the burden on the practitioner of imposing typically unverifiable sparsity or discrete heterogeneity conditions on $\bm{\beta}_o$ while maintaining the natural requirement that $\bm{\beta}_o$ be defined on a bounded support. GPE is thus agnostic of and robust to different configurations of $\bm{\beta}_o$ -- it accommodates various degrees of (non)-sparsity and heterogeneity. Thanks to a simple data-driven selection rule, GPE achieves parsimony and adaptability to different empirically relevant configurations.
This paper develops and applies GPE to the linear model with exogenous covariates in an $i.i.d.$ setting. To enhance GPE's usefulness to the practitioner, some important extensions remain. It will be interesting to consider GPE under endogeneity with possibly many (weak) instruments, clustered or weakly dependent data, and non-linear models with both smooth and non-smooth objective functions.