Estimation and Inference in High-Dimensional Panel Data Models with Interactive Fixed Effects

\heading{Estimation and Inference in}{High-Dimensional Panel Data Models}{with Interactive Fixed Effects}

\authors {Maximilian R\"ucker\footnotemark[1]}{Ulm University} {Michael Vogt\footnotemark[2]}{Ulm University}

\authors {Oliver Linton\footnotemark[3]}{University of Cambridge} {Christopher Walsh\footnotemark[4]}{Newcastle University}

\footnotetext[1]{Address: Institute of Statistics, Department of Mathematics and Economics, Ulm University, Helmholtzstrasse 20, 89081 Ulm, Germany. Email: [email removed].} \footnotetext[2]{Corresponding author. Address: Institute of Statistics, Department of Mathematics and Economics, Ulm University, Helmholtzstrasse 20, 89081 Ulm, Germany. Email: [email removed].} \footnotetext[3]{Address: Faculty of Economics, University of Cambridge, Austin Robinson Building, Sidgwick Avenue, Cambridge, CB3 9DD, UK. Email: [email removed].} \footnotetext[4]{Address: Newcastle University Business School, 5 Barrack Road, Newcastle upon Tyne, NE1 4SE, UK. Email: [email removed].} \setcounter{footnote}{4}

Key words: panel data; interactive fixed effects; CCE estimator; high-dimensional model; lasso; desparsified lasso.

JEL classifications: C13; C23; C55.

Introduction

Nowadays, economic panel data sets are often “high-dimensional”: they contain a wide variety of time-varying characteristics and controls whose number is quite substantial in comparison to the sample size and may even exceed it. For example, in the case of low-frequency financial panel data, there is a rapidly evolving literature on the so-called “factor zoo”, which involves a large number of observed firm-specific characteristics that have been proposed as potential drivers of stock risk premia. Harvey-et-al2015 document $382$ such factors published in top journals and also point to the ongoing exponential growth in their number. A microeconomic example where the covariate dimensionality is pressing is presented in Belloni2016 who attempt to determine the social costs of gun ownership. Their analysis involves estimating the effect of gun prevalence on crime rates in a fixed-effect panel model. Their data set comprises information on $988$ explanatory variables, while the overall sample size is $nT=3705$ (with $n=195$ being the cross-section dimension and $T=19$ the time series length). Another example from macroeconomics concerns the determinants of economic growth. LuSu2016 estimate the effect of various possible determinants on the GDP growth rate in a panel with $n=108$ countries over $T=36$ years. As they point out, the survey by Durlauf2005 lists $145$ potential determinants of economic growth. Additionally allowing for interaction terms and nonlinear transformations of these variables, one easily arrives at a situation where the number of available covariates is comparable to or even exceeds the sample size.

Estimating high-dimensional panel models where the number of explanatory variables $p$ is large relative to the sample size $nT$ is far from trivial. Standard methods and theory from high-dimensional statistics are mostly tailored to a simple cross-sectional i.i.d.\ data structure. Panel data, in contrast, comprise complicated dependence structures that need to be taken into account: they usually exhibit non-negligible correlation in the time (and cross-sectional) dimension. Moreover, in order to account for unobserved heterogeneity, models with intricate error structures involving fixed-effect or factor-type components are considered rather than models with simple i.i.d.\ errors. These complicated data structures are presumably the reason why the literature on econometric methods for high-dimensional panels is quite limited. Below, we give a brief overview of the existing literature and of how our contribution fits into it.

The main aim of this paper is to develop estimation and inference methods for the high-dimensional panel data model with interactive fixed effects:

for $1\leq t\leq T$ and $1\leq i\leq n$, where $i$ is the cross-section index and $t$ the time series index, $Y_{it}$ is a real-valued response variable, $X_{it}$ is a vector of $p$ regressors and $\beta $ is the unknown parameter vector of length $p$. We allow the number of potential regressors $p$ to be very large but impose sparsity on $\beta$ in the sense that the number of non-zero elements of $\beta$ is relatively small compared to $p$. The error structure of the model comprises two components: a standard idiosyncratic error term $\varepsilon _{it}$ and the interactive fixed effects component $\gamma _{i}^{\top }F_{t}$, where $F_{t}$ is a vector of unobserved factors and $\gamma _{i}$ is a vector of unobserved factor loadings. The regressors $X_{it}$ are allowed to be correlated with the factor structure, which induces endogeneity issues in model (ref). The interactive fixed effects in (ref) allow to model unobserved heterogeneity in a quite flexible manner, in particular, much more flexibly than standard fixed effects $a_{i}$ and $b_{t}$ in a model of the form $Y_{it}=\beta^{\top}X_{it}+a_{i}+b_{t}+\varepsilon _{it}$.

In the traditional low-dimensional case where the number of regressors $p$ is small and fixed, model (ref) has been analyzed extensively in the literature. The most popular estimator of $\beta$ in this traditional setting is the common correlated effects (CCE) estimator of Pesaran2006. In this paper, we develop an estimator that can be regarded as an extension of the CCE method to the case where there are many potential covariates or controls. As the original CCE method, our approach is based on the following strategy: we “project away” the unobserved factors, i.e., we (approximately) eliminate them by applying a particular projection device to the response and the covariates. To estimate $\beta$, we then apply $\ell_1$-penalized least squares methods, i.e., lasso methods to the projected regression. We call our estimator a high-dimensional CCE (HD-CCE) estimator.

As detailed in Section (ref), the original CCE approach breaks down completely as soon as $p>T$ (and performs very poorly already for $p$ somewhat smaller than $T$), thus imposing very strong restrictions on the number of explanatory variables $p$. Our approach, in contrast, works for $p$ in a very wide range: we can deal with the standard \textquotedblleft low-dimensional\textquotedblright\ case where $p$ is small and fixed, the \textquotedblleft moderately high-dimensional\textquotedblright\ case where $p$ is fairly large but still smaller than the sample size $nT$ and the \textquotedblleft truly high-dimensional\textquotedblright\ case where $p$ exceeds the sample size $ nT$. Precise conditions on the size of $p$ in comparison to $n$ and $T$ are provided in the context of our theoretical results in Section (ref). For our estimator to work with $p$ in such a wide range, we require a projection device that approximately eliminates the factors no matter whether $p$ is small or big. To construct such a device, we make use of methods from high-dimensional factor analysis Fan2013 that are based on singular value decompositions of high-dimensional covariance matrices and principal components thresholding. Related principal components based methods have been used in the CCE context before Juodis2021b, however, for very different purposes and only in the low-dimensional case with $p$ small and fixed. Notably, the ability of our estimator to deal with both low- and high-dimensional situations does not come without cost: in contrast to the original CCE approach, we require an estimate of the number of factors. We propose a simple estimation procedure which can be regarded as a formalization of scree plots that are very common in applied factor analysis.

In the theoretical part of the paper, we derive the convergence rate of our HD-CCE estimator. We further establish an inference procedure for scalar parameters of interest. As usual in high-dimensional statistics, we need to desparsify or debias our lasso-type estimator in order to perform inference. We show that the desparsified version of our procedure is asymptotically normal and provide consistent standard errors that can be used for confidence intervals or hypothesis tests. We note that in contrast to most of the literature on panel models with interactive fixed effects, our methods and theory are not restricted to the large-$T$-case where both $n$ and $T$ tend to infinity but we also cover the small-$T$-case where $n$ tends to infinity and $T$ is a fixed natural number. This makes our methods applicable in a very wide range of application contexts. We provide numerical evidence on the performance of our methods by Monte Carlo experiments and illustrate the usefulness of our methods by an application to financial panel data.

Literature review

In the low-dimensional case with $p$ small and fixed, panel data models with interactive fixed effects are well understood. Since its introduction, the CCE estimator has become a standard tool for their analysis, giving rise to a whole new strand of the literature with numerous extensions such as Kapetanios2011, ChudikPesaranTosetti2011, PesaranTosetti2011 , ChudikPesaran2015, Westerlund2018, Westerlund2019, Juodis2021a, BrownSchmidtWooldridge2021 and Juodis2021b to name just a few. There are several alternatives to the CCE estimator which can be used to estimate the parameter vector $\beta$ in the low-dimensional case. The most important one is a least squares approach originally studied in Bai2009 and theoretically further explored in MoonWeidner2015 among others. The philosophy behind this approach is quite different from that of the CCE method: whereas the CCE approach eliminates the factors and the corresponding loadings by a suitable transformation of the model, the least squares approach treats them as additional parameters to be estimated. One disadvantage of the least squares approach is that the criterion function to be minimized is not convex. Hence, to compute the estimator, one needs to solve a non-convex optimization problem. Recently, least squares estimation with nuclear norm penalization has been proposed to overcome this problem. The resulting estimator minimizes a convex criterion function and can thus be efficiently computed by standard methods from convex optimization. It has, however, the disadvantage that its convergence rate is fairly slow in general. Recent studies on nuclear norm penalized estimators for panel data models with interactive fixed effects include ChernozhukovHansenLiaoZhu2018, BeyhumGautier2019 and MoonWeidner2019. A state-of-the-art review of methods for fixed effects panels, including interactive fixed effects and other variants, can be found in Bonhomme2024.

Whereas panel models with interactive fixed effects are well studied in the low-dimensional case, they are largely unexplored in high dimensions. Indeed, the literature on high-dimensional panels is rather limited in general. High-dimensional panel models with random and fixed effects have been considered in Kock2013, Kock2016, Belloni2016 and KockTang2019: Kock2013 derives theory for bridge estimators in both random and fixed effects models, while Kock2016 analyzes a model with a hybrid error structure that is in-between random and fixed effects. Belloni2016 introduce the so-called cluster-lasso to estimate the unknown parameters in a model with an individual fixed effect. Econometric methods for high-dimensional panel models with interactive fixed effects have been developed in LuSu2016 and BelloniChenPadillaWang2019: LuSu2016 extend the least squares method of Bai2009 to a high-dimensional dynamic panel model by adding a group-lasso penalty. However, they only consider a situation where $p$ grows fairly slowly with the sample size. BelloniChenPadillaWang2019 develop nuclear norm penalized estimation methods for high-dimensional quantile panel regression. A high-dimensional version of the panel partial factor model, which is closely related to panel models with interactive fixed effects, is investigated in HansenLiao2019. Cheng-et-al2024 study another closely related model framework, specifically, a high-dimensional panel regression model for financial data with interactive fixed effects where the factor loadings are driven nonparametrically by observed stock-specific characteristics or covariates. In their model, the covariates are assumed to be weakly dependent across both cross section and time series, which is incompatible with the type of factor structure in the covariates that we assume and exploit in this paper. To the best of our knowledge, CCE-type approaches suited to high dimensions have not been developed at all in the literature so far.

Structure of the paper

The model framework which underlies our theoretical analysis is introduced in detail in Section (ref), while identification issues are discussed in Section (ref). The HD-CCE estimator and its desparsified version are derived step by step in Section (ref). Section (ref) is dedicated to the practical implementation of our estimators, in particular, to the choice of the involved tuning parameters. The main theoretical results are laid out in Section (ref). We provide a simulation study in the supplementary material and illustrate our methods by an analysis of the “factor zoo” in Section (ref). A brief overview of the simulation study can be found in Section (ref).

R code

Our methods are implemented in the R package hdcce which can be downloaded from https://github.com/RueckerM/hdcce. Moreover, replication files are available at https://github.com/RueckerM/hdcce-ReplicationFiles.

Notation

Matrices are denoted by bold letters, whereas scalars and vectors are printed in normal font. For a vector $v = (v_1,\ldots,v_q)^\top \in \mathbb{R}^q$ and a set $S \subseteq \{1,\ldots,q\}$, we let $v_S = (v_i: i \in S)$ be the vector which consists of the entries $ v_i$ with $i \in S$ only. In addition, we sometimes write $v_{-i}$ to denote the vector $v$ without the $i$-th component. We let $\|v\| = (\sum_i v_i^2)^{1/2}$ denote the Euclidean norm of $v$, $\|v\|_1 = \sum_i|v_1|$ its $ \ell_1$-norm, and $\|v\|_{\infty} = \max_i|v_i|$ its $\ell_\infty$-norm. For a generic matrix $\boldsymbol{A}^{T \times p}$, we denote the row vectors by $A_t$ and the column vectors by $A_{(j)}$, that is, $\boldsymbol{A} = (A_1 \ldots A_T)^\top = (A_{(1)} \ldots A_{(p)})$. Moreover, the matrix $\boldsymbol{A}$ without the $t$-th row is denoted by $\boldsymbol{A}_{-t}$ and that without the $j$-th column by $\boldsymbol{A}_{(-j)}$. The symbols ${\color{black}{\psi}}_{\min}(\boldsymbol{A})$ and ${\color{black}{\psi}} _{\max}(\boldsymbol{A})$ are used to denote the minimal and maximal eigenvalue of a square matrix $\boldsymbol{A} \in \mathbb{R}^{q \times q}$. In addition, we sometimes write ${\color{black}{\psi}}_1(\boldsymbol{A}) \ge {\color{black}{\psi}}_2(\boldsymbol{A}) \ge \ldots \ge {\color{black}{\psi}}_q(\boldsymbol{A})$ to denote the eigenvalues of $\boldsymbol{A}$ (in decreasing order). For a general (not necessarily square) matrix $\boldsymbol{A} = (a_{ij})$, $\| \boldsymbol{A} \|$, $\|\boldsymbol{A} \|_1 $, $\|\boldsymbol{A}\|_\infty$ and $\|\boldsymbol{A}\|_{\text{max}}$ are its spectral norm, $\ell_1$-norm, $\ell_\infty$-norm and elementwise norm, respectively. In particular, $\| \boldsymbol{A} \| = {\color{black}{\psi}}_{\max}^{1/2}(\boldsymbol{A}^\top \boldsymbol{A})$, $\| \boldsymbol{A}\|_1 = \max_j \sum_i |a_{ij}|$, $\|\boldsymbol{A}\|_\infty = \max_i \sum_j |a_{ij}|$ and $\|\boldsymbol{A}\|_{\max} = \max_{ij} |a_{ij}|$. The symbol $\boldsymbol{A} ^{-}$ stands for the generalized inverse of a matrix $\boldsymbol{A}$ and the symbol $\boldsymbol{I}_q$ for the $q \times q$ identity matrix. Sometimes, we also write $ \boldsymbol{I}$ instead of $\boldsymbol{I}_q$ for short. Finally, the indicator function is denoted by $1(\cdot)$ and the cardinality of a set $S$ by $|S|$.

Model framework

We observe a sample of panel data $\{ (Y_{it}, X_{it}): 1 \le t \le T, \, 1 \le i \le n \}$ with real-valued random variables $Y_{it}$ and $\mathbb{R}^p$-valued random vectors $X_{it} = (X_{it,1},\ldots,X_{it,p})^\top$, where $n$ is the cross-section dimension and $T$ the time series length. We consider the following two scenarios:

We regard both $T$ and $p$ as a function of $n$, that is, $T=T(n)$ and $p=p(n)$. Hence, asymptotic statements are to be understood in the sense that $n \to \infty$ (and $T=T(n) \to \infty$ in the large-$T$-case). The dimension $p$ of the random vector $X_{it}$ is allowed to be large, potentially much larger than $n$ and $T$. Put differently, we allow $p$ to grow with $n$ (and $T$). The only restriction is that $p$ does not grow too quickly compared to $n$ (and $T$). Hence, the methods and theory of this paper are valid for any $p$ which is not too large compared to $n$ (and $T$). In particular, we cover both the traditional low-dimensional case where $p$ is small and fixed and the high-dimensional case where $p$ grows potentially much faster than $n$ (and $T$). Precise conditions on the size of $p$ compared to $n$ and $T$ are provided in Section (ref).

We consider a high-dimensional version of the linear panel data model with interactive fixed effects analyzed in Pesaran2006. The model has the form

for each cross-sectional unit $i$, where $Y_i = (Y_{i1},\ldots,Y_{iT})^\top \in \mathbb{R}^T$ is the response vector, $\beta = (\beta_1,\ldots,\beta_p)^\top$ is the unknown parameter vector, $\boldsymbol{X}_i = (X_{i1} \ldots X_{iT})^\top \in \mathbb{R}^{T \times p}$ is the regressor matrix, $\varepsilon_i = (\varepsilon_{i1},\ldots,\varepsilon_{iT})^\top \in \mathbb{R}^T$ is the vector of idiosyncratic errors with $\mathbb{E}[\varepsilon_{it}] = 0$ for all $i$ and $t$, and $\boldsymbol{F} \gamma_i$ is the interactive fixed effects part of the error. More specifically, $\boldsymbol{F} = (F_1 \ldots F_T)^\top \in \mathbb{R}^{T \times K}$ with $F_t = (F_{t,1},\ldots,F_{t,K})^\top$ is a matrix of unobserved factors and $\gamma_i = (\gamma_{i,1},\ldots,\gamma_{i,K})^\top$ is a vector of (unknown) individual-specific factor loadings. Throughout the paper, we treat the factors $F_t$ as non-random parameters. Put differently, we implicitly condition on the factors $F_1,\ldots,F_T$ in our theoretical analysis as is common in the literature MoonWeidner2015. The factor loadings $\gamma_i$, in contrast, are considered to be random. The regressors in (ref) are supposed to have the structure

where $\boldsymbol{\Gamma}_i \in \mathbb{R}^{p \times K}$ is a matrix of individual-specific factor loadings and $\boldsymbol{Z}_i = (Z_{i1} \ldots$ $\ldots Z_{iT})^\top \in \mathbb{R}^{T \times p}$ represents the idiosyncratic part of the regressors with $\mathbb{E}[Z_{it}] = 0$ for all $i$ and $t$. This structure implies that the regressors $\boldsymbol{X}_i$ are in general correlated with the unobserved part of equation (ref), $e_i = \boldsymbol{F}\gamma_i + \varepsilon_i$, via the interactive fixed effects.

The main difference of model (ref)--(ref) from Pesaran's original model is that we allow the dimension $p$ of the regressors $X_{it} = (X_{it,1},\ldots,X_{it,p})^\top$ to be large, possibly much larger than the overall sample size $nT$. Without structural constraints on the parameter vector $\beta$, model (ref)--(ref) is not estimable in general. As usual in the literature on high-dimensional statistics, we impose a sparsity constraint on $\beta$. In particular, we assume that the set $S=\{j: \beta_j \ne 0\}$ of non-zero components of $\beta$ has cardinality $s := |S|$ considerably smaller than the sample size $nT$. Hence, only a small subset of regressors is active, that is, enters the model with a non-zero coefficient. Precise conditions on the size of the sparsity index $s$ are provided in Section (ref). Notably, there have been some attempts to perform estimation and inference in high-dimensional models with non-sparse structures in recent years ZhuBradic2018,SilinFan2022. However, even though the assumption of sparsity is not harmless Kolesar2025, we here follow the main bulk of the literature on high-dimensional statistics and work under a sparsity constraint. In contrast to the number of regressors $p$, the number of unknown factors $K$ is assumed to be fixed in magnitude as in Pesaran's model. Assuming that $K$ is comparably small makes sense as $K$ plays a role analogous to the number of active regressors $s$ rather than the total number of regressors $p$. We do not assume that $K$ is known a priori, and determine it from the data.\footnote{Following Pesaran2006, one may additionally include observed factors in model (ref)--(ref) and allow for heterogeneous parameter vectors $\beta_i = \beta + \eta_i$ with i.i.d.\ disturbances $\eta_i$. In particular, as long as the random disturbances $\eta_i$ produce sparse parameter vectors $\beta_i$ and we are in the large-$T$-case, it is possible to estimate the individual $\beta_i$'s. However, if the $\beta_i$'s are non-sparse, it will in general only be possible to estimate the (sparse) mean vector $\beta$.}

If our interest focuses on point estimation of $\beta$, the above model description is fully sufficient. If the aim is to perform inference, in contrast, we need additional structure. Suppose in particular we want to compute confidence bands for the coefficient $\beta_j$ of the $j$-th regressor. To be able to do so, we additionally impose the nodewise regression equation

where $X_{i(j)}$ is the $j$-th column of the matrix $\boldsymbol{X}_i$, $\boldsymbol{X}_{i(-j)}$ is the matrix $\boldsymbol{X}_i$ without the $j$-th column, ${\color{black}{\theta}}$ is a sparse parameter vector, $u_i$ denotes the idiosyncratic error term with $\mathbb{E}[u_i] = 0$ and $\nu_i \in \mathbb{R}^K$ is a vector of factor loadings. The quantities ${\color{black}{\theta}}$, $\nu_i$ and $u_i$ depend on $j$. For convenience, however, we suppress this dependence in the notation. According to (ref), the $j$-th regressor can be represented as a sparse linear function of the other regressors plus an interactive fixed effects error structure. Such a nodewise regression equation is very common in high-dimensional inference, both when desparsified lasso techniques vandeGeer2014 and double selection techniques Belloni2014 are used. Notably, it is no problem to satisfy both the regressor equation (ref) and the nodewise equation (ref) for the $j$-th regressor in our framework. In particular, if the $j$-th regressor is modelled via the nodewise equation (ref), it also fulfills (ref): $X_{i(j)} = \boldsymbol{F} \Gamma_{i,j} + Z_{i(j)}$ with $\Gamma_{i,j} := \boldsymbol{\Gamma}_{i,-j}^\top \theta + \nu_i$ and $Z_{i(j)} := \boldsymbol{Z}_{i(-j)} \theta + u_i$, which follows immediately upon plugging equation (ref) for all but the $j$-th regressor, i.e., the equation $\boldsymbol{X}_{i(-j)} = \boldsymbol{F} \boldsymbol{\Gamma}_{i,-j}^\top + \boldsymbol{Z}_{i(-j)}$ into (ref).

A list of the technical conditions that we impose on the model components in equations (ref), (ref) and (ref) to derive our theoretical results can be found in Section (ref).

Identification

Model (ref)--(ref) contains the following unobserved components: the parameter vector $\beta$, the factor structure $\Theta_{\text{fac}} = \{ \boldsymbol{F}, \{ \boldsymbol{\Gamma}_i, \gamma_i \}_{i=1}^n \}$ consisting of the factors and their loadings, and the idiosyncratic structure $\Theta_{\text{idio}} =\{ \boldsymbol{Z}_i, \varepsilon_i\}_{i=1}^n$ consisting of the idiosyncratic part of the regressors and the idiosyncratic errors. Importantly, the parameter vector $\beta$, the factor structure $\Theta_{\text{fac}}$ and the idiosyncratic structure $\Theta_{\text{idio}}$ are in general not identified. Put differently, the parameter vector $\beta$, the factor structure $\Theta_{\text{fac}}$ and the idiosyncratic structure $\Theta_{\text{idio}}$ which satisfy model equations (ref)--(ref) and the technical assumptions of Section (ref) are in general not unique.

In what follows, we show that the parameter vector $\beta$ and the number of factors $K$ are identified if certain additional constraints are imposed. That the factor structure $\Theta_{\text{fac}}$ (apart from $K$) and the idiosyncratic structure $\Theta_{\text{idio}}$ remain unidentified is no problem at all for our methods and theory. For our theoretical arguments to work, it suffices to consider some factor structure $\Theta_{\text{fac}}$ and some idiosyncratic structure $\Theta_{\text{idio}}$ such that the model equations and the technical conditions are fulfilled. Which version is considered does not matter.

We start with identification of $K$. As usual, we normalize the factors $F_t$ to be orthonormal:

Given this normalization, we impose the following assumption on the mean loading matrix $\boldsymbol{\Gamma} = \mathbb{E}[\boldsymbol{\Gamma}_i] \in \mathbb{R}^{p \times K}$:

(ref) is a standard condition in the literature on high-dimensional approximate factor models; see e.g.\ Fan2013 and BaiLiao2016. By imposing it, we focus on the case of strong factors.\footnote{It is in principle possible to weaken (ref). In particular, our theory does not require all eigenvalues of $\boldsymbol{\Gamma}^\top \boldsymbol{\Gamma}$ to be of order $p$ as assumed in (ref). Instead, we could allow the eigenvalues to be of different order as long as their orders are large enough, in particular, larger than $p \sqrt{\log p} / \sqrt{n}$. Such a generalization of condition (ref) would however influence the convergences rates of our HD-CCE estimator.} Under (ref), the eigenvalues of $\boldsymbol{\Gamma}^\top \boldsymbol{\Gamma} / p$ are strictly positive for all $p$, which implies that the matrix $\boldsymbol{\Gamma}$ has full rank $K$ for all $p$. An analogous full-rank condition is required in the original CCE approach of Pesaran2006.\footnote{Note that Pesaran2006 also treats the rank-deficient case. However, as shown in WesterlundUrbain2013, his results only hold if $\gamma_i$ and $\boldsymbol{\Gamma}_i$ are uncorrelated. Hence, the full-rank condition on $\boldsymbol{\Gamma}$ is indeed required in the original CCE approach unless one is willing to make the very strong assumption that $\gamma_i$ and $\boldsymbol{\Gamma}_i$ are uncorrelated.} Under (ref) and (ref), we can prove the following identification result.

The proof of this as well as the subsequent results on identification can be found in the supplementary material.

We next turn to identification of $\beta$. In the high-dimensional case with $p$ potentially larger than the full sample size $nT$ itself, there is of course no way to identify $\beta$ in general. However, we can get identification if we restrict attention to para\-meter vectors $\beta$ with certain properties. Specifically, we focus on vectors $\beta$ which are $s$-sparse, that is, which have at most $s$ non-zero components. As we will see, under certain constraints, there is a unique $s$-sparse parameter vector $\beta$ which satisfies model (ref). In order to formulate the precise identification result, we introduce some notation: Let $\mathcal{L}_{\boldsymbol{F}} = \{\boldsymbol{F} v: v \in \mathbb{R}^K\}$ be the column space of the matrix $\boldsymbol{F}$ and $\boldsymbol{\Pi} = \boldsymbol{I} - \boldsymbol{F} (\boldsymbol{F}^\top \boldsymbol{F})^{-1} \boldsymbol{F}^\top$ the projection matrix onto the orthogonal complement of $\mathcal{L}_{\boldsymbol{F}}$. Since $\boldsymbol{\Pi} \boldsymbol{F} = \boldsymbol{0}$ by construction, applying $\boldsymbol{\Pi}$ to the model equation in (ref) yields $\boldsymbol{\Pi} Y_i = \boldsymbol{\Pi} \boldsymbol{X}_i \beta + \boldsymbol{\Pi} \varepsilon_i$. Stacking the projected model equations $\boldsymbol{\Pi} Y_i = \boldsymbol{\Pi} \boldsymbol{X}_i \beta + \boldsymbol{\Pi} \varepsilon_i$ for all $i$, we obtain the model

where \[ Y^\perp =

, \ \ \boldsymbol{X}^\perp =

, \ \ \varepsilon^\perp =

. \] In order to identify the $s$-sparse parameter vector $\beta$, we impose a restricted eigenvalue (or compatibility) condition on the design matrix $\boldsymbol{X}^\perp$ in model (ref). Such a condition is very common in high-dimensional statistics BuehlmannvandeGeer2011 and can be formulated as follows.

We assume that with probability tending to $1$, the design matrix $\boldsymbol{X}^\perp$ satisfies the $\textnormal{RE}(I,\varphi)$ condition for all $I \subseteq \{1,\ldots,p\}$ with $|I| \le 2s$. More formally:

Under (ref), the parameter vector $\beta$ is identified in the following sense.

How reasonable are the restricted eigenvalue conditions on $\boldsymbol{X}^\perp$ in (ref)? It can be shown that (ref) is implied by an analogous assumption on the idiosyncratic matrix $\boldsymbol{Z} = (\boldsymbol{Z}_1^\top \ldots \boldsymbol{Z}_n^\top)^\top$ from equation (ref). Specifically, Lemmas (ref) and (ref) in the supplementary material show that (ref) is implied by the following condition:

As the matrix $\boldsymbol{Z}$ does not depend on the factors $\boldsymbol{F}$, it has a completely standard structure and can be regarded as an “ordinary” design matrix in a setting with sample size $nT$ and dimension $p$. Hence, imposing a restricted eigenvalue condition on $\boldsymbol{Z}$ is as restrictive or unrestrictive as imposing such a condition on the design matrix in a plain vanilla high-dimensional linear model. Notably, it is possible to verify that $\boldsymbol{Z}$ fulfills (ref) under certain distributional assumptions. Theorem 1 in RaskuttiWainwrightYu2010, for example, shows that (ref) is satisfied if the random vectors $Z_{it}$ are independent across $i$ and $t$ and $Z_{it} \sim N(0,\boldsymbol{\Lambda})$ with ${\color{black}{\psi}}_{\min}(\boldsymbol{\Lambda}) \ge c > 0$ and $\max_{1 \le j \le p} \boldsymbol{\Lambda}_{jj} \le C < \infty$. This result remains to hold true when the variables $Z_{it}$ are non-Gaussian with sufficiently light tails; see e.g.\ Theorem 7 in JavanmardMontanari2014.

Estimation and inference

A very popular technique to estimate the parameter vector $\beta$ in the low-dimen\-sional case is the common correlated effects (CCE) approach of Pesaran2006. In the high-dimensional case, however, this estimation technique breaks down and straightforward extensions are not possible. In this section, we construct a novel estimator which does work in high dimensions. As it is similar in spirit to the CCE approach, we call it a high-dimensional CCE estimator, or HD-CCE estimator for short. The section is structured as follows: First, we outline the general strategy to estimate $\beta$ which underlies both our and the CCE approach. We then explain why the CCE estimator collapses in high dimensions. Next, we introduce our estimation approach and give some heuristic discussion why it works. Finally, we explain how to perform inference based on the HD-CCE estimator.

A general estimation strategy

A general strategy to estimate $\beta$ in the panel data model (ref)--(ref) with interactive fixed effects is to eliminate or “project away” the unknown factors from the model equation by a suitable transformation and then to apply regression techniques to the transformed data.

To formalize this idea, we first consider the oracle case where the factors $\boldsymbol{F}$ are observed. In this case, the model equation $Y_i = \boldsymbol{X}_i \beta + \boldsymbol{F} \gamma_i + \varepsilon_i$ can be regarded as a partitioned regression model, where the factors $\boldsymbol{F}$ are additional regressors and the design matrix is given by $(\boldsymbol{X}_i \ \boldsymbol{F})$. The factors can be eliminated as follows: As already defined above, let $\mathcal{L}_{\boldsymbol{F}} = \{ \boldsymbol{F} v: v \in \mathbb{R}^K \}$ be the column space of the factor matrix $\boldsymbol{F}$ and \[ \boldsymbol{\Pi} = \boldsymbol{I} - \boldsymbol{F} (\boldsymbol{F}^\top \boldsymbol{F})^{-1} \boldsymbol{F}^\top \] the projection matrix onto the orthogonal complement of $\mathcal{L}_{\boldsymbol{F}}$. Since $\boldsymbol{\Pi} \boldsymbol{F} = \boldsymbol{0}$ by construction, we can pre-multiply the model equation by $\boldsymbol{\Pi}$ to get that $\boldsymbol{\Pi} Y_i = \boldsymbol{\Pi} \boldsymbol{X}_i \beta + \boldsymbol{\Pi} \boldsymbol{F} \gamma_i + \boldsymbol{\Pi} \varepsilon_i = \boldsymbol{\Pi} \boldsymbol{X}_i \beta + \boldsymbol{\Pi} \varepsilon_i$, thus “projecting away” the factors $\boldsymbol{F}$. An estimator of $\beta$ can be obtained by applying regression techniques for high-dimensional linear models to the transformed data $\{ (\boldsymbol{\Pi} Y_i, \boldsymbol{\Pi} \boldsymbol{X}_i): 1 \le i \le n \}$. Specifically, running a lasso regression on the transformed data leads to the estimator \[ \widehat{\beta}_{\color{black}{\lambda}}^{\text{oracle}} \in \underset{b \in \mathbb{R}^p}{\text{argmin}} \bigg\{ \frac{1}{nT} \sum_{i=1}^n \big\| \boldsymbol{\Pi} Y_i - \boldsymbol{\Pi} \boldsymbol{X}_i b \big\|^2 + {\color{black}{\lambda}} \|b\|_1 \bigg\}, \] where ${\color{black}{\lambda}} > 0$ is the penalty constant of the lasso. If $p$ is much smaller than the sample size $nT$ (in particular, in the low-dimensional case with fixed $p$), there is of course no need to work with the lasso. One may rather set the penalty constant ${\color{black}{\lambda}}$ to $0$ and use the least squares estimator $\widehat{\beta}_0^{\text{oracle}}$.

Obviously, the oracle estimator $\widehat{\beta}_{\color{black}{\lambda}}^{\text{oracle}}$ is not feasible in practice: since the factors $\boldsymbol{F}$ are not observed, the projection matrix $\boldsymbol{\Pi} = \boldsymbol{I} - \boldsymbol{F} (\boldsymbol{F}^\top \boldsymbol{F})^{-1} \boldsymbol{F}^\top$ and thus the estimator $\widehat{\beta}_{\color{black}{\lambda}}^{\text{oracle}}$ cannot be computed. To obtain a feasible estimator of $\beta$, we need to replace the unknown matrix $\boldsymbol{\Pi}$ by a proxy. The construction of such a proxy in high dimensions turns out to be quite intricate. This is the main technical challenge we need to deal with.

Breakdown of the CCE estimator in high dimensions

Before we construct a proxy of $\boldsymbol{\Pi}$ in high dimensions, we review the traditional low-dimensional case where (i) the number of regressors $p$ is a fixed natural number, (ii) $p$ is small in the sense that $p < T$, and (iii) the number of factors $K$ is not larger than $p$, that is, $K \le p$.

The CCE approach of Pesaran2006 provides an elegant way to proxy $\boldsymbol{\Pi}$ in this low-dimensional case. For simplicity, we only use the regressors $X_{it}$ for the construction (and thus ignore the responses $Y_{it}$). This gives a clearer picture of the approach and does not affect our argumentation. For a generic random variable $R_{it}$, let $\overline{R}_t = n^{-1} \sum_{i=1}^n R_{it}$ be its cross-sectional average. The CCE approach proxies the projection matrix $\boldsymbol{\Pi} = \boldsymbol{I} - \boldsymbol{F} (\boldsymbol{F}^\top \boldsymbol{F})^{-1} \boldsymbol{F}^\top$ by \[ \overline{\boldsymbol{\Pi}} = \boldsymbol{I} - \overline{\boldsymbol{X}} (\overline{\boldsymbol{X}}^\top \overline{\boldsymbol{X}})^{-} \overline{\boldsymbol{X}}^\top, \] where $\overline{\boldsymbol{X}} = (\overline{X}_1 \ldots \overline{X}_T)^\top$ is the matrix containing the cross-sectional averages $\overline{X}_t = (\overline{X}_{t,1},\ldots,\overline{X}_{t,p})^\top$ of the regressor variables. Under suitable regularity conditions, it can be shown that $\overline{\boldsymbol{\Pi}} Y_i \approx \overline{\boldsymbol{\Pi}} \boldsymbol{X}_i \beta + \overline{\boldsymbol{\Pi}} \varepsilon_i$ in the low-dimensional case. Hence, pre-multiplying the model equation by $\overline{\boldsymbol{\Pi}}$ approximately eliminates the factors. We may thus use $\overline{\boldsymbol{\Pi}}$ as an observable proxy of $\boldsymbol{\Pi}$ and estimate $\beta$ by applying least squares methods to the sample of transformed data $\{ (\overline{\boldsymbol{\Pi}} Y_i, \overline{\boldsymbol{\Pi}} \boldsymbol{X}_i): 1 \le i \le n \}$.

Why does the CCE approach not work in the high-dimensional case where $p$ is large? In particular, why not simply estimate $\beta$ by applying lasso rather than least squares techniques to the sample of transformed data $\{ (\overline{\boldsymbol{\Pi}} Y_i, \overline{\boldsymbol{\Pi}} \boldsymbol{X}_i): 1 \le i \le n \}$? The problem is that the CCE proxy $\overline{\boldsymbol{\Pi}}$ breaks down completely in high dimensions. To see this, consider the following situation:

In this situation, the column space of $\overline{\boldsymbol{X}}$ is con\-si\-der\-ably larger than the column space of $\boldsymbol{F}$. In particular, the columns of $\overline{\boldsymbol{X}}$ span the whole space $\mathbb{R}^T$. As a consequence, $\overline{\boldsymbol{\Pi}} = \boldsymbol{I} - \overline{\boldsymbol{X}} (\overline{\boldsymbol{X}}^\top \overline{\boldsymbol{X}})^{-} \overline{\boldsymbol{X}}^\top$ is the projection matrix onto the orthogonal complement of $\mathbb{R}^T$, which is the linear space consisting of the null vector only. Put differently, $\overline{\boldsymbol{\Pi}}$ is the null matrix (that is, the matrix with the entry $0$ everywhere), which is obviously an extremely poor proxy of the projection matrix $\boldsymbol{\Pi}$. These observations point to a general shortcoming of the CCE approach which is well-known in the literature Karabiyik2017: If $p$ is comparably large, the column space of $\overline{\boldsymbol{X}}$ tends to be much larger than the column space of $\boldsymbol{F}$, implying that $\overline{\boldsymbol{\Pi}}$ is a poor proxy of $\boldsymbol{\Pi}$. In the worst case scenario, the columns of $\overline{\boldsymbol{X}}$ span the whole space $\mathbb{R}^T$, which means that $\overline{\boldsymbol{\Pi}} = \boldsymbol{0}$. This worst case occurs whenever $\overline{\boldsymbol{X}} \in \mathbb{R}^{T \times p}$ has full rank $T$. Importantly, this may already happen when $p \ge T$. Hence, the CCE approach runs into trouble not only in the high-dimensional case where $p$ is much larger than $n$ and $T$, but already when $p$ has size comparable to $T$. The larger $p$, the more likely it is that the matrix $\overline{\boldsymbol{X}}$ has rank $T$. Hence, in high dimensions, the proxy $\overline{\boldsymbol{\Pi}}$ of the CCE approach is not reliable and can be expected to break down frequently.

Definition of the HD-CCE estimator

We now construct a proxy of the unknown projection matrix $\boldsymbol{\Pi}$ which does work in high dimensions and build an estimator of $\beta$ based on it.

Step 1: Estimation of the unknown number of factors $\boldsymbol{K}$

Compute the $p \times p$ matrix $\widehat{\boldsymbol{\Sigma}} = T^{-1} \sum_{t=1}^T \overline{X}_t \overline{X}_t^\top$ from the cross-sectional averages $\overline{X}_t$ and perform an eigendecomposition of $\widehat{\boldsymbol{\Sigma}}$, which yields the eigenvalues $\widehat{{\color{black}{\psi}}}_1 \ge \widehat{{\color{black}{\psi}}}_2 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_p \ge 0$ and the corresponding orthonormal eigenvectors $\widehat{U}_1,\ldots,\widehat{U}_p$. Estimate the unknown number of factors $K$ by \[ \widehat{K} = \sum_{j=1}^p 1\big(\widehat{{\color{black}{\psi}}}_j \ge \tau\big), \] where $\tau = \tau_{n,T}$ is a threshold parameter that is of slightly smaller order than $p$. Precise technical conditions on $\tau$ can be found in Section (ref) and rules for selecting $\tau$ in practice are discussed in Section (ref).

Step 2: Approximation of the unknown projection matrix $\boldsymbol{\Pi}$

Let $\widehat{\boldsymbol{U}} = (\widehat{U}_1 \ldots \widehat{U}_{\widehat{K}})$ be the matrix of eigenvectors of $\widehat{\boldsymbol{\Sigma}}$ that correspond to the $\widehat{K}$ largest eigenvalues $\widehat{{\color{black}{\psi}}}_1 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_{\widehat{K}}$ and define $\widehat{\boldsymbol{W}} = \overline{\boldsymbol{X}} \widehat{\boldsymbol{U}}$. Approximate $\boldsymbol{\Pi}$ by \[ \widehat{\boldsymbol{\Pi}} = \boldsymbol{I} - \widehat{\boldsymbol{W}} (\widehat{\boldsymbol{W}}^\top \widehat{\boldsymbol{W}})^{-} \widehat{\boldsymbol{W}}^\top. \]

Step 3: Estimation of $\boldsymbol{\beta}$

Run a lasso regression on the transformed data sample $\{ (\widehat{Y}_i, \widehat{\boldsymbol{X}}_i): 1 \le i \le n \}$, where $\widehat{Y}_i = \widehat{\boldsymbol{\Pi}} Y_i$ and $\widehat{\boldsymbol{X}}_i = \widehat{\boldsymbol{\Pi}} \boldsymbol{X}_i$. Specifically, define the lasso estimator of $\beta$ by \[ \widehat{\beta}_{\color{black}{\lambda}} \in \underset{b \in \mathbb{R}^p}{\text{argmin}} \bigg\{ \frac{1}{nT} \sum_{i=1}^n \big\| \widehat{Y}_i - \widehat{\boldsymbol{X}}_i b \big\|^2 + {\color{black}{\lambda}} \|b\|_1 \bigg\}, \] where ${\color{black}{\lambda}} > 0$ is the penalty constant of the lasso.

Heuristic idea behind the HD-CCE estimator

We now give some heuristic arguments why our estimation approach works in high dimensions. We in particular explain why the matrix $\widehat{\boldsymbol{\Pi}}$ defined in Step 2 of the algorithm provides a good approximation to the unknown projection matrix $\boldsymbol{\Pi}$ even when $p$ is very large. Since the heuristics are essentially the same for large and small $T$, we restrict attention to the large-$T$-case.

Our estimation algorithm is based on the following observation: The cross-sectional averages $\overline{X}_t = n^{-1} \sum_{i=1}^n X_{it}$ satisfy a high-dimensional approximate factor model of the form

The error terms $u_t = (u_{t,1},\ldots,u_{t,p})^\top$ in this model are negligible in the sense that $u_{t,j} = o_p(1)$ for any $t$ and $j$ as $n \to \infty$. This directly follows from the fact that under our regularity conditions, $\overline{\Gamma}_j = \Gamma_j + o_p(1)$ and $\overline{Z}_{t,j} = o_p(1)$ for any $t$ and $j$ as $n \to \infty$, where $\Gamma_j$ and $\overline{\Gamma}_j$ denote the $j$-th row of $\boldsymbol{\Gamma}$ and $\overline{\boldsymbol{\Gamma}}$, respectively. Hence, it holds that $\overline{X}_t \approx \boldsymbol{\Gamma} F_t$, or put differently, $\overline{\boldsymbol{X}} \approx \boldsymbol{F} \boldsymbol{\Gamma}^\top$, which means that the variables $\overline{X}_t$ approximately follow a factor model.

In Step 1 of the estimation algorithm, we exploit this observation as follows: As $\overline{X}_t$ satisfies (ref), the matrix $\overline{\boldsymbol{\Sigma}} = \mathbb{E}[T^{-1} \sum_{t=1}^T \overline{X}_t \overline{X}_t^\top]$ is closely related to a high-dimensional covariance matrix in an approximate factor model. Such covariance matrices tend to have spiked eigenvalues as observed and exploited e.g.\ in Fan2013. We thus expect the eigenvalues of $\overline{\boldsymbol{\Sigma}}$ to be spiked as well. More formally, we can show that under our assumptions, the first $K$ eigenvalues of $\overline{\boldsymbol{\Sigma}}$ are (at least) of order $p$ (in the sense of being bounded from below by $cp$ for some positive constant $c$ and sufficiently large $n$), whereas the others are of (much) smaller order (in the sense of being $o(p)$). The eigenvalues $\widehat{{\color{black}{\psi}}}_1 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_p$ of the estimator $\widehat{\boldsymbol{\Sigma}} = T^{-1} \sum_{t=1}^T \overline{X}_t \overline{X}_t^\top$ can be shown to behave similarly: whereas the $K$ largest eigenvalues are of order $p$, the others are of considerably smaller order. This suggests to estimate $K$ by thresholding the eigenvalues of $\widehat{\boldsymbol{\Sigma}}$. In particular, we may work with the estimator $\widehat{K} = \sum_{j=1}^p 1(\widehat{{\color{black}{\psi}}}_j \ge \tau)$ introduced in Step 1 of the algorithm.

In Step 2 of the algorithm, we exploit the observation that $\overline{X}_t$ satisfies an approximate factor model as follows: Let $\widehat{\boldsymbol{U}} = (\widehat{U}_1 \ldots \widehat{U}_{\widehat{K}})$ be the matrix of eigenvectors of $\widehat{\boldsymbol{\Sigma}}$ that correspond to the $\widehat{K}$ largest eigenvalues $\widehat{{\color{black}{\psi}}}_1 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_{\widehat{K}}$. Since $\overline{\boldsymbol{X}} \approx \boldsymbol{F} \boldsymbol{\Gamma}^\top$, it holds that \[ \widehat{\boldsymbol{\Sigma}} = \frac{\overline{\boldsymbol{X}}^\top \overline{\boldsymbol{X}}}{T} \approx \boldsymbol{\Gamma} \Big(\frac{\boldsymbol{F}^\top \boldsymbol{F}}{T} \Big) \boldsymbol{\Gamma}^\top = \boldsymbol{\Gamma} \boldsymbol{\Gamma}^\top, \] where we have used that $\boldsymbol{F}^\top \boldsymbol{F}/T = \boldsymbol{I}_K$ by (ref). Let $\boldsymbol{\Gamma} = \boldsymbol{U} \boldsymbol{D} \boldsymbol{V}^\top$ be the singular value decomposition of $\boldsymbol{\Gamma}$, where the matrices $\boldsymbol{U} \in \mathbb{R}^{p \times K}$ and $\boldsymbol{V} \in \mathbb{R}^{K \times K}$ have orthonormal columns and $\boldsymbol{D}$ is a diagonal matrix which contains the singular values on its main diagonal. With this decomposition, we further obtain that \[ \widehat{\boldsymbol{\Sigma}} \approx \boldsymbol{\Gamma} \boldsymbol{\Gamma}^\top = \boldsymbol{U} \boldsymbol{D}^2 \boldsymbol{U}^\top. \] This suggests that the matrix $\widehat{\boldsymbol{U}}$ of the first $\widehat{K}$ eigenvectors of $\widehat{\boldsymbol{\Sigma}}$ can be regarded as an estimator of the matrix $\boldsymbol{U}$ whose columns are the first $K$ eigenvectors of $\boldsymbol{\Gamma} \boldsymbol{\Gamma}^\top$. So far, we have seen that $\widehat{\boldsymbol{U}} \approx \boldsymbol{U}$ and $\overline{\boldsymbol{X}} \approx \boldsymbol{F} \boldsymbol{\Gamma}^\top$, which taken together yields that

Since $\boldsymbol{V} \boldsymbol{D}$ is invertible under the full-rank condition on $\boldsymbol{\Gamma}$ in (ref), the $K$ columns of the matrix $\boldsymbol{W} := \boldsymbol{F} \boldsymbol{V} \boldsymbol{D}$ span the same linear space as those of $\boldsymbol{F}$. Consequently,

Moreover, since $\boldsymbol{W} \approx \widehat{\boldsymbol{W}} := \overline{\boldsymbol{X}} \widehat{\boldsymbol{U}}$ by (ref), a good proxy of the projection matrix $\boldsymbol{\Pi}$ should be given by \[ \widehat{\boldsymbol{\Pi}} = \boldsymbol{I} - \widehat{\boldsymbol{W}} (\widehat{\boldsymbol{W}}^\top \widehat{\boldsymbol{W}})^{-} \widehat{\boldsymbol{W}}^\top, \] which is the proxy defined in Step 2 of the algorithm.

From the heuristic discussion so far, it follows that $\widehat{\boldsymbol{\Pi}} \boldsymbol{F} \approx \boldsymbol{\Pi} \boldsymbol{F} = \boldsymbol{0}$. Hence, applying the matrix $\widehat{\boldsymbol{\Pi}}$ to the model equation $Y_i = \boldsymbol{X}_i \beta + \boldsymbol{F} \gamma_i + \varepsilon_i$ leads to the transformed (approximate) model equation $\widehat{\boldsymbol{\Pi}} Y_i \approx \widehat{\boldsymbol{\Pi}} \boldsymbol{X}_i \beta + \widehat{\boldsymbol{\Pi}} \varepsilon_i$ for each $i$. Stacking these equations for all $i$, we obtain the (approximate) high-dimensional linear panel regression model \[ \widehat{Y} \approx \widehat{\boldsymbol{X}} \beta + \widehat{\varepsilon} \quad with \quad \widehat{Y} =

, \ \widehat{\boldsymbol{X}} =

and \ \widehat{\varepsilon} =

, \] which does not have any interactive fixed effects in the errors. To obtain an estimator of $\beta$, we apply standard techniques from high-dimensional linear regression to this transformed model. Specifically, we work with lasso techniques, which leads to the estimator $\widehat{\beta}_{\color{black}{\lambda}}$ defined in Step 3 of the algorithm.

A least squares version of the HD-CCE estimator

So far, our discussion has concentrated on the high-dimensional case where $p$ is large and may even exceed the sample size $nT$. However, the CCE approach does not only break down in this high-dimensional setting. It rather becomes unreliable as soon as $p \ge T$. This is particularly problematic when the time series length $T$ is fairly short as often happens in microeconomic applications. In this case, the number of available regressors $p$ easily exceeds $T$, which means that we are faced with the following situation:

where the symbol $a \ll b$ is here used informally to express that $a$ is considerably smaller than $b$.

In the situation given by (ref), the CCE method is essentially inapplicable. Our estimator, in contrast, works perfectly fine. It is also possible to replace it by a least squares version since there is no need to use the lasso when $p \ll nT$. This is done as follows: We construct $\widehat{K}$ and $\widehat{\boldsymbol{\Pi}}$ exactly as described in the first two steps of the estimation algorithm. However, instead of using the lasso in the third step, we apply least squares to the transformed data $\{ (\widehat{Y}_i, \widehat{\boldsymbol{X}}_i): 1 \le i \le n \}$ with $\widehat{Y}_i = \widehat{\boldsymbol{\Pi}} Y_i$ and $\widehat{\boldsymbol{X}}_i = \widehat{\boldsymbol{\Pi}} \boldsymbol{X}_i$. This yields the least-squares-type estimator \[ \widehat{\beta}_{\text{LS}} \in \underset{b \in \mathbb{R}^p}{\text{argmin}} \bigg\{ \frac{1}{nT} \sum_{i=1}^n \big\| \widehat{Y}_i - \widehat{\boldsymbol{X}}_i b \big\|^2 \bigg\}, \] which is nothing else than the lasso $\widehat{\beta}_{\color{black}{\lambda}}$ with ${\color{black}{\lambda}} = 0$.

It depends of course on the specific sizes of $n$, $T$ and $p$ whether it makes more sense to use the lasso $\widehat{\beta}_{\color{black}{\lambda}}$ (with some ${\color{black}{\lambda}} > 0$) or the least squares version $\widehat{\beta}_{\text{LS}}$. If $p$ is only slightly larger than $T$ in the situation given by (ref), there is only a small number of regressors in the model and one may prefer to use the least squares estimator $\widehat{\beta}_{\text{LS}}$. This in particular has the advantage that we do not have to select the penalty parameter ${\color{black}{\lambda}}$. If $p$ is substantially larger than $T$, that is, if there is a comparably large number of regressors in the model, one may prefer to use the lasso instead for the following reasons: The least squares estimator can be expected to be outperformed by penalized least squares methods such as the lasso. Moreover, since the lasso performs not only estimation but also variable selection, it produces results that are easier to interpret.

The desparsified HD-CCE estimator

As is well known, the lasso -- and thus in particular our HD-CCE estimator -- has a very complicated limiting distribution which is hardly tractable. For this reason, it cannot be used for statistical inference in practice. A common way to circumvent this issue is to desparsify or debias the lasso; see vandeGeer2014, JavanmardMontanari2014, ZhangZhang2014 and Belloni2014. In what follows, we demonstrate how desparsified lasso techniques can be applied to our HD-CCE estimator. To do so, we focus on the following inference problem: we want to compute (asymptotic) confidence bands for the coefficient $\beta_j$ of the $j$-th regressor. Our approach to solve this inference problem is as follows.

Step 1: Estimation of the projection matrix

For inference purposes, we need to construct the proxy of the projection matrix $\boldsymbol{\Pi}$ slightly differently than we did for estimation purposes. In particular, we need to replace the matrix $\overline{\boldsymbol{X}} = n^{-1} \sum_{i=1}^n \boldsymbol{X}_i$ by $\overline{\boldsymbol{X}}_{(-j)}$ which results from eliminating the $j$-th column of $\overline{\boldsymbol{X}}$. This helps us to control certain bias terms in the asymptotic theory. Once this replacement is done, the construction proceeds as before: (i) Compute the $(p-1) \times (p-1)$ matrix $\widetilde{\boldsymbol{\Sigma}} = \overline{\boldsymbol{X}}_{(-j)}^\top \overline{\boldsymbol{X}}_{(-j)}/T $ with eigenvalues $\widetilde{{\color{black}{\psi}}}_1 \geq \ldots \geq \widetilde{{\color{black}{\psi}}}_{p-1} \geq 0$ and corresponding eigenvectors $\widetilde{U}_{(1)}, \ldots, \widetilde{U}_{(p-1)}$. (ii) Estimate $\boldsymbol{\Pi}$ by

where $\widetilde{\boldsymbol{W}} = \overline{\boldsymbol{X}}_{(-j)} \widetilde{\boldsymbol{U}}$ and $\widetilde{\boldsymbol{U}} = (\widetilde{U}_{1} \ldots \widetilde{U}_{\widehat{K}})$ with $\widehat{K}$ as defined before.\footnote{It is possible to replace the estimator $\widehat{K}$ by $\widetilde{K} = \sum_{\ell=1}^{p-1} 1\big(\widetilde{{\color{black}{\psi}}}_\ell \ge \tau\big)$ with an appropriately chosen threshold sequence $\tau = \tau_{n,T}$. However, there is no need to do so from a theoretical point of view.} Notably, the matrix $\widetilde{\boldsymbol{\Sigma}}$ depends on $j$. The same holds for the quantities based on it such as $\widetilde{\boldsymbol{\Pi}}$, $\widetilde{\boldsymbol{W}}$ and $\widetilde{\boldsymbol{U}}$ as well as further expressions defined in the subsequent steps. For simplicity of notation, we however suppress the dependence on $j$ throughout.

Step 2: Estimation of the parameter vectors $\boldsymbol{\beta}$ and $\boldsymbol{{\color{black}{\theta}}}$

We estimate $\beta$ as before by applying lasso techniques to the projected sample of data $\{ (\widetilde{Y}_i, \widetilde{\boldsymbol{X}}_i): 1 \le i \le n \}$ with $\widetilde{Y}_i = \widetilde{\boldsymbol{\Pi}} Y_i$ and $\widetilde{\boldsymbol{X}}_i = \widetilde{\boldsymbol{\Pi}} \boldsymbol{X}_i$. This yields the estimator

Analogously, we estimate the parameter vector ${\color{black}{\theta}}$ in the nodewise equation by

where $\widetilde{X}_{i(j)} = \widetilde{\boldsymbol{\Pi}} X_{i(j)}$, $\widetilde{\boldsymbol{X}}_{i(-j)} = \widetilde{\boldsymbol{\Pi}} \boldsymbol{X}_{i(-j)}$ and ${\color{black}{\kappa}}$ is the penalty constant of the lasso.

Step 3: Desparsifying the HD-CCE estimator

Let $\widetilde{{\color{black}{\Delta}}}_i = \widetilde{X}_{i(j)} - \widetilde{\boldsymbol{X}}_{i(-j)} \widetilde{{\color{black}{\theta}}}_{{\color{black}{\kappa}}}$ be the residual vector from the nodewise lasso regression and write $\widetilde{{\color{black}{\Delta}}} = (\widetilde{{\color{black}{\Delta}}}_1^\top,\ldots,\widetilde{{\color{black}{\Delta}}}_n^\top)^\top$. Following the strategy in vandeGeer2014, we define the desparsified HD-CCE estimator of $\beta_j$ by

with $\widetilde{Y} = (\widetilde{Y}_1^\top, \dots, \widetilde{Y}_n^\top)^\top$, $\widetilde{\boldsymbol{X}} = (\widetilde{\boldsymbol{X}}_1^\top \dots \widetilde{\boldsymbol{X}}_n^\top)^\top$ and $\widetilde{X}_{(j)} = ( \widetilde{X}_{1(j)}^\top, \dots, \widetilde{X}_{n(j)}^\top)^\top$. Notably, this estimator is closely related (but not identical) to the double lasso ChernozhukovHansen2022, which amounts to a least squares regression of lasso residuals from the main equation on lasso residuals from the nodewise equation.

Step 4: Inference with the desparsified HD-CCE estimator

To perform inference with the desparsified HD-CCE estimator, we consider the statistic \[ \mathbb{T}_j = \frac{\widetilde{{\color{black}{\Delta}}}^\top \widetilde{X}_{(j)}}{\mathcal{N}}(\widetilde{b}_{j} - \beta_j) \quad \text{with} \quad \mathcal{N} = \sqrt{\sum_{t,t'=1}^T \bigg\{ \sum_{i=1}^n \widetilde{{\color{black}{\Delta}}}_{it} \widetilde{{\color{black}{\Delta}}}_{it'} \mathbb{E}[ \varepsilon_{it} \varepsilon_{it'} ] \bigg\}}. \] This statistic can be treated as approximately standard normal, which is formally justified in Section (ref). Hence,

with $q_\kappa$ the $\kappa$-quantile of the standard normal distribution is an approximate confidence band of level $(1-\alpha)$ for $\beta_j$, that is, $\mathbb{P}( \beta_j \in \mathbb{C}_{j,\alpha}) \approx 1 - \alpha$. As the normalization term $\mathcal{N}$ in the definition of $\mathbb{T}_j$ is not available in practice, we need to replace it by an estimator $\widetilde{\mathcal{N}}$. There are different ways to do so:

In our R package, the normalizations $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{IID}}$, $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{HET}}$ and $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{HAC}}$ are available as different options. In the supplement, we run a number of simulation exercises to explore the performance of the statistic $\mathbb{T}_j$ with these normalizations. In the theoretical analysis of Section (ref), we focus on the i.i.d.\ case and thus on the normalization $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{IID}}$.

Implementation

The HD-CCE estimator $\widehat{\beta}_{\color{black}{\lambda}}$ depends on two tuning parameters: the threshold parameter $\tau$ for the estimation of $K$ and the penalty parameter ${\color{black}{\lambda}}$ of the lasso. The desparsified HD-CCE estimator $\widetilde{b}_j$ additionally involves the penalty parameter ${\color{black}{\kappa}}$ from the nodewise lasso regression. We now discuss how to select these tuning parameters in practice.

Choice of $\boldsymbol{\tau}$

Our estimator of $K$ is defined as $\widehat{K} = \sum_{k=1}^p 1(\widehat{{\color{black}{\psi}}}_k \ge \tau)$, where $\widehat{{\color{black}{\psi}}}_1 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_p$ are the eigenvalues of $\widehat{\boldsymbol{\Sigma}}$ in descending order. It can be shown formally that the eigenvalues $\widehat{{\color{black}{\psi}}}_k$ are of order $p$ for $k \le K$ but of much smaller order for $k > K$. Hence, to ensure that $\widehat{K}$ is a consistent estimator of $K$, we need to choose $\tau$ such that it separates the “large” eigenvalues of order $p$ (that is, those with $k \le K$) from the “small” ones (that is, those with $k > K$). As a practical rule-of-thumb, we regard an eigenvalue $\widehat{{\color{black}{\psi}}}_k$ as “small” if $\widehat{{\color{black}{\psi}}}_k/\widehat{{\color{black}{\psi}}}_1 < \alpha$ with some small $\alpha$ (such as $\alpha = 0.05$ or $\alpha = 0.01$). Put differently, we regard $\widehat{{\color{black}{\psi}}}_k$ as “small” if it is less than $100\cdot\alpha\%$ of the largest eigenvalue $\widehat{{\color{black}{\psi}}}_1$ in size. This rule-of-thumb results in the choice $\tau = \alpha \widehat{{\color{black}{\psi}}}_1$.\footnote{From a theoretical point of view, we need to let $\alpha = \alpha_{n,T}$ slowly go to $0$ with increasing sample size to make sure that $\tau = \alpha \widehat{{\color{black}{\psi}}}_1$ is of somewhat smaller order than $p$ and thus produces a consistent estimator $\widehat{K}$ of $K$. In practice, however, the sample size is fixed, implying that $\alpha$ is a fixed number as well. We thus do not reflect the dependence of $\alpha$ on $n$ and $T$ in the notation.}

The estimator $\widehat{K}$ is closely related to a simple graphical tool that is frequently used in factor analysis: a scree plot which depicts the eigenvalues $\widehat{{\color{black}{\psi}}}_1 \ge \ldots \ge \widehat{{\color{black}{\psi}}}_p$ in descending order. Typically, a large gap or elbow becomes visible in such a plot which allows to distinguish the large eigenvalues from the small ones. The estimator $\widehat{K}$ formalizes this graphical tool by thresholding the eigenvalues. There are many alternatives to the estimator $\widehat{K}$. Determining the number of factors is a well-understood problem in factor analysis. See for example Kapetanios2010 and Onatski2010 as well as Chapter 6 in Jolliffe2002 for an overview of common approaches.

Choice of $\boldsymbol{{\color{black}{\lambda}}}$ for the HD-CCE estimator

We choose the penalty parameter $\lambda$ by a version of cross-validation, the details of which are explained below. Another possibility is to adapt selection methods that are based on the effective noise of the lasso LedererVogt2021 to the setting at hand. Yet another possibility is to adapt the method of Belloni2016. This would, however, require to estimate the factors $F_t$ and the loadings $\gamma_i$, which goes a bit against the philosophy of our approach to eliminate or “project away” the factors rather than estimate them. Generally speaking, it is highly non-trivial to derive theory for data-driven selection of the lasso's tuning parameter already in a plain-vanilla linear model with i.i.d.\ cross-sectional data; see ChetverikovLiaoChernozhukov2021 for cross-validated lasso and LedererVogt2021 for effective noise based methods. We thus take a pragmatic approach to the problem of selecting ${\color{black}{\lambda}}$ in this paper: As in most other theoretical treatments of the lasso in the literature, we regard the penalty parameter ${\color{black}{\lambda}}$ as a deterministic quantity that converges to $0$ at an appropriate rate when deriving our theory. In the empirical part of the paper, we choose ${\color{black}{\lambda}}$ by the following version of $L$-fold cross-validation: Divide the sample $\{(\widehat{Y}_i,\widehat{X}_i): i=1,\ldots,n\}$ of the projected data into $L$ folds $\mathcal{F}_1,\ldots,\mathcal{F}_{L}$, where $\mathcal{F}_{\ell} = \{(\widehat{Y}_i,\widehat{X}_i): i \in \mathcal{I}_\ell \}$ with $\mathcal{I}_\ell = \{ (\ell-1) \lfloor n/L \rfloor + 1,\ldots, \ell \lfloor n/L \rfloor\}$ for $\ell=1,\ldots,L-1$ and $\mathcal{I}_L = \{ (L-1) \lfloor n/L \rfloor + 1,\ldots, n\}$. Then run standard $L$-fold cross-validation over a grid of $\lambda$-values.\footnote{In our R package, we use the grid chosen by the cross-validation function of the glmnet package.}

Choice of $\boldsymbol{{\color{black}{\lambda}}}$ and $\boldsymbol{{\color{black}{\kappa}}}$ for the desparsified HD-CCE estimator

We follow the selection strategy advocated in dezeure2015nodewisechoice. Specifically, the penalty parameter ${\color{black}{\lambda}}$ is chosen by cross-validation as before and the nodewise penalty constant ${\color{black}{\kappa}}$ is selected by the following procedure dezeure2015nodewisechoice:

Data normalization

As the standard lasso, our HD-CCE estimator is not invariant to the scaling of the regressors. In the literature on the lasso, it is common practice to normalize the regressors prior to estimation to have empirically zero mean and unit variance and then to return the estimated coefficients on the original scale. This is, for example, the baseline procedure when fitting the lasso with the very popular R package glmnet. We essentially follow this convention in our R package hdcce: we compute the HD-CCE estimator (and the nodewise estimator for its desparsified variant) from a normalized version of the projected regressors $\widehat{\boldsymbol{\Pi}} \boldsymbol{X}_i$, but we output the coefficient estimates on the original scale. Specifically, we normalize the projected regressors to have empirical variance $1$ (but we do not centre them as the projection approximately centres them anyway). For simplicity, this normalization step is not reflected in our theory, but it is possible to adjust the theory accordingly.

Theoretical results

Assumptions for the analysis of the HD-CCE estimator

The components of model (ref)--(ref) are assumed to satisfy the following regularity conditions:

In the large-$T$-case, we additionally assume that the model variables form weakly dependent time series processes that satisfy the following mixing conditions:

(ref)--(ref) are very similar to the assumptions in Pesaran2006. However, unlike there, we do not impose any linearity or stationarity assumptions on the involved time series. Notably, the final requirement in (ref) according to which $\max_{1 \le k \le K} \{ T^{-1} \sum_{t=1}^T |F_{t,k}|^\nu \} \le C < \infty$ is rather mild. If $\{F_{t,k}: t =1,\ldots,T\}$ were a time series of weakly dependent random variables with sufficiently many moments, then standard concentration bounds would imply that $T^{-1} \sum_{t=1}^T |F_{t,k}|^\nu \le C < \infty$ with probability approaching $1$. Hence, if we think of our factors as realizations of such time series, the final requirement in (ref) will be fulfilled with high probability. Combined with (ref), the moment conditions on the entries of the loading matrices $\boldsymbol{\Gamma}_i$ in (ref) imply that the factors are pervasive (i.e., each factor must drive a sufficiently large number of regressors), which is a common assumption in the literature Onatski2012. It is in principle possible to drop (ref) in the large-$T$-case and to do without any conditions on the time series dependence of the model variables as in the small-$T$-case. However, then we could not fully account for the time series information in the data. As a consequence, we would obtain a slower convergence rate for our estimator of $\beta$. For simplicity, the mixing coefficients in (ref) are assumed to decay to zero exponentially fast. It is possible though to allow for sufficiently fast polynomial decay instead.

Besides the conditions (ref)--(ref) on the model components, we need some restrictions on the dimension $p$ and the sparsity index $s$. In the large-$T$-case, we impose the following conditions on the dimension para\-meters $n$, $T$, $p$, $s$ and $K$:

(ref) essentially says that $p$ is not allowed to grow too quickly in comparison to $n$ and $T$. To better understand the restrictions on $p$, let us consider the special case $n=T$. In this case, the two restrictions of (ref) simplify to $(nT)^{({\color{black}{\nu}}/4)-1} \gg p$. Hence, how fast $p$ can grow in comparison to the sample size $nT$ depends on how many moments ${\color{black}{\nu}}$ the model variables have. If all moments exist, ${\color{black}{\nu}}$ can be chosen as large as desired and $p$ can grow as any polynomial of $nT$. If ${\color{black}{\nu}}$ is quite small in contrast, say ${\color{black}{\nu}} = 8 + \delta$ for some small $\delta > 0$, then $p$ can only grow slightly faster than the sample size $nT$. (ref) imposes constraints on the growth of the sparsity index $s$, that is, on the number of non-zero components of $\beta$. As one can see, $s$ is restricted to grow slightly more slowly than $\min\{n,T\}$. In the special case $n=T$, in particular, $s$ can only grow slightly more slowly than $\sqrt{nT}$. In the small-$T$-case, our conditions on the dimension parameters $n$, $T$, $p$, $s$ and $K$ are as follows:

(ref) puts restrictions on the growth of $p$. Analogously to the large-$T$-case, the more moments ${\color{black}{\nu}}$ exist, the faster $p$ is allowed to grow in comparison to $n$. In particular, if all moments exist, then $p$ can grow as any polynomial of $n$. (ref) imposes constraints on the growth of the sparsity index $s$. As can be seen, the more moments ${\color{black}{\nu}}$ exist, the faster $s$ is allowed to increase. In particular, if all moments exist, then $s$ can grow almost as fast as $\sqrt{n}$. In contrast, if only a few moments exist, say ${\color{black}{\nu}} = 8 + \delta$ for some small $\delta > 0$, then $s$ must grow considerably more slowly than $\sqrt{n}$.

All in all, the above conditions on the dimensions $n$, $T$ and $p$ allow us to deal with a wide range of scenarios (as long as the tails of the model variables are not too thick, i.e., as long as sufficiently many moments ${\color{black}{\nu}}$ exist). In particular, we can deal with “standard low-dimensional” scenarios where $p$ is small and fixed, with “moderately high-dimensional” scenarios where $p$ is fairly large but still smaller than the sample size $nT$ and with “truly high-dimensional” scenarios where $p$ exceeds the sample size $nT$.

Convergence rate of the HD-CCE estimator

We now derive the convergence rate of the HD-CCE estimator $\widehat{\beta}_{\color{black}{\lambda}}$. To formulate the result, we let $\{h_n\}$ be any sequence of positive real numbers which slowly diverges to infinity. For instance, we may choose $h_n = C \log \log n$ with some constant $C > 0$.

The proof of Theorem (ref) is provided in the technical appendices. We briefly give some remarks on the derived convergence rates.

Assumptions for the analysis of the desparsified HD-CCE estimator

In addition to (ref)--(ref), we impose the following conditions on the components in model equations (ref), (ref) and (ref):

Assumptions (ref) and (ref) require the error terms $\varepsilon_{it}$ and $u_{it}$ to be independent both across $i$ and $t$. We impose these rather strong independence conditions to avoid certain complications in the extremely technical derivation of the limit distribution of the desparsified HD-CCE estimator. However, we conjecture that it is possible to weaken these conditions (in particular, to allow for weak dependence across $t$ and non-identical distributions across $i$). In the simulations, we assess the performance of our inference methods in situations where the errors are not i.i.d. The assumptions on the factor loadings $\nu_i$ in (ref) are rather mild.

Conditions (ref)--(ref) and (ref)--(ref) on the dimension para\-meters $n$, $T$, $p$, $s$ and $K$ are sufficient to guarantee a good convergence behaviour of the HD-CCE estimator. For inference purposes, however, we require additional constraints on them. Specifically, in the large-$T$-case, we assume the following:

The first part of (ref) requires that $T$ diverges somewhat more slowly than $n$. Hence, for our inference theory, we need to restrict the growth of $T$ more strongly than for the convergence analysis of the HD-CCE estimator. The second part of (ref) is fulfilled, e.g., if $p$ grows at most polynomially in $nT$ and the number of moments ${\color{black}{\nu}}$ of the model variables is sufficiently large. Given (ref), the growth conditions on the sparsity indices $s$ and $\|{\color{black}{\theta}}\|_0$ in (ref) are satisfied if $s$ and $\| {\color{black}{\theta}} \|_0$ diverge a bit more slowly than $\sqrt{nT}$. In the small-$T$-case, we impose additional constraints similar to those in (ref):

Condition (ref) is satisfied if $p$ grows at most polynomially in $n$ and the number of moments ${\color{black}{\nu}}$ is large enough. Moreover, if ${\color{black}{\nu}}$ can be chosen as large as desired (meaning that all moments of the model variables exist), the sparsity indices $s$ and $\|{\color{black}{\theta}}\|_0$ can grow almost as fast as $\sqrt{n}$ under condition (ref).

We finally need to strengthen conditions (ref) and (ref) a bit. We in particular impose the following constraints additional to them:

As in the analysis of the HD-CCE estimator, we could impose the restricted eigenvalue conditions on the matrices $\boldsymbol{X}^\perp$ and $\boldsymbol{X}_{(-j)}^\perp$ (or on the idiosyncratic parts $\boldsymbol{Z}$ and $\boldsymbol{Z}_{(-j)}$) rather than directly on $\widetilde{\boldsymbol{X}}$ and $\widetilde{\boldsymbol{X}}_{(-j)}$ at the cost of additional technical arguments. However, as these technical arguments are completely analogous to those in the analysis of the HD-CCE estimator, we work with assumptions (ref) and (ref) for the sake of simplicity.

Limit distribution of the desparsified HD-CCE estimator

We now show that the desparsified HD-CCE estimator is asymptotically normal both in the large-$T$ and the small-$T$-case. As before, we let $\{h_n\}$ be any sequence of positive real numbers which slowly diverges to infinity (e.g., $h_n = C \log \log n$ with some constant $C > 0$).

The proof of Theorem (ref) is given in the technical appendices.

Simulations

In the supplementary material, we evaluate our estimation and inference methods by Monte Carlo experiments. In the first part of the simulation study, we examine the estimation performance of the HD-CCE estimator (and its least squares version) in various low- and high-dimensional scenarios. The Monte Carlo experiments show that the estimator performs well, being almost as accurate as certain oracle methods which presuppose knowledge of the true projection matrix $\boldsymbol{\Pi}$ and/or the true active set $S = \{j: \beta_j \ne 0 \}$. In the second part, we evaluate the size and power properties of the inference procedures that are based on the desparsified HD-CCE estimator. The simulation results demonstrate accurate size control as well as good power numbers in both low and high dimensions. Finally, in the third part, we run various robustness checks: we demonstrate that our methods are robust to moderate overestimation of $K$, we investigate what happens when some factors are much stronger than others, and we demonstrate that our inference methods are robust to heteroskedasticity and serial correlation of the idiosyncratic errors.

Empirical study

We apply our methods to a financial dataset that has been the subject of much previous work. The goal is to identify firm characteristics that are informative about the firm's stock return. We suppose that the monthly excess stock return $R_{it}$ of firm $i$ at time $t$ satisfies the model equation

for $1 \le i \le n$ and $1 \le t \le T$, where the latent factors $F_{t,k}$ capture the comovement of returns, $\mu_i$ is a stock-specific mean and {$C_{i,{t-1}}$} is a vector of “slowly moving” firm characteristics observed at time point $t-1$. Using the notation $\gamma_{i,0} := \mu_i$ and $F_{t,0} := 1$ for all $t$, model (ref) can be reformulated as

and thus corresponds to our main model equation (ref) with $Y_{it}=R_{it}$ and $X_{it}=C_{i,{t-1}}$. Daniel1997 called (ref) the characteristic-based pricing model and interpret $\mu_{it}:=\mu_i + \beta^{\top}C_{i,{t-1}}$ as the expected return, although they specified only a scalar attribute $C_{i,{t-1}}$ at a time and work with observed Fama French factors. Chen2023 work with a similar model and consider a large number of firm-specific characteristics, but their methodology is portfolio sorting based. Our work is also closely related to Green2017 whose data we use. Model (ref) is in the spirit of much recent “machine learning” work in finance, see e.g.\ Gu et al.\ (2020, eq.\ 22 and 23) \nocite{Gu2020} and Nagel (2021, eq.\ 3.1 and 3.2) \nocite{Nagel2021} and the discussion therein.

We apply model (ref) to a sample of large cap stocks {($n=29$)} from April 2017 to March 2022 {($T=60$)} that are or recently were constituent stocks of the Dow Jones Industrial Average. The names and ticker symbols of the stocks are given in Table (ref). The firm characteristics we use in the application including a brief description are listed in Table (ref). They comprise a subset of the characteristics collected by Green2017 for the period 1980--2014, which have been extended to 2022 by Shaoran Li, who kindly shared the data with us. We refer to the appendix of Green2017 for details on their precise definitions. In comparison to the full set of $102$ characteristics collected by Green2017, we drop two characteristics that are constant over the sample period, nine characteristics that are calculated from past stock returns, and one characteristic that is zero for most firms, leaving us with a total of $p=90$ characteristics.\footnote{The two constant characteristics were ipo, indicating an initial stock issue, and sin, indicating stocks of firms that operate in a sinful industry. Our sample contained neither initial stock issues nor firms from a sinful industry. The nine characteristics that are calculated from past returns include so-called momentum measures that correspond to the cumulative past returns over a given period (mom1m, mom6m, mom12m, mom36m), the change in certain momentum measures (\textit{chmom}), past returns at specific time points (\textit{ear}), the maximum daily return (\textit{maxret}) and the variation of daily returns (\textit{retvol}) over the past month, as well as an illiquidity measure (\textit{ill}). Finally, dividend initiation (\textit{divi}) is zero for nearly all firms not allowing to run the nodewise regressions to calculate the desparsified estimate for this characteristic.} The data are sourced from CRSP, Compustat and I/B/E/S, missing values are imputed using Freyberger2024, the industry adjusted size measure \textit{mve_ia} is scaled by $10^{-4}$ and the weighted number of days without a trade (\textit{zerotrade}) is scaled by $10^{5}$ to bring their values more in line with the other characteristics. Green2017 do Fama Macbeth regressions (averaging across time the OLS coefficients of cross-sectional regressions on all characteristics) using data on all common stocks on the NYSE, AMEX, and NASDAQ with available Compustat accounting data over the period 1980--2014 and the subperiods 1980--2003 and 2003--2014, that is, they had both large $n$ and large $T$ and $p$ smaller than either dimension. They found that over the full period only 12 of the considered characteristics consistently predicted returns for non-microcap stocks. However, return predictability sharply declined in 2003, reducing independent predictors to just two characteristics after that year. We consider only large cap stocks and only a five year period to mitigate nonstationarity issues, which is consistent with a lot of standard practice in this literature.

In what follows, we estimate the parameter vector $\beta$ in model (ref) by our methods and compute pointwise confidence bands for the coefficients $\beta_j$. Notably, the firm characteristics $C_{i,t-1}$ are lagged by one time period, which requires some slight modifications of our methodology. We first present and interpret the estimation results and then provide the details on how to modify and implement our methods.

Estimation results

Our estimates of the coefficients $\beta_j$ and their confidence bands are presented in Figures (ref) and (ref). For better interpretability, we report the scaled coefficients $\widehat{\beta}_{\lambda,j}^{\text{sc}} := \widehat{\beta}_{\lambda,j} \cdot \widehat{s}_j$, where $\widehat{s}_j^2$ is the empirical variance of the $j$-th firm characteristic. We thus normalize the firm characteristics to have empirical variance $1$ and report the coefficients on the resulting scale. Figure (ref) shows the non-zero coefficients $\widehat{\beta}_{\lambda,j}^{\text{sc}}$ ordered by their absolute size, while Figure (ref) shows all estimated coefficients $\widehat{\beta}_{\lambda,j}^{\text{sc}}$ including the zero ones along with their 95% confidence bands.\footnote{In Figure (ref), it is possible that a coefficient $\widehat{\beta}_{\lambda,j}^{\text{sc}}$ is not included in its confidence band. The reason is as follows: the confidence band is computed from the desparsified, i.e., the debiased version of the HD-CCE estimator and the bias of the HD-CCE estimator may be so strong that it lies outside the band. In particular, a coefficient $\widehat{\beta}_{\lambda,j}^{\text{sc}}$ may be equal to zero, while the corresponding confidence interval does not include zero. In Figure (ref), this happens in four cases: chnanalyst, ms, orgcap and rd_mve.}

We share some findings with Green2017. Namely, we find only a small number of characteristics $j$ with a significant non-zero coefficient $\widehat{\beta}_{\lambda,j}^{\text{sc}}$. However, the characteristics we find hardly overlap with those identified by Green2017, albeit they used quite different methodology and their data covered an earlier time period. Nevertheless, the signs of our signi\-fi\-cant non-zero coefficients are in line with prior literature in most cases. The size effect ($mve$) is strongly negative Banz1981 with a positive offset from the industry-adjusted size variable $mve\_ia$; both effects are statistically significant at the $1\%$ level and are amongst the variables with the largest coefficients. Leverage ($lev$) has a negative and statistically significant effect on excess returns at the $5\%$ level, while illiquidity ($baspread$) has a positive effect and is significant at the $1\%$ level, consistent with much earlier work on liquidity premiums Amihud1986. Finally, $beta$ has a positive effect on excess returns, being significant at the $5\%$ level.

Implementation details

We observe the data sample $\{ (Y_{it},C_{it}): 1 \le t \le T, \, 1 \le i \le n \}$ and consider the model consisting of the two equations $R_{it} = \beta^\top C_{i,t-1} + \sum_{k=0}^K \gamma_{i,k} F_{t,k} + \varepsilon_{it}$ and $C_{i,t-1} = \boldsymbol{\Gamma}_i F_{t-1} + Z_{i,t-1}$ for $2 \le t \le T$ and $1 \le i \le n$, where $F_t = (F_{t,0}, F_{t,1},\ldots,F_{t,K})^\top$ and $F_{t,0} := 1$ for all $t$ is a time-constant factor. These equations can equivalently be expressed as

for $1 \le i \le n$, where $\boldsymbol{F} \in \mathbb{R}^{T \times (K+1)}$ is the factor matrix with $t$-th row $F_t$, $\boldsymbol{C}_i = (C_{i1} \ldots C_{iT})^\top \in \mathbb{R}^{T \times p}$ is the matrix of firm characteristics and for a general matrix $\boldsymbol{A} \in \mathbb{R}^{T \times q}$, we let $\boldsymbol{A}_{k:\ell}$ be the matrix which results from $\boldsymbol{A}$ by deleting the rows $1,\ldots,k-1$ and $\ell+1,\ldots,T$. Notably, the error structure in (ref) depends on the factor matrix $\boldsymbol{F}_{2:T}$, whereas the firm characteristics in (ref) depend on the factor matrix $\boldsymbol{F}_{1:(T-1)}$ lagged by one time period. To take this into account, we slightly modify the projection matrix $\widehat{\boldsymbol{\Pi}}$. Specifically, we define $\widehat{\boldsymbol{\Pi}}$ to be the projection onto the orthogonal complement of the linear space $\mathcal{L} = \{ (\boldsymbol{1}, \widehat{\boldsymbol{W}}_{1:(T-1)}, \widehat{\boldsymbol{W}}_{2:T}) \, v: v \in \mathbb{R}^{2\widehat{K}+1} \}$ which is spanned by the constant factor $\boldsymbol{1} = (1,\ldots,1)^\top \in \mathbb{R}^{T-1}$ and the columns of $\widehat{\boldsymbol{W}}_{1:(T-1)}$ and $\widehat{\boldsymbol{W}}_{2:T}$. Here, $\widehat{\boldsymbol{W}}$ is defined exactly as in Section (ref) with $\overline{\boldsymbol{X}} = \overline{\boldsymbol{C}}$ and $\overline{\boldsymbol{C}} = n^{-1} \sum_{i=1}^n \boldsymbol{C}_i$ the cross-sectional average of the matrices $\boldsymbol{C}_i$. With this choice of the projection matrix $\widehat{\boldsymbol{\Pi}}$, we compute the HD-CCE estimator exactly as described in Section (ref). The desparsified HD-CCE estimator is modified analogously. Our theory can be easily adapted to these modifications.

We implement our estimators as follows: The tuning parameter $\tau$ is chosen as suggested in Section (ref), in particular, $\tau=\alpha\widehat{\psi}_1$ with $\alpha=0.01$. This results in the estimate $\widehat{K}=2$. The penalty constant ${\color{black}{\lambda}}$ of the HD-CCE estimator $\widehat{\beta}_\lambda$ is chosen by a leave-one-firm-out version of the cross-validation procedure outlined in Section (ref) (i.e., $L = n$). Moreover, for each $j$, the penalty constant $\kappa$ of the nodewise lasso $\widetilde{\theta}_\kappa$ is chosen as detailed in Section (ref) (with leave-one-firm-out cross-validation). The confidence bands depicted in Figure (ref) are computed according to (ref), where we estimate $\mathcal{N}$ by $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{HET}}$. We use the heteroskedasticity robust normalization $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{HET}}$ rather than $\widetilde{\mathcal{N}}^{\hspace{1pt} \text{HAC}}$ because the idiosyncratic errors $\varepsilon_{it}$ in equation (ref) can be interpreted as idiosyncratic return shocks that are uncorrelated over time. The stars $^*$, $^{**}$ and $^{***}$ after the variable names in Figures (ref) and (ref) indicate whether a variable is significant at the $10\%$ , $5\%$ and $1\%$ level, respectively, i.e., whether the corresponding confidence band excludes zero. In order to make the estimated coefficients (and the corresponding confidence bands) better interpretable, we scale them by the empirical variances of the firm characteristics, i.e., we report the scaled versions $\widehat{\beta}_{\lambda,j}^{\text{sc}}$ and the correspondingly scaled bands.

Concluding remarks

In this paper, we have developed new estimation and inference methods for high-dimensional panel data models with interactive fixed effects. Our methods rely on the following general idea: rather than estimating the unobserved factor structure, we eliminate the factors from the model by a projection. Our methods can thus be regarded as a high-dimensional analogue of the CCE approach which is frequently used in the standard low-dimensional case. The projection device of the CCE approach breaks down completely in high dimensions and a simple fix is not possible. One of the main contributions of the paper is to come up with a projection device which works in both low and high dimensions. This device can be combined with high-dimensional regression techniques to obtain an estimator of the unknown parameter vector. We have focused on (desparsified) lasso techniques, but it is in principle possible to work with other techniques (Dantzig selector, SCAD, etc.). In our theoretical analysis, we have derived the convergence rate of our HD-CCE estimator and the asymptotic distribution of its desparsified version. Of course, this is only a first step towards a comprehensive theory. In what follows, we briefly outline some avenues for future research.

Variable selection

Another interesting issue besides parameter estimation and inference is variable selection. As the original lasso itself, our HD-CCE estimator will be selection consistent only under extremely strong conditions. To get better variable selection properties, it should be possible to combine our HD-CCE approach with adaptive lasso techniques.

Overestimating the number of factors $\boldsymbol{K}$

In our theory, we have assumed that we can consistently estimate the unknown number of factors $K$ (and we have provided a consistent estimator $\widehat{K}$). We conjecture that the theory can be extended to estimators $\widehat{K}$ with the property that $K \le \widehat{K} \le K_{\max}$ with probability tending to $1$, where $K_{\max}$ is a given upper bound on the number of factors. Such an extension would be formal proof that (moderate) overestimation of the number of factors is indeed unproblematic, as suggested by the simulation evidence in the supplementary material.

Dealing with nonlinear transformations

Suppose we observe a sample of panel data $\{(Y_{it}, X_{it}^{\text{raw}}): 1 \le t \le T, \, 1 \le i \le n \}$, where $X_{it}^{\text{raw}} = (X_{it,1}^{\text{raw}}, \ldots, X_{it,p_0}^{\text{raw}})^\top$ is a vector of $p_0$ directly observed variables. Rather than only using the raw variables as regressors in the model, we would also like to include interactions $X_{it,j}^{\text{raw}} \cdot X_{it,k}^{\text{raw}}$ and nonlinear transformations such as polynomials $(X_{it,j}^{\text{raw}})^q$. Collecting all of the resulting regressors -- the raw, the interacted and the transformed variables -- in a long vector $X_{it} = (X_{it,1},\ldots,X_{it,p})^\top$, we consider the high-dimensional model $Y_{it} = \beta^\top X_{it} + \gamma_i^\top F_t + \varepsilon_{it}$. As before, we assume that the observed variables $X_{it}^{\text{raw}}$ satisfy an approximate factor model of the form (ref), i.e., $X_{it}^{\text{raw}} = \boldsymbol{\Gamma}_i F_t + Z_{it}$. However, this factor structure is in general not preserved when the variables $X_{it}^{\text{raw}}$ are transformed nonlinearly. Hence, we cannot assume that all covariates $X_{it,j}$ satisfy an approximate factor model and thus have to modify our HD-CCE estimation algorithm from Section (ref).

We suggest to proceed as follows: we run Steps 1 and 2 of the algorithm on the raw variables $X_{it}^{\text{raw}}$ only, i.e., we construct $\widehat{K}$ and the projection matrix $\widehat{\boldsymbol{\Pi}}$ on the basis of $X_{it}^{\text{raw}}$. We then apply Step 3 with the thus constructed projection matrix. In this way, we can easily accommodate interactions and nonlinear transformations. Our theoretical results on the HD-CCE estimator from Theorem (ref) should still be valid after this modification, provided that the projected design matrix $\boldsymbol{X}^\perp$ satisfies a restricted eigenvalue (RE) condition as detailed in Section (ref). However, showing that $\boldsymbol{X}^\perp$ indeed satisfies such an RE condition is extremely difficult, in particular, much more difficult than in the case analyzed in this paper. Moreover, our desparsification strategy needs to be adapted properly. A natural way would be to partition the covariates $X_{it,j}$ for $j=1,\ldots,p$ into groups. In particular, neglecting interactions, we may put all transformations of a given variable $X_{it,\ell}^{\textnormal{raw}}$ into one group and make inference on the resulting groups of coefficients. It is, however, not straightforward to extend our desparsified HD-CCE procedure to do so. The main issue is this: even if $X_{it,\ell}^{\textnormal{raw}}$ satisfies a nodewise regression equation with an interactive fixed effects error structure similar to (ref), nonlinear transformations of $X_{it,\ell}^{\textnormal{raw}}$ will in general not do so.

Acknowledgements

We thank the Editor and three referees for helpful comments that greatly improved the paper. We thank Alexei Onatski and Hashem Pesaran for insightful discussions and Shaoran Li for supplying the data used in Section (ref). Financial support by the DFG (German Research Foundation) -- project number 501082519 -- is gratefully acknow\-ledged. Computing time granted on the supercomputer CLAIX at RWTH Aachen as part of the NHR4CES infrastructure is thankfully acknowledged as well.

{ {0.1em} }

\allowdisplaybreaks[3]

Throughout the appendices, we let $c$ and $C$ denote generic positive constants that may take a different value on each occurrence. The symbols $c_j$ and $C_j$ with subscript $j$ (which may be either a natural number or a letter) are specific constants that are defined in the course of the appendices. Unless stated differently, the constants $c$, $C$, $c_j$ and $C_j$ depend neither on the dimensions $n$, $T$, $p$ nor on the sparsity index $s$. To emphasize that they do not depend on any of these parameters, we sometimes refer to them as absolute constants.

\setcounter{equation}{0}