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.
17,913 characters · 4 sections · 30 citation commands
Factor-Augmented Machine Learning Panel Regressions
\noindentKeywords: factor-augmented panel regressions; sparse-group LASSO; MIDAS; mixed-frequency data; high-dimensional forecasting; cross-sectional dependence.
\noindentJEL classification: C23; C32; C53; C55; C58.
\etocsettocdepth.toc{none}
Modern nowcasting and forecasting problems often rely on high-dimensional panel datasets with many cross-sectional units, a large number of predictors, and mixed-frequency data arriving in real time. Such panels typically exhibit strong cross-sectional dependence, reflecting common shocks in macroeconomic and financial time series. At the same time, the predictive component is often both sparse and dense: a few observed predictors may be important, after controlling for common factors generating broad comovement across units and variables.
This paper develops asymptotic theory for a factor-augmented sparse-group LASSO framework for nowcasting and forecasting with mixed-frequency panels. The previous literature on mixed-frequency regressions includes the low-dimensional MIDAS regressions introduced in ghysels2006predicting, ghysels2007midas, and andreou2010regression.
There is also a large literature on high-dimensional machine learning methods for i.i.d. data, e.g. bickel2009simultaneous, belloni2014inference, and cai2022sparse.
babii2022machine proposed to use the sparse-group LASSO of simon2013sparse for high-dimensional time series and developed the asymptotic theory. mogliani2021bayesian proposed a Bayesian version of the group LASSO for MIDAS regressions. Existing high-dimensional panel regression methods, such as babii2022machinepanel, allow for high-dimensional covariates under approximate sparsity but do not exploit factor structure to capture cross-sectional dependence in the errors. Related sparse-plus-dense and factor-based panel approaches, see hansen2016factor, ruecker2022estimation, and fan2022bridging, do not cover the mixed-frequency nowcasting problem with MIDAS aggregation and sparse-group regularization. The closest work to this paper is beyhum2024factor, which studies factor-augmented sparse MIDAS regressions for time series, but not panel data regressions.
\paragraph*{Notation}
For an integer $N\in\mathbb{N}$, let $[N]=\{1,\ldots,N\}$. For a finite set $C$, let $|C|$ denote its cardinality. For a vector $v\in\mathbb{R}^p$, let $|v|_q=(\sum_{j=1}^p|v_j|^q)^{1/q}$ denote the $\ell_q$-norm with $|v|_\infty = \max_{j\in[p]}|v_j|$ and $|v|_0=\sum_{j=1}^p\mathbf{1}_{v_j\ne 0}$. A group is defined as a set of indices $G\subset[p]$. Let $v_G\in\mathbb{R}^p$ be a vector such that $(v_G)_i=v_i$ if $i\in G$ and $(v_G)_i=0$ otherwise. For vectors $v,w\in\mathbb{R}^n$, let $\langle v,w\rangle_n=n^{-1}\langle v,w\rangle$, where $\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product, and let $\|\cdot\|_n=\sqrt{\langle\cdot,\cdot\rangle_n}$ be the associated norm. For a matrix $A\in\mathbb{R}^{n \times m}$, let $\|A\|_{\rm op}=\sup_{|v|_2=1}|Av|_2$. For a collection of groups $\mathcal{G}=\{G_1,\dots,G_K\}$, let $\|v\|_{2,1}=\sum_{k=1}^K|v_{G_k}|_2$ and $v_{\mathcal{G}}$ be the vector such that $(v_{\mathcal{G}})_i=v_i$ if $i\in\mathcal{G}$ and $(v_{\mathcal{G}})_i=0$ otherwise. The quantity $n_1\vee n_2$ is the maximum of $n_1$ and $n_2$. For a matrix $A$, let $\sigma_r(A),\sigma_{\min}(A)$ and $\sigma_{\max}(A)$ denote the $r^{\rm th}$ largest, smallest, and largest singular values of $A$, respectively.
Consider the following factor-augmented predictive regression:
where $q_{i,t}\in\mathbb{R}^p$ includes all observed covariates, $f_{i,t}\in\mathbb{R}^R$ are latent factors, and $a_{i,t}\in\mathbb{R}$ is a misspecification error due to approximate sparsity or misspecified factor structure. The latent factors $f_{i,t}$ and the misspecification error $a_{i,t}$ can generate the cross-sectional dependence in errors.\footnote{This covers, for example, the interactive-effects specification $f_{i,t}=\Lambda_i g_t$, where $g_t\in\mathbb{R}^r$ is a common factor and $\Lambda_i\in\mathbb{R}^{R\times r}$ is a unit-specific loading matrix. Cross-sectional dependence then arises through the common factor; see also babii2025tensor for a more refined tensor factor model analysis.} Therefore, the model represents a flexible specification with approximately sparse signals, dense components, and cross-sectional dependence.
Covariates may be measured at higher frequencies than the outcome $y_{i,t+1}$. In this case we will transform the high-frequency data into Equation (ref) using the MIDAS approach of babii2022machine. Consider a panel of $K_x$ covariates for $N$ units observed over $T$ low-frequency periods, with $m_x$ high-frequency observations for each low-frequency period:\footnote{For simplicity of presentation, we assume that the number of high-frequency observations is the same across all low-frequency periods. This assumption can be easily relaxed, and the observations can also be allowed to extend forward into the nowcasting period.} \[ \left\{x_{i,t-(j-1)/m_x,k}^H:\ i\in[N],\ t\in[T],\ j\in[m_x],\ k\in[K_x]\right\}. \] Suppose this panel follows the approximate factor model:
where $b_k\in\mathbb{R}^{K_f}$ is the loading vector, $f_{i,t-(j-1)/m_x}^H\in\mathbb{R}^{K_f}$ is the high-frequency factor, and $u_{i,t-(j-1)/m_x,k}^H$ is the idiosyncratic component. A second panel of $K_z$ high-frequency covariates \[ \left\{z_{i,t-(j-1)/m_z,k}^H:\ i\in[N],\ t\in[T],\ j\in[m_z],\ k\in[K_z]\right\} \] need not follow a factor model and will enter Equation (ref) only through the approximately sparse coefficients.
For each unit $i$ and covariate $k$, collect the high-frequency observations in matrices $X_{i,k}^H :=(x_{i,t-(j-1)/m_x,k}^H)_{t\in[T],\,j\in[m_x]}$ and $Z_{i,k}^H:=(z_{i,t-(j-1)/m_z,k}^H)_{t\in[T],\,j\in[m_z]}$. These high-frequency observations can be mapped with MIDAS weights into the regression model in Equation (ref) as follows. Let $w_l:[0,1]\to\mathbb{R}$ for $l\in[L]$ be a dictionary of functions used to approximate the MIDAS weights. Define the weighting matrices $W^x:=\left(w_l((j-1)/m_x)/m_x\right)_{j\in[m_x],l\in[L]}$ and $W^z:=\left(w_l((j-1)/m_z)/m_z\right)_{j\in[m_z],l\in[L]}.$ Then $X_{i,k}^H W^x\in\mathbb{R}^{T\times L}$ and $Z_{i,k}^H W^z\in\mathbb{R}^{T\times L}$ are the MIDAS-weighted versions of high-frequency observations. Define
where $p_x=LK_x$ and $p_z=LK_z$. Let $x_{i,t}\in\mathbb{R}^{p_x}$ and $z_{i,t}\in\mathbb{R}^{p_z}$ be the $t^{\text{th}}$ rows of $\mathbf{x}_i$ and $\mathbf{z}_i$ transposed, respectively. The aggregated covariates are collected in a vector $q_{i,t}:=(1,z_{i,t}^\top,x_{i,t}^\top)^\top\in\mathbb{R}^p$, where $p=p_z+p_x+1$, and we assume that the initial observations of lagged variables are available. Then $q_{i,t}$ corresponds to the vector of covariates in Equation (ref).
For estimation purposes, using matrix notation, define $\mathbf{q}_i:=(q_{i,1},\dots,q_{i,T})^\top\in\mathbb{R}^{T\times p}$, $\mathbf{y}_i:=(y_{i,2},\dots,y_{i,T+1})^\top\in\mathbb{R}^{T}$, and stack $\mathbf{y}:=(\mathbf{y}_1^\top,\dots,\mathbf{y}_N^\top)^\top\in\mathbb{R}^{NT}$, $\mathbf{Q}:=(\mathbf{q}_1^\top,\dots,\mathbf{q}_N^\top)^\top\in\mathbb{R}^{NT\times p}$, and $\mathbf{X}:=(\mathbf{x}_1^\top,\dots,\mathbf{x}_N^\top)^\top\in\mathbb{R}^{NT\times p_x}$. Let also $\mathbf{A}\in\mathbb{R}^{NT}$ and $\mathbf{E}\in\mathbb{R}^{NT}$ be the stacked versions of $a_{i,t}$ and $\varepsilon_{i,t+1}$. Using matrix notation, Equation (ref) becomes
The factor structure of covariates, see Equation (ref), implies that the dense block of covariates can be written as
where $\mathbf{F}\in\mathbb{R}^{NT\times R}$ are the (MIDAS-aggregated) factors with $R=LK_f$, $B\in\mathbb{R}^{p_x\times R}$ are the factor loadings, and $\mathbf{U}\in\mathbb{R}^{NT\times p_x}$ are the idiosyncratic components. The MIDAS aggregation introduces a group structure in the regression model, where each original predictor contributes a group of $L$ dictionary coefficients; see babii2022machine. The group structure with $K$ groups is described as a partition $\{G_1,\dots,G_K\}$ of $[p]=\{1,\dots,p\}$.
We extract factors from $\mathbf{X}$ using PCA. Let the columns of $\hat{\mathbf{F}}/\sqrt{NT}$ be the eigenvectors associated with the largest $R$ eigenvalues of $\mathbf{X}\mathbf{X}^\top$, so that $\frac{1}{NT}\hat{\mathbf{F}}^\top\hat{\mathbf{F}}=I_R$. The factor-augmented sg-LASSO estimator is
where $|d|_1 = \sum_{j=1}^p |d_j|$ and $\|d\|_{2,1} = \sum_{k=1}^K |d_{G_k}|_2$ are the $\ell_1$ and $\ell_{2,1}$ norms, respectively, and $\lambda_1,\lambda_2\geq0$ are tuning parameters. The factor-augmented sg-LASSO thus captures both the dense component of the covariates through the factors and the approximately sparse component through the penalized regression on $\mathbf{Q}$. The $\ell_1$ LASSO penalty encourages sparsity at the coordinate level, while the $\ell_{2,1}$ group LASSO penalty encourages sparsity at the group level.
The estimator in Equation (ref) combines a dense component, represented by the factors estimated from $\mathbf{X}$, with a sparse component, represented by the penalized regression on $\mathbf{Q}$; see also beyhum2024factor. We consider an asymptotic regime in which $N,T\to\infty$, the number of covariates $p\to\infty$ may grow with $N$ and $T$, and the number of factors $R$ remains fixed. The fixed-$R$ assumption can also be relaxed; see babii2025tensor, beyhum2022factor, and freeman2023linear.
Let $S=\{j\in[p]:\delta_j\ne0\}$ be the active coordinates in the true parameter $\delta\in\mathbb{R}^p$ and let $\mathcal A=\{k\in[K]:\delta_{G_k}\ne0\}$ be the active group indices. Put $s=|S|$, $m=|\mathcal A|$, and for an integer $\ell\geq1$, let
be the total number of elements in the largest $\ell$ groups. Let $P_{\hat{\mathbf{F}}}=(NT)^{-1}\hat{\mathbf{F}}\hat{\mathbf{F}}^\top$, $M_{\hat{\mathbf{F}}}=I_{NT}-P_{\hat{\mathbf{F}}}$, and $\tilde{\mathbf{Q}}=M_{\hat{\mathbf{F}}}\mathbf{Q}$.
Assumption (ref) requires that regression errors are sub-Gaussian conditionally on the covariates. It can be relaxed in several different ways, e.g. winsorizing the data or robustifying the least-squares objective function with Huber or median-of-means losses; see lugosi2019mean for more details. Assumption (ref) is a restricted eigenvalue condition and allows for singularity of the design matrix. Assumption (ref) requires the oracle rates for tuning parameters. It is also known that the LASSO estimator with tuning parameters selected by cross-validation achieves similar results as the oracle LASSO estimator; see chetverikov2021cross.
The following result holds:
Theorem (ref) states the oracle inequality for the factor-augmented sg-LASSO estimator. The proof of this and other results can be found in the Appendix. Next, we deduce the convergence rate of the factor-augmented sg-LASSO estimator under additional assumptions. Let $u_{i,t}\in\mathbb{R}^{p_x}$ be a column vector corresponding to the row for observation $(i,t)$ in $\mathbf{U}$.
Assumption (ref) imposes a set of restrictions needed to characterize the estimation error in factors; see bai2003inferential. Some of the more high-level conditions in this assumption can also be reduced to a set of more primitive conditions; see e.g. beyhum2024factor and references therein.