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.
64,645 characters · 8 sections · 41 citation commands
\newtheorem{corollary}{Corollary} \newtheorem{definition}{Definition} \newtheorem{lemma}{Lemma} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newtheorem{theorem}{Theorem} \newtheorem{assumption}{Assumption} \newtheorem{example}{Example}
\numberwithin{corollary}{section} \numberwithin{definition}{section} \numberwithin{equation}{section} \numberwithin{proposition}{section} \numberwithin{remark}{section} \numberwithin{theorem}{section}
\allowdisplaybreaks[4]
With the emergence of the data-rich world, numerous applications in many fields, including business and economics, focus on panel data regressions within a high-dimensional framework, where the number of variables can be very large and potentially exceed the sample size. For example, in the realm of asset pricing, academic researchers and financial analysts have put forth over five hundred risk factors and individual firm characteristics (which continue to grow) to explain the cross-sectional relationship of stock return (feng2020taming). Given the necessity for innovative statistical methodologies to dissect such data, there is an increasing interest in discussing regularized estimation of high-dimensional (HD) panel data models (e.g., vogt2022cce,babii2022machine,belloni2023high). Nonetheless, achieving robust estimation and inference using HD panel data remains a challenge, particularly when idiosyncratic errors present non-Gaussianity, time series autocorrelation (TSA), and cross-sectional dependence (CSD). Our main goal is to enrich the relevant literature (e.g., FLY2015, KP2019,gupta2023robust) with a complete toolkit for conducting robust inference about the parameters of interest involved in HD panel data models associated with complex dependence structure.
That said, in this paper, (i) we consider HD panel data models in which the idiosyncratic errors exhibit TSA, CSD, heteroscedasticity, as well as heavy--tailed behaviors; (ii) we debias the least absolute shrinkage and selection operator (LASSO) estimator to achieve asymptotic normality, and explore the estimation of HD long--run covariance matrix using a thresholded estimator; and (iii) we accommodate such scenarios that the number of regressors may increase more rapidly than the sample size.
To achieve these goals with the aforementioned data features, we establish two Nagaev--types of concentration inequalities, one for partial sum and the other for quadratic form, subject to a set of easily verifiable conditions. Leveraging these two inequalities, we derive a non-asymptotic bound for the LASSO estimator, achieve asymptotic normality via the nodewise LASSO regression, and obtain a sharp convergence rate for the thresholded heteroskedasticity and autocorrelation consistent (HAC) estimator. Furthermore, we considerably generalize the main ideas to investigate a class of widely used HD panel data models (e.g., vogt2022cce, belloni2023high) associated with interactive effects.
Up to this point, it is worth commenting on some key references, and outlining our contributions to the relevant literature accordingly.
The rest of this paper is organized as follows. Section (ref) presents the main results of this paper, and Section (ref) showcases their practical relevance by considering a HD panel data model with interactive effects. Section (ref) conducts Monte Carlo simulations to examine the theoretical findings. Section (ref) consists of the empirical results by studying a HD asset pricing model. Section (ref) concludes. Appendix (ref) verifies several assumptions in the paper. A numerical implementation procedure is presented in Appendix (ref). Many other technical details are given in Appendix (ref)--(ref) and online Appendices B and C.
Before proceeding further, we introduce a set of notations which are repeatedly used throughout. For any $\mathbf{x} \in \mathbb{R}^{n}$, $|\mathbf{x}|_p\coloneqq (\sum_{i=1}^{n}|x_i|^p)^{1/p}$, in which $x_i$ stands for the $i^{th}$ element of $\mathbf{x}$; for a random vector $\mathbf{v}$, $\|\mathbf{v}\|_p\coloneqq (E|\mathbf{v}|_{p}^p )^{1/p}$; $|\cdot|$ denotes the absolute value of a scalar or the cardinality of a set; for an $m\times n$ matrix $\mathbf{A} = (A_{ij})_{ i\leq m, j\leq n}$, $|\mathbf{A}|_{p}$ denotes the induced $\ell_p$ matrix norm such that $|\mathbf{A}|_{p} \coloneqq \max_{\mathbf{x} \neq \mathbf{0}}|\mathbf{A}\mathbf{x}|_p/|\mathbf{x}|_p$; $|\mathbf{A}|_{\mathrm{max}} \coloneqq\max_{ i\leq m,j\leq n}|A_{ij}|$ denotes the matrix element-wise max norm; $|\mathbf{A}|_F$ denotes the Frobenius norm; $|\mathbf{A}|_{*} \coloneqq \sum_{k=1}^{\min(m,n)}\psi_{k}(\mathbf{A})$ denotes the nuclear norm, where $\psi_{k}(\cdot)$ is the $k^{th}$ largest singular value of a matrix; $\psi_{\mathrm{max}}(\mathbf{A})$ and $\psi_{\mathrm{min}}(\mathbf{A})$ denote the largest and smallest singular values of $\mathbf{A}$, respectively; for two conformable matrix $\mathbf{A}$ and $\mathbf{B}$, let $\mathbf{A}\circ \mathbf{B}$ be the Hadamard product; $\mathbf{I}_a$ stands for an $a\times a$ identity matrix; let $\mathbf{1}_a$ stand for a $a\times 1$ vector of ones; let $\to_P$, $\to_D$ and $=_D$ denote convergence in probability, convergence in distribution and equal in distribution, respectively; for two constants $a$ and $b$, $a \asymp b$ stands for $a=O(b)$ and $b=O(a)$; for a given positive integer $L$, let $[L]=\{1,2,\ldots, L\}$.
We start this section by introducing a set of HD panel data. In what follows, $i\in [N]$ and $t\in [T]$ always index individuals and time periods respectively, and both $(N,T)$ may diverge. Additionally, let $j\in [d]$ index the variables, and $d$ may diverge along with $(N,T)$. Accordingly, for a given time period $t\in [T]$ and a given variable $j\in [d]$, we denote the following $N\times 1$ vector:
where $\{\pmb{\varepsilon}_t\}_{t \in \mathbb{Z}}$ are a sequence of independent and identically distributed (i.i.d.) random vectors, and each $u_{ji}(\cdot)$ is a measurable function. Define the coupled version of $\mathcal{F}_t$ in order to measure dependence: $\mathcal{F}_t^*= (\pmb{\varepsilon}_t, \ldots ,\pmb{\varepsilon}_{1},\pmb{\varepsilon}_{0}^*,\pmb{\varepsilon}_{-1},\ldots)$, where $\pmb{\varepsilon}_{0}^*$ is an i.i.d. copy of $\{\pmb{\varepsilon}_t\}$. Let us introduce Assumption 1 below.
Assumption (ref) generalizes the functional dependence measure of wu2005nonlinear to a panel data setting, allowing for a wide range of data generating processes (DGPs). For example, one often chooses $\mathbf{w} = (1,0,\ldots, 0)^\top$ and $\mathbf{w} = (\frac{1}{\sqrt{N}},\ldots, \frac{1}{\sqrt{N}})^\top$ in some settings (e.g., Pesaran2006; bai2009panel). This condition regulates the decay rate of temporal dependence after taking cross--sectionally weighted averages. By doing so, $u_{ji}(\mathcal{F}_t)$'s are allowed to exhibit cross--sectional dependence of various unknown forms. In Example (ref) of Appendix (ref), we show the flexibility of Assumption (ref).
Assumption (ref) is sufficient for us to establish two Nagaev--types of sharp concentration inequalities, allowing for correlation over both dimensions and non-Gaussianity. These results will be used to study LASSO and HD long-run covariance estimation respectively, and they are also generally applicable to deal with some other scenarios.
Notably, $\mu_u$ is guaranteed to be fixed under Assumption (ref), and $M$ can be a sufficiently large constant. With this lemma in hand, we are ready to present the benchmark model and investigate its estimation and inference.
In this section, we consider the following benchmark model:
where $\mathbf{X}_t$ is a $N\times d$ matrix such that
Accordingly, $\mathbf{y}_{t}=[y_{1t},\ldots, y_{Nt}]^\top$, $\mathbf{e}_{t}=[e_{1t},\ldots, e_{Nt}]^\top$, and $\pmb{\beta}_0=[\beta_{0,1},\ldots,\beta_{0,d}]^\top$. Moreover, we let $J\coloneqq\{ j\in [d]\mid \beta_{0,j} \neq 0\}$ and $J^c \coloneqq [d] \setminus J$, and accordingly denote $s\coloneqq |J| $ and $\beta_{\min}\coloneqq \min_{j\in J} |\beta_{0,j}|$.
Model (ref) has a simple form for us to understand how dependence over both dimensions influences typical HD analysis, such as babii2022machine and baek2023local. It is straightforward to allow addictive fixed effects by using the usual demean procedure. In Section (ref), we extend the main ideas to discuss a class of HD panel data models associated with interactive fixed effects, which as mentioned in the introduction has drawn attentions (e.g., vogt2022cce,belloni2023high, and references therein) in recent years.
Before proceeding further, we introduce some suitable regularity conditions on the regressors and errors permitting non-Gaussianity, TSA, CSD, and heteroscedasticity.
Note that Assumption (ref).1 is similar to Assumption (ref). Assumption (ref).2 formulates the compatibility condition, and enriches the literature by regulating the singular values of the second moment of $\mathbf{x}_{it}$. Typically, the literature imposes such a condition on a sample covariance matrix (e.g., Assumption A2 in chernozhukov2021lasso).
Note also that Assumption (ref) allows that the elements of $\mathbf{x}_{it}$ can be both serially dependent on the time--series dimension and cross-sectionally dependent on the cross--sectional dimension. While we do not consider the case where there are dynamic features involved in our panel data, it is possible to modify our method and theory to allow for strong cross--sectional dependence and lagged dependent variables in $\mathbf{x}_{it}$ under strict exogeneity conditions on the errors. We plan to leave such topics for future research.
We now come back to investigate the model (ref). The first two steps are about penalized estimation and variable selection.
\hrule
\hrule
Given $g_j = \frac{1}{|\widehat{\beta}_j|}$ or $g_j = \frac{\omega^*}{\max(|\widehat{\beta}_j|,\ \omega^*)}$ with a threshold $\omega^*$, model (ref) becomes either an adaptive LASSO (zou2006adaptive) or a conservative LASSO (caner2018asymptotically). In connection with Lemma (ref), we are now able to establish the following results.
Under the condition:
Lemma (ref) presents the non-asymptotic error bounds of Step 1, and implies consistent variable selection of Step 2. Notably, Lemma (ref) does not require specific tail behavior, and model (ref) infers that $d$ may grow at a polynomial order of the sample size.
If we impose further, for example, an exponential tail assumption, $d$ can even diverge at an exponential rate as in van2014asymptotically, wherein an i.i.d. assumption and an exponential tail are adopted to facilitate the development. There is a trade-off between the tail behavior of the error components and the divergence rate of $d$.
To proceed, we start to approximate the following two matrices:
The second one is a HD long-run covariance matrix permitting heavy-tailed errors. In what follows, we approximate these two matrices using node-wise LASSO and thresholded HAC covariance estimation.
We further introduce some notation. For $\forall j \in [d]$, let $\mathbf{X}_j$ and $\mathbf{X}_{-j}$ respectively be the $j^{th}$ column of $\mathbf{X}$ and the sub-matrix of $\mathbf{X}$ with the $j^{th}$ column removed. Accordingly, we let
and it is easy to check that $|\pmb{\gamma}_{j}|_2 <\infty$ and $E(\mathbf{X}_{-j}^\top\pmb{\eta}_j) = \mathbf{0}$. We then run the following node-wise LASSO regression to debias the LASSO estimate and perform inference:
where $\{\widetilde{\omega}_j\}_{j=1}^{d}$ are a sequence of penalty terms at the order $\widetilde{\omega}_j \asymp \sqrt{\log (d)/(NT)}$. Finally, we propose the debiased LASSO estimator in Step 3 below.
\hrule
\hrule
Here, $\widehat{\pmb{\Omega}}_x$ serves as an asymptotic approximation to the inverse of $\mathbf{X}^\top\mathbf{X}/(NT)$. To study Step 3, we impose the following conditions.
Assumption 3.1 basically requires $\frac{s\, \log(d)}{\sqrt{NT}}\rightarrow 0$, while Assumption 3.2 imposes some conditions on two population matrices. Using them, we establish Theorem (ref), which applies to all elements of $\pmb{\beta}_0$, including those 0's, after choosing certain $\rho$.
Before we make the asymptotic normality in ((ref)) feasible for inferential purposes, we propose to estimate $\pmb{\Theta}$, and then develop a thresholded HAC covariance estimator.
\hrule
\hrule
The conditions on $a(\cdot )$ of Assumption (ref).1 are satisfied by a number of commonly used kernels, such as the Bartlett and Parzen kernels. The last two conditions of Assumption (ref).1 ensure the spectral norm consistency of our thresholded HAC covariance matrix estimator. On top of Assumption (ref), Assumption (ref).2 further requires an algebraic decay rate of the auto-covariance matrix, which is used to derive the order of a bias term involved in truncating the long-run covariance matrix. Assumption (ref).3 controls the order of elements in $\pmb{\Theta}_{\ell}$, and allows for the presence of many “small” but nonzero elements.
Note that the technical conditions (i.e., $E(e_{it} \mid \mathbf{X}) = 0$, $\lim\sup_{T\rightarrow \infty} s^2\sqrt{\ell \log (d)/T}<\infty$ and $\frac{d^3T\log T}{(T\ell \log d)^{q/4} }\to 0$) assumed in Theorem (ref) are only used to ensure that the estimation errors in $\{\widehat{e}_{it}\}$ are negligible. Theorem (ref) immediately infers that
Up to this point, our investigation about model (ref) is completed. Again, the above investigation considerably extends and enriches the relevant literature by allowing for the idiosyncratic errors to exhibit TSA, CSD, heteroskedasticity, as well as heavy-tailed behavior.
In this section, we generalize the above investigation for model (ref) to a class of widely used interactive fixed--effect models of the form:
which admits the following matrix representation:
where $\pmb{\Xi}_0 = \mathbf{F}_0\pmb{\Lambda}_0^\top$, $\mathbf{F}_0 = [\mathbf{f}_{01},\ldots,\mathbf{f}_{0T}]^\top$, $\pmb{\Lambda}_0 = [\pmb{\lambda}_{01},\ldots,\pmb{\lambda}_{0N}]^\top$, and $\mathbf{e}$ is defined accordingly. Note that $\pmb{\Xi}_0$ admits the singular value decomposition (SVD): $$\pmb{\Xi}_0 = \mathbf{U}_0\mathbf{D}_0\mathbf{V}_0^\top,$$ where $\mathbf{U}_0$ and $\mathbf{V}_0$ are $T\times T$ and $N\times N$ respectively. Let $\mathbf{U}_{0,[r]} \in \mathbb{R}^{T\times r} $ and $\mathbf{V}_{0,[r]} \in \mathbb{R}^{N\times r}$ be the sub-matrices of singular vectors associated with the largest $r$ singular values of $\pmb{\Xi}_0$, where $r$ is assumed to be fixed and finite throughout the rest of this paper.
The latent factor structure induces a nonconvex regression problem, and we propose to deal with it via a nuclear norm constraint (e.g., moon2018nuclear,belloni2023high). We again propose a multiple--step procedure to investigate (ref): (i) In Step 1, we use an $\ell_1$-nuclear norm penalized estimation method to obtain consistent estimator, which however suffers from substantial shrinkage bias (see Proposition (ref) for details); and (ii) In Step 2, we propose an iterative estimation method that removes the shrinkage bias in order to obtain robust estimation and inference.
\hrule
\hrule
To proceed, we need to introduce the following assumptions.
Assumption (ref).1 requires the idiosyncratic error matrix to have an operator norm of order $\max(\sqrt{N},\sqrt{T})$, and nests a class of high--dimensional MA($\infty$) processes as special cases (e.g., see Example (ref) of Appendix (ref) for details). Assumption (ref).2 is often referred to as the restricted strong convexity condition in the literature (e.g., negahban2011estimation,miao2023high). The first part of Assumption (ref).3 is commonly used in the relevant literature, and the second part of Assumption (ref).3 requires that the eigenvalues of $\pmb{\Sigma}_{\lambda}^{1/2}\pmb{\Sigma}_{f}\pmb{\Sigma}_{\lambda}^{1/2}$ are distinct in order to identify the corresponding eigenvectors.
The first condition of Assumption (ref).1 requires that the non-zero elements of $\pmb{\beta}_{0}$ are not too small. The second condition of Assumption (ref).1, together with the first condition, ensures the consistency of variable selection. The third condition of Assumption (ref).1 imposes an extra restriction on the diverging rate of the number of nonzero elements in $\pmb{\beta}_0$, which is used to validate the compatibility condition. Assumption (ref).2 requires the matrix to be positive definite, and is standard.
The following proposition establishes the consistency of model selection and an asymptotic distribution.
{
}
Finally, we deal with the two biases and establish valid inference in Step 3.
\hrule
\hrule
Before we prove the theoretical results in Appendix B, we evaluate the finite--sample performance of the proposed estimation and inferential procedure by simulated and real datasets.
A detailed numerical implementational procedure is given in Appendix (ref). Using it, we evaluate the results of Section (ref) by considering the following data generating process: $$ \textbf{DGP1}: \quad y_{it} = \alpha_i + \mathbf{x}_{it}^\top\pmb{\beta}_{0} + e_{it},\quad \mathbf{e}_t = \rho_e\mathbf{e}_{t-1} + \pmb{\Sigma}_{e}^{1/2}\pmb{\varepsilon}_{e,t}, $$ where $\alpha_i = \frac{1}{T}\sum_{t=1}^{T}(x_{1,it} + x_{2,it})$ is an individual fixed--effect, $\pmb{\Sigma}_{\varepsilon} = \{\delta_{\varepsilon_e}^{|i-j|}\}_{i,j\in [N]}$, $\pmb{\varepsilon}_{e,t}$ follows from an $N$-dimensional $t$-distribution with a degree freedom of 5, and $\rho_e, \delta_{\varepsilon_e}\in\{0.2,0.5\}$, corresponding to low and moderate dependence in the dynamics of error innovations.
DGP1 also allows for both the serial and cross--sectional correlations in $\mathbf{x}_{it}$ as follows: $\mathbf{x}_{l,t} = \rho_x + 0.2\mathbf{x}_{l,t-1} + \pmb{\Sigma}_{x}^{1/2}\pmb{\varepsilon}_{l,t} \ \ \text{for} \ \ l\in [d]$, where $\pmb{\Sigma}_{x} = \{0.2^{|i-j|}\}_{i,j\in [N]}$, $\rho_x \sim N(0,1)$, $\pmb{\varepsilon}_{l,t}$ follows from an $N$-dimensional $t$-distribution with a degree freedom of 5, and $\{\pmb{\varepsilon}_{l,t}\}_{l}$ is mutually independent for $l \in [d]$. Here, we let $d \in\{50,500\}$, $\beta_{0,j} = 0.2+0.1j$ for $j\in [5]$ and $\beta_{0,j} = 0$ for $j\ge 6$. When running regression, we deal with the fixed effects $\alpha_i$ by removing the time mean on both sides of DGP1. Therefore, we use the demeaned data to conduct the estimation procedure of Section (ref).
For DGP1, we consider two sets of sample sizes, which are $N, T \in (20,30,40)$ and $N,T \in (50,100,200,400)$. For each pair of $(N,T)$, we conduct $1000$ replications. In addition, we measure the accuracy of LASSO estimates by the root mean squared error (RMSE) $\text{RMSE}(\widehat{\pmb{\beta}}) \coloneqq \sqrt{\frac{1}{1000}\sum_{j=1}^{1000}|\widehat{\pmb{\beta}}^{(j)} - \pmb{\beta}_0|_F^2},$ where $\widehat{\pmb{\beta}}^{(j)}$ is the estimate of $\pmb{\beta}_0$ at the $j^{th}$ replication. In order to evaluate the finite sample performance of our estimation and inferential procedure, we calculate the empirical coverage rates (ECR) for the nonzero elements in $\pmb{\beta}_0$, i.e., $\beta_{0,1}$-$\beta_{0,5}$ based on Steps 3 & 4 of Section (ref). We then take the average across these elements for ease of presentation. For comparison, we also report the empirical coverage rates (denoted by ECR2) using the HAC estimator considered in babii2022machine, which is not robust in the presence of CSD. We also report the ratio of sign consistency (RSC) of the adaptive LASSO procedure (zou2006adaptive), i.e., the ratio of $\{\widehat{\pmb{\beta}}_{w}^{(j)} =_s \pmb{\beta}_0\}_{j=1}^{1000}$. Tables (ref)-(ref) report these results for the cases with $d = 50$ and $d = 500$ respectively.
In view of Tables (ref) and (ref), a few facts emerge. First, as expected, RMSE of the LASSO estimator decreases as both $N$ and $T$ increase. Second, as expected, RMSE increases if either $\rho_e$ or $\delta_{\varepsilon_e}$ increases. Third, when $\rho_e = 0.2$, ECR is very close to its nominal level even when $T=20$, although the distortion in ECR increases with the increase of $\rho_e$ and alters with $\delta_{\varepsilon_e}$ slightly. Fourth, when $\rho_e = 0.5$, ECR converges to its nominal level only with an increase in $T$. This is consistent with our theoretical prediction that the estimation error of long-rung covariance matrix is independent of $N$. Fifth, ECR2 is always below its nominal level and the distortion in ECR2 increases with the increase of $\delta_{\varepsilon_e}$. This is not surprising since the HAC estimator in babii2022machine is not robust to the presence of CSD and just includes a proportion of the asymptotic variance. Finally, the adaptive LASSO procedure can correctly identify the sparsity pattern as long as the sample size is not so small. When the sample size is relatively small, RSC converges to one rapidly as the sample size increases. In Tables (ref) and (ref), we find similar patterns. Interestingly, for $d = 500$ and $\rho_e = 0.2$, ECR tends to be larger than its nominal level when $T$ is relatively small. In addition, RMSE for $d = 500$ is slightly larger than that for $d = 50$, which is consistent with that RMSE should be proportion to $\sqrt{\log d}$ according to Lemma (ref).1.
We then evaluate the findings of Section (ref). Consider the following DGP: $$ \textbf{DGP2}: \quad y_{it} = \mathbf{x}_{it}^\top\pmb{\beta}_{0} + \pmb{\lambda}_{0i}^\top \mathbf{f}_{0t} + e_{it},\quad \mathbf{e}_t = \rho_e\mathbf{e}_{t-1} + \pmb{\Sigma}_{e}^{1/2}\pmb{\varepsilon}_{e,t}, $$ where $\pmb{\lambda}_{0i} = [x_{1,i1},x_{2,i1}]^\top$ and $\mathbf{f}_{0t} = [x_{3,1t},x_{4,1t}]^\top$. Here, $\mathbf{x}_{it}$, $\pmb{\beta}_0$ and $\{e_{it}\}$ are generated in exactly the same way as in DGP1. We compute the RMSE of the first and second stage estimators for DGP2, denoted them by $\text{RMSE1}$ and $\text{RMSE2}$ respectively. To evaluate the finite sample performance of our inference procedure, we compute the empirical coverage rates (ECR) of the non-zero elements in $\pmb{\beta}_0$, i.e., $\beta_{0,1}$-$\beta_{0,5}$. We then take the average across these elements. We also report the average shares of the relevant variables included (TPR, true positive rate) and the average shares of the irrelevant variables included (FPR, false positive rate) of the conservative LASSO procedure, as well as the exact estimation rate (EER) of the number of factors by using the singular value thresholding procedure. Tables (ref)-(ref) show the results of DGP2 for $d = 50$ and $d = 500$ respectively.
Tables (ref) and (ref) reveal some notable points. First, the first-stage $\ell_1$-nuclear norm penalized estimator has a larger RMSE due to the regularization, which slowly vanishes as the sample size increases. On the other hand, the second-stage has much smaller RMSE, which decreases as both $N$ and $T$ increase. Second, for $\rho_e = 0.5$, the finite sample coverage probabilities are smaller than their nominal level (95%) when $T$ is small, but are quite close to 95% as $T$ increases. Third, CSD and TSA in $\{e_{it}\}$ significantly affect the accuracy of our model selection procedure. For each pair of $(N,T)$, larger $\rho_e$ and $\delta_{\varepsilon_e}$ infer larger FPR. {Nevertheless, our procedure can pick up all the relevant variables when either $N$ or $T$ is larger than 40, and FPR converges to zero quickly when either $N$ or $T$ is large.} Finally, the proposed singular value thresholding procedure can correctly determine the number of factors {when the sample size is not too small}.
Tables (ref) and (ref) report the simulation results of DGP2 for $d = 500$. We find similar patterns. In an analogous way to DGP1, RMSE1 for $d = 500$ is slightly larger than that for $d = 50$, which is consistent with our theoretical prediction. Notably, increasing $d$ does not significantly affect the estimation accuracy of the number of factors.
Firm level characteristics which potentially can predict future stock returns have drawn considerable attention in the literature of asset pricing (see, KELLY2019501, ChenZimmermann2021, belloni2023high, and many references therein, for example). In this section, we aim to select the firm characteristics that provide incremental information about U.S. monthly stock returns and to estimate how selected characteristics affect expected returns. Importantly, we highlight the necessity of accounting for cross-sectional dependence and demonstrate its implications on inferring significant firm characteristics.
We collect the return data of different firms of S&P 500 from Center for Research in Security Prices (CRSP) (available at \url{www.crsp.org}), and match them with firm characteristics assembled by ChenZimmermann2021, which are available at \url{www.openassetpricing.com}. These data are recorded monthly, and we pay attention to the period from Jan of 1990 to Dec of 2023 which in total gives 408 time periods. After matching “permno" code in both datasets and removing firm characteristics with missing values, we end up with 161 firms (i.e., $N=161$) and 60 characteristics (i.e., $d=60$). We present the variable names with their means and standard deviations in Table (ref), and refer interested readers to ChenZimmermann2021 for detailed definitions of these variables. As these variables are measured in different units and have significant differences in terms of their standard deviations, we normalize each variable (including stock return) to ensure mean 0 and standard deviation 1 before regression.
We begin by presenting the results using model (ref) with fixed effects as in the simulation studies. To run estimation, we always regress the return of firm $i$ at time $t$ ($y_{it}$) on the firm level characteristics of firm $i$ at time period $t-1$ ($\mathbf{x}_{i,t-1}$). The estimation procedure is identical to those presented in Section (ref), so we do not repeat the details here. To show the presence of the weak CSD of the estimation residuals, we conduct the CD test\footnote{The asymptotic distribution of the CD test follows the standard normal distribution, so at the 5% significance level, the critical values are $\pm 1.96$. We refer interested readers to Pesaran2004 for more details.} of Pesaran2004 on the residuals, and obtain the test statistic 61.32 which is a sign of the existence of weak CSD. Additionally, we run Jarque-Bera test based on the estimation residuals for each $i$, and find that we reject the null for all $i$'s, which infers that no individual has normally distributed residuals. Therefore, it shows the necessity of accounting for non-Gaussian error process in theory.
The weighted Lasso procedure identifies five significant determinants in future stock returns, which are Earnings announcement return, Pastor-Stambaugh liquidity beta, Idiosyncratic risk (3 factor), Maximum return over month, Realized volatility, respectively. In contrast, by examining the estimated confidence intervals for each element in the debiased Lasso estimates, we find that the predictor `Maximum return over month' is insignificant with a slope coefficient of -0.132 and a 95% confidence interval of (-0.278, 0.014). It is not surprising that the predictor `Maximum return over month' survives in the weighted Lasso selection procedure as the value `-0.1321' is quite large in this study. In addition, the estimated confidence intervals indicate five other (marginally) significant predictors, which may be due to the well-known multiple-hypothesis-testing problem in this case. Note that even when there is no truly significant firm characteristics, we expect to identify about three significant predictors ($60\times 5\%$) due to pure sampling variation. Then, four firm level characteristics are selected, and are reported in Table (ref) above.
The above selection result shares certain similarity with belloni2023high. For example, both studies acknowledge the importance of liquidity beta. While they focus on the selection at different quantile, we emphasize the importance of accounting for dependence of error terms. For the sake of comparison, we also examine the estimated confidence intervals using the HAC estimator considered in babii2022machine, which is not robust to the presence of CSD. These estimated confidence intervals indicate that there are 21 significant firm characteristics, which is unusually large. We do expect that these estimated confidence intervals are quite narrow as the HAC estimator considered in babii2022machine only just include a proportion of the asymptotic variance with the presence of CSD. This point can be also seen in the constructed confidence intervals reported in Table (ref). Overall, the above estimation and selection results highlight the necessity of accounting for cross-sectional dependence when selecting significant firm characteristics.
{
}
We next consider the model (ref) with unobservable common factors, of which the latter one is also considered in KELLY2019501 but under the finite dimensional framework. We also conduct the CD test to detect the presence of the weak CSD of the estimation residuals, and obtain the test statistic of $-11.78$. For both models, we end up with $T=407$. This result indicates the presence of weak CSD even after accounting for unobservable common factors (strong CSD).
{
}
Again, since our weighted LASSO procedure may select irrelevant variables in finite sample, we further eliminate irrelevant firm level characteristics based on the constructed confidence intervals. Then, four firm level characteristics are selected, and are reported in Table (ref) above. Compared to the estimates associated with the selected variables using the fixed effects model, interestingly, we find that after accounting for the common factors, the estimates tend to have narrower confidence intervals overall. This point is important for detecting the significant firm characteristics because the literature (e.g., rapach2013) shows that certain predictors do provide useful signals for forecasting stock returns but these signals are of small magnitudes and are hidden by the large uncertainty of error innovations. Instead, the large uncertainty of error innovations can be captured by using the method of interactive fixed effects. It then shows the necessity of including the factor structure practically, and so demonstrates the practical relevance of Section (ref).
In this paper, we propose a robust inferential procedure for the proposed HD panel data models. Specifically, (i) we pay attention to non-Gaussian, serially and cross-sectionally correlated and heteroskedastic error processes; (ii) we develop an estimation method for high-dimensional long-run covariance matrix using a thresholded estimator; and (iii) we allow for the number of regressors to grow faster than the sample size. In order to establish the corresponding theory, we derive two Nagaev-types of sharp concentration inequalities, one for a partial sum and the other for a quadratic form, subject to a set of easily verifiable conditions. Leveraging these two inequalities, we develop a non-asymptotic bound for the LASSO estimator, establish an asymptotic normality via the nodewise LASSO regression, and derive a sharp convergence rate for the thresholded HAC estimator.
We believe that our study provides the relevant literature with a complete toolkit for conducting estimation and inference for the parameters of interest within a class of HD panel data settings. We also demonstrate the practical relevance of these estimation and inferential methods by investigating a class of HD panel data models associated with interactive effects. Moreover, we conduct extensive numerical studies using simulated and real data examples.
{ {2.pt plus 0ex} }