Factor Models with Sparse VAR Idiosyncratic Components

Factor Models with Sparse VAR Idiosyncratic Components

Introduction

In the past twenty years, factor models have emerged as a major tool for the analysis and forecast of high-dimensional time series. Such models are characterized by the assumed existence of a specific decomposition of the high-dimensional vector of time series into two mutually orthogonal unobserved components. The common component, driven by a finite number of (possibly dynamic) factors, represents the comovements among the series and it is of reduced rank. The weakly cross-correlated idiosyncratic component (cf. generalized factor model) represents individual features of the series. \footnote{Conversely to the generalized factor model, the exact factor model assumes no cross-sectional dependence in the idiosyncratics, thus working with the assumption of a diagonal cross-autocovariance matrix.} The ramification of the literature on factor models is mostly due to the way this decomposition is characterized. Traditionally, the common component is identified as those eigenvalues of the covariance matrix which diverge while the idiosyncratics are those which stay bounded. forni2000generalized and forni2015dynamic assume the explosive eigenvalues to be also dynamic, reflecting both a contemporaneous and lagged effect of the common components on the series. stock2002forecasting,bai2002determining,bai2019rank assume explosive static eigenvalues (i.e. only contemporaneous effects of the factors) along with a finite-dimensional factor space. One thing is certain: the two components are radically different objects that need to be treated differently to better capture their respective dependence structures. Assuming away any cross-sectional dependence of the idiosyncratic components can be misleading. Some consider them as simple univariate autoregressive processes forni2005generalized or even as white noise processes and drop them when producing forecasts lam2012factor. This however neglects the predictive power that idiosyncratic components can have, thus resulting in less accurate forecasts. In this paper we lean forward towards a reduced rank plus sparse characterization of the factor model decomposition by assuming the idiosyncratic component to follow a high-dimensional vector autoregressive (VAR) model. This allows cross-sectional and time dependence in the idiosyncratic term. In the first estimation step, we employ principal component analysis (PCA) to estimate the factors. In a second step, high-dimensional penalized VAR through the (adaptive) lasso is used in order to estimate the idiosyncratic components. By so doing we combine a “dense" modeling approach for the factors with a “sparse" modeling approach for the idiosyncratic component.\footnote{The terms dense and \textit{sparse} are used to distinguish estimation approaches that require or not the structural assumption of a sparse coefficient vector to perform some dimensionality reduction giannone2017economic. PCA is therefore a dense approach while lasso is a sparse one.} Thus, allowing for a more refined disentangling of the dependence structure of the two components. We show the consistent estimation of both the sparse VAR model driving the idiosyncratic components and the factors, as both the cross-sectional and time dimensions grow large. When estimating the sparse VAR for the idiosyncratic component, a naive approach would be to simply plug in the standard rates derived for the factor estimation. This however leads to a suboptimal rate. Instead, an important contribution of our work is deriving detailed expressions of the occurring errors. This enables us to obtain tighter rates for the second step and also employ a semi-parametric estimator for the inverse of the spectral density matrix. We discuss the implications of our proposed framework for forecasting, factor-augmented regression, bootstrap of factor models, and the estimation of time series dependence networks. We also propose a joint information criterion that combines the approach of bai2002determining with an extra penalty allowing for simultaneous lag-length estimation of the VAR model. The benefit of our proposed procedure is confirmed through extensive simulations where different levels of sparsity, number of factors, lag-length of the VARs, and idiosyncratic covariance matrix are considered for different sample sizes and dimensions. We also compare our combined procedure with the standard high-dimensional forecasting methods which fully rely on either a sparse or a dense procedure. There already exists applications in the literature combining (dynamic) factor models with sparse vector autoregressive models, see e.g., barigozzi2017network,barigozzi2019nets and more recently barigozzi2022fnets. However, barigozzi2017network,barigozzi2019nets do not present theoretical results about the combined approach and the framework considered in barigozzi2022fnets differs in important aspects from the one considered here, see Section (ref) and the discussion after Assumption (ref) for details. In the non-dynamic idiosyncratics set-up, kneip2011factor,fan2020factor combine factors with regularized models. Since regularized methods such as the lasso have difficulties with strongly correlated regressors, especially in the context of model selection, they aim to decorrelate the regressors by adjusting for the factors. Furthermore, fan2021bridging provide hypothesis tests to check whether after removing factors (as well as trends in a first step) the regressors possess some pre-defined weakly correlated structure or not. fan2020factor,fan2021bridging allow for time-dependent regressors, however, they do not consider nor allow that the idiosyncratic part follows a sparse vector autoregressive model where the cross-sectional sparsity can grow with the sample size. In the context of high-dimensional VAR models, another approach is to consider the slope matrices as a combination of a low-rank matrix and a sparse matrix as done in basu2019low. The low-rank part takes here a similar role as the common component of the factor model and it is estimated by nuclear-norm regularization. In the context of high-dimensional VAR models with strong cross-sectional correlated noise, a combination of low-rank plus sparse has been also explored by lin2019approximate and more recently in miao2022high. However, these approaches differ from ours in terms of model and estimation approach. A more detailed discussion can be found in Remark (ref) in Section (ref). The remainder of the paper is organized as follows: Section (ref) introduces the factor model with sparse VAR idiosyncratic components and reports few standard assumptions defining its behavior. Section (ref) is devoted to describing the two-step procedure used to estimate the factor model with sparse VAR idiosyncratic components and prove its consistency. Theorem (ref) derives a representation of the idiosyncratic components estimation error while Theorem (ref) is the main result establishing bounds for the estimation error for the second step of the estimation procedure i.e., for the lasso on the sample estimates of the idiosyncratic component. The same two-step procedure with mild additional assumptions can also be employed in estimating the spectral density of the process and Theorem (ref) derives the relative estimation error bounds. { Section (ref) discusses the implications of our model for forecasting, factor augmented regression, bootstrapping factor models, and the estimation of time series dependence networks. The latter point is illustrated by estimating networks based on the partial coherence for the FRED-MD data set.} Section (ref) considers the problems of: estimating the number of factors, determining the lag-length in the VAR, and tuning the penalty parameter for the lasso. Section (ref) reports simulation results for our proposed method under different VAR data generating processes in terms of design and sparsity. Finally, Section (ref) concludes. \bigbreak A few words on notation. Throughout the paper we use boldface characters to indicate vectors and boldface capital characters for matrices. For any $n$-dimensional vector $\boldsymbol{x}$, we let ${\left\lVert\boldsymbol{x}\right\rVert}_p = \left(\sum_{i=1}^n |x_i|^p \right)^{1/p}$ denote the $\ell_p$-norm and $\bm e_j=(0,\ldots,0,1,0,\ldots, 0)^\top$ denotes a unit vector of appropriate dimension with the one appearing in the $j$th position. Furthermore, for a $r\times s$ matrix $\BS A=(a_{i,j})_{i=1,\ldots,r, j=1,\ldots,s}$, $\|\BS A\|_1=\max_{1\leq j\leq s}\sum_{i=1}^r|a_{i,j}|=\max_j \| \BS A \bm e_j\|_1$, $\|\bm A\|_\infty=\max_{1\leq i\leq r}\sum_{j=1}^s|a_{i,j}|=\max_{i} \| \bm e_i^\top \BS A\|_1$ and $\|\BS A\|_{\max}=\max_{i,j} |\bm e_i^\top \BS A \bm e_j|$. $\bm A^i$ denotes the $i$th matrix power of $\bm A$ and $\bm A^{(i)}$ refers to the $i$th element of a sequence of matrices. We denote the largest absolute eigenvalue of a square matrix $\BS A$ by $\sigma_{\max}(\BS A)$ and let $\|\BS A \|_2^2=\sigma_{\max}(\BS A \BS A^\top)$. We denote the smallest eigenvalue of a matrix $\bm A$ by $\sigma_{\min}(A)$. For any index set $S \subseteq \{1, \ldots, n\}$, let $\boldsymbol{x}_{S}$ denote the sub-vector of $\boldsymbol{x}_t$ containing only those elements $x_i$ such that $i \in S$. $\|\BS x\|_0$ denotes the number of non-zero elements of $\BS x$.

The Model

We work {with a generalized factor model} where both factors and idiosyncratic components are allowed to be {(second order)} stationary stochastic processes and the loadings are static. {To elaborate,} let $\boldsymbol{x}_t=(x_{1,t},\ldots,x_{N,t})^\top$, $t=1,\ldots, T$, be a $N\times T$ rectangular data array representing a finite realization of an underlying real-valued stochastic process $\{x_{i,t}\}.$ Assume that for each $t$, $\boldsymbol{x}_t$ can be decomposed into a sum of a common component $\boldsymbol{\chi}_t=(\chi_{1,t}, \dots, \chi_{N,t})^{\top}$ and an idiosyncratic component $\boldsymbol{\xi}_t=(\xi_{1,t}, \dots, \xi_{N,t})^{\top}$, both latent {and mutually orthogonal at all leads and lags}. Then, the factor model decomposition takes the following usual form

{The common component has reduced rank i.e.,} $\BS \chi_{t}$ is driven linearly by an $r$-dimensional vector of common factors $\boldsymbol{f}_t=(f_{1,t},\ldots,f_{r,t})^{\top}$, where $r$ is considered as fixed as both the cross sectional dimension $N$ and the time series dimension $T$ grow large and $r\ll N$. The common components $\chi_{i,t}$ can {then be represented by} the following linear combination

where $\ell_{i,s}$ are denoted as loadings. Note that $\chi_{i,t}$ is uniquely defined. But since for {any} rotation matrix $\bm H$, $\chi_{i,t}=\boldsymbol{\Lambda}^{\top}_i \bm H \bm H^{-1}\boldsymbol{f}_t$ is a valid linear combination as well, $\boldsymbol{\Lambda}^{\top}_i, \boldsymbol{f}_t$ are only identified up to some arbitrary rotation.

We do not assume that the factors $\boldsymbol{f}_t$ nor the idiosyncratic component $\boldsymbol{\xi}_t$ are independent and identically distributed but we allow them to be stationary stochastic processes. We consider that the factors are given by a one-sided linear process, see Assumption (ref) below. This includes the cases that the factors are driven by a stable vector autoregressive (VAR) model. Additionally, the idiosyncratic component $\boldsymbol{\xi}_t$ {is allowed to be weakly cross-correlated and we consider this} to follow a sparse VAR model of order $p$ as

for $\bm v_t$ being a white noise process, $\bm A^{(j)}$ the sparse slope matrices and $\bm B^{(j)}$ the moving average matrices of the vector moving average (VMA)$(\infty)$-representation of the VAR$(p)$ model; see Assumption (ref) and (ref) below for details on sparsity and moment conditions. This includes as special case that the idiosyncratic components are driven by individual univariate autoregressive processes or even i.i.d..

{Let us note that even if the factors are dynamic, the relationship between $\boldsymbol{x}_t$ and $\boldsymbol{f}_t$ is assumed here to be static. This type of factor model is widely used in practice and it differs from the framework of forni2000generalized which assumes a pervasive dynamics of the common factors where $\boldsymbol{x}_t$ is set to also depend on $\boldsymbol{f}_t$ with lags in time. In several cases, though it is possible to transform a dynamic relationship into a stacked static relationship, see among others Section 2.1.2 in stock2016dynamic. Especially, if the assumption of a finite-dimensional span of the common component is used, then one can cast a dynamic representation into a static one bai2007determining. Alternatively, block-VAR filtering of $\bm x_t$ as proposed in forni2015dynamic also allows to turn lagged loadings into static ones.}

In the following Assumptions (ref), (ref), and (ref), the sparsity and stability conditions, the factors, moment conditions, and loadings are further specified. \bigbreak

Note that Assumption (ref),(i) is quite general as the sparsity is row-wise and {is allowed} to grow with the sample size. Let us emphasize that the weaker assumption of approximate sparsity instead of exact sparsity (i.e., $q=0$), is used throughout. In the context of forecasting, the assumption of a stable and row-wise sparse VAR model is standard in the literature of sparse VAR models, see among others kock2015oracle,han2015direct,masini2019regularized. When the focus is on estimating the dependency structure, e.g., spectral density matrices, additional column-wise sparsity seems unavoidable, see krampe2020statistical for a discussion of different sparsity concepts for VAR models. Here, we only require additional column-wise sparsity in Section (ref), where spectral density estimation is discussed. Row-wise sparsity (with or without additional column sparsity) includes as special case the univariate autoregressive model for each idiosyncratic component. The latter is generally not allowed, if the sparsity condition is specified for the entire matrix, e.g., $\sum_{s=1}^p \sum_{i,j=1}^N |\BS A_{i,j}^{(s)}|^q \leq k$, see Section 2 in krampe2020statistical for further discussion.

{The common and idiosyncratic components are identified based on the diverging behavior of the eigenvalues of the covariance matrix. This implies restrictions with respect to the matrix-norm $\|\cdot\|_2$ but not with respect to $\|\cdot\|_\infty$. Hence, as $\bm \xi$ is an idiosyncratic component, $\|\BS \Gamma_\xi(0)\|_2$ is bounded but $\|\BS \Gamma_\xi(0)\|_\infty$ can still grow with dimension. Consequently, when modeling the idiosyncratic component with a VAR model such behavior should not be ruled out by over-restricting the VAR slope matrices. Since row-wise sparsity of the slope matrices allows that $\|\boldsymbol{\frak{A}}\|_\infty$ can grow with sparsity parameter $k$, Assumption (ref) allows growth in $\|\BS \Gamma_\xi(0)\|_\infty$ as it is specified by the (possibly growing) parameter $k_\xi$. }

{If one would work within the framework of forni2017dynamic, their assumptions on the serial and cross-sectional dependence of the idiosyncratic terms are more restrictive. To be specific, Assumption 4 in forni2017dynamic impose $\|\cdot\|_1$- and $\|\cdot\|_\infty$-boundedness on the idiosyncratic coefficient matrices. This would imply that $k_\xi$ is bounded and in our case would also imply the slope VAR matrices to be bounded in matrix-norm $\|\cdot\|_\infty$. Hence, this would restrict the sparsity of the VAR slope matrices to be (more or less) fixed which would be less general and not desirable. Let us mention that this framework of dynamic factor models with boundedness conditions on the idiosyncratic coefficient matrices is considered in barigozzi2022fnets. An extension of the forni2017dynamic assumption to allow for growing sparsity is beyond the scope of this paper and is left for future research.}

In the established literature on factor estimation, a common assumption is to restrict the growth of the linear dependence of the idiosyncratic component. For instance, Assumption 3C in bai2003 states that the absolute sum of all covariances of the idiosyncratic component grow with order $N$. Here, we quantify the linear dependence of the idiosyncratic component by the condition $\|\Gamma_\xi(0)\|_\infty=\max_i \sum_{j=1}^N |\mbox{Cov}(x_{i,t},x_{j,t})|\leq k_\xi$ and the object $k_\xi$. Since $\|\Gamma_\xi(0)\|_\infty\leq \sqrt{N} \|\Gamma_\xi(0)\|_2$, $\sqrt{N}$ is an upper bound for the growth rate of $k_\xi$. As discussed previously, a bounded $k_\xi$ can be too restrictive hence we consider that $k_\xi$ grows moderately. We do not specify here a rate for $k_\xi$ but a rate smaller than $\sqrt{N}$ seems most realistic and would be more in line with established assumptions in the factor literature. The reason for this is that a growth rate of $\sqrt{N}$ would allow the absolute sum of all covariances could grow with a rate $N^{3/2}$. This would violate, for instance, the previously mentioned Assumption 3C in bai2003. Note further that if the maximal absolute row sum of the covariance of the idiosyncratic component is growing way faster than the average row sum, we may end up in the context of weak factors, see among others onatski2012asymptotics. Note that using the VMA$(\infty)$-representation, Assumption (ref) gives the upper bound $k^2 \|\BS \Sigma_v\|_\infty$ for $k_\xi$, where $\BS \Sigma_v=\mbox{Var}(\BS v_t)$ is the variance matrix of the residuals of the idiosyncratic component.

Assumption (ref) and (ref) imply that $\{\bm \xi_t\}$ is stationary and let the autocovariance function be given by $\boldsymbol{\Gamma}_{\xi}(s-t)=\mbox{Cov}(\boldsymbol{\xi}_{s},\boldsymbol{\xi}_t)$. Furthermore, Assumption (ref) implies that the factors are also a stationary process such that $\{\boldsymbol{x}_t\}$ itself is indeed stationary. In order to quantify the dependence of stochastic processes, we use the concept of functional dependence, see wu2005nonlinear. Since this is only necessary for the proofs, we do not introduce the notation here and refer to Remark (ref) in the appendix.

The moment condition in Assumption (ref) refers to the situation in which only a finite number of moments, here $\zeta$, are finite. Hence, we do not assume sub-Gaussian processes or similar, which is often assumed for sparse VAR processes, see basu2015,kock2015oracle,han2015direct. For sub-Gaussian processes, the polynomial terms depending on $\zeta$ would vanish, which would result in tighter error bounds obtained later on. The reason for only assuming $\zeta$ finite moments is to be more in line with the classical factor literature, see among others bai2003,stock2002forecasting,forni2000generalized,forni2017dynamic. E.g., bai2003 derived inferential results for factor models under $8$th finite moments of the idiosyncratic part and $4$th finite moments of the factors. Note that the filter in Assumption (ref) can be the one-sided representation of a stable VARMA model as in Assumption 2 in forni2017dynamic. Assumption (ref) is a standard assumption in the context of strong factor models, see stock2002forecasting,bai2003. It implies that each of the factors provides a non-negligible contribution to the variance of each component of $\{\BS x_t\}$. We would like to point out here that the time and cross-sectional dependence of the idiosyncratic component is only limited by assuming that it follows a sparse VAR model. Furthermore, it is not clear if assuming a sparse VAR model for the idiosyncratic part is a special case of the assumptions to time and cross-section dependence in the factor literature, see among others Assumption C in bai2003. The reason for this is that the sparsity is not fixed but it can grow with the sample size. Nevertheless, the error bounds obtained later on requires that the sparsity cannot grow too fast with increasing dimension.

Estimation

In this section we propose a two-step approach to estimate a factor model with sparse VAR idiosyncratic components and prove its consistency. For this, let $\boldsymbol{x}_t,$ $t=1,\dots,T$ be some observations and let $\boldsymbol{X}=\boldsymbol{\chi}+\boldsymbol{\Xi}$ denote the $T\times N$ matrix form of (ref). Furthermore, $\boldsymbol{\Lambda}$ denotes the $N\times r$ matrix of loadings and $\boldsymbol{F}$ denotes the $T \times r$ matrix of factors such that $\boldsymbol{\chi}=\boldsymbol{F} \boldsymbol{\Lambda}^\top$ is the matrix counterpart of (ref). Then, an estimation of the factor decomposition can be obtained by using Principal Components Analysis (PCA), see among others bai2003,bai2020simpler. The number of factors $r$ is considered as known here. {Note that the number of factors can be determined by various approaches (see Section (ref) for further discussion) and it can be estimated with probability tending to one, see among others bai2002determining.} To elaborate with the estimation, let $$\boldsymbol{X}/\sqrt{NT}=\BS U_{NT} \BS D_{NT} \BS V_{NT}^\top,$$ denote a singular value decomposition of $\boldsymbol{X}/\sqrt{NT}$ such that $\BS D_{NT}$ is a diagonal matrix with the singular values arranged in descending order on its diagonal. $\BS U_{NT}$ and $\BS V_{NT}$ are the corresponding left and right singular vectors, respectively. This can be further written as $$\BS U_{NT} \BS D_{NT} \BS V_{NT}^\top=\BS U_{NT,r} \BS D_{NT,r} \BS V_{NT,r}^\top+\BS U_{NT,N-r} \BS D_{NT,N-r} \BS V_{NT,N-r}^\top,$$ where $\BS D_{NT,r}$ is a diagonal matrix with the first $r$ largest singular values, $d_{NT,1},\dots,d_{NT,r}$, arranged in descending order on its diagonal, $\BS D_{NT,N-r}$ is a diagonal matrix with the remaining $N-r$ largest singular values, and $\BS U_{NT,r},\BS U_{NT,N-r},\BS V_{NT,r},\BS V_{NT,N-r}$ are the corresponding left and right singular vectors. Then, the estimators of a rotated version of $\boldsymbol{F}$ and $\boldsymbol{\Lambda}$ are given by $$\hat{\boldsymbol{F}}=\sqrt{T} \bm U_{NT,r} \text{ and } \hat{\boldsymbol{\Lambda}}=\sqrt{N} \BS V_{NT,r} \BS D_{NT,r},$$ such that $\hat{\boldsymbol{\chi}}=\hat{\boldsymbol{F} }\hat{\boldsymbol{\Lambda}}^\top$ and $\hat {\BS \xi}=\BS x_t-\hat{\boldsymbol{\chi}} $. This uses the normalization $\hat {\boldsymbol{F}}^\top \hat{\boldsymbol{F}}/T=\bm I_r$ and $\hat{\boldsymbol{\Lambda}}^\top \hat{\boldsymbol{\Lambda}}$ is a diagonal matrix. Consider the estimated idiosyncratic components $\hat \bm \xi$. As it is assumed that $\{\BS \xi_t\}$ follows a sparse vector autoregressive model, we estimate this sparse VAR on $\hat \bm \xi$ by regularized methods such as the (adaptive) lasso. This idea leads to the following two-step estimation procedure:

By estimating the factors $\boldsymbol{f}_t$ through standard Principal Components Analysis (PCA) and the sparse VAR models of the idiosyncratic components $\boldsymbol{\xi}_{t}$ via sparse penalized regression techniques, we combine a dense estimation approach with a sparse one. This can possibly better capture and disentangle both the dependence coming from the diverging eigenvalues of $\mathbb{E}(\boldsymbol{X}\boldsymbol{X}^{\top})$, i.e., the factors, as well as the dependence coming from the non-diverging eigenvalues of $\mathbb{E}(\boldsymbol{X}\boldsymbol{X}^{\top})$, i.e., the idiosyncratic components. The estimation of factors and loadings via PCA is a well established method in the literature, see among others stock2002forecasting,bai2003, and the common and idiosyncratic component can be estimated with rate $O_P(\max(1/\sqrt{T},1/\sqrt{N}))$. Since it is no different in this setting, we focus our presentation on the second estimation step. For a sparse stationary VAR model, deviation bounds and restricted eigenvalue conditions can be established, see among others basu2015,kock2015oracle. Given these, the consistency of the lasso can be derived and, under additional Gaussianity assumption, one obtains a rate for strict sparsity $k$ of $O_P(k \sqrt{\log(N)/T})$. However, as the idiosyncratic component $\{\BS \xi_t\}$ is not observed in our setting and hence needs to be estimated, the regression in Step 2 is performed only with the estimated idiosyncratic component. Consequently, the aforementioned results cannot be applied here. Before analyzing the second step, we in fact need to quantify the estimation error $\BS w_t:=\boldsymbol{\hat \xi_t}-\BS \xi_t\equiv \BS F \BS \Lambda^{\top}-\hat{\BS F} \hat{\BS \Lambda}^{\top}$ arising from the first step. For the consistency of the lasso, this means quantifying the estimation error $\BS w_t$ in quantities such as $\|1/T \sum_{t=1}^T (\BS \xi_t+\BS w_t)(\BS \xi_t+\BS w_t)^\top\|_{\max}$. If we simply apply the rate derived in the literature for approximate factor models, see among others stock2002forecasting,bai2003 which derive $\BS w_t=O_P(\max(1/\sqrt{T},1/\sqrt{N}))$, we would obtain $\|1/T \sum_{t=1}^T (\BS \xi_t+\BS w_t)(\BS \xi_t+\BS w_t)^\top\|_{\max}=\|1/T \sum_{t=1}^T \BS \xi_t\BS \xi_t^\top\|_{\max}+O_P(\max(1/\sqrt{T},1/\sqrt{N}))$. This may lead to a rate for the second step of $O_P(k (\sqrt{\log(N)/T}+k\max(1/\sqrt{N},1/\sqrt{T}))$. However, this can be improved if we analyze the estimation error $\BS w_t$ more closely. For this, we follow the idea of the decomposition in eq. (6) in bai2020simpler. To elaborate, we have $1/(NT) \BS X \BS X^\top \hat{\BS F}=\hat{\BS F} \BS D_{NT,r}^2$. Plugging in (ref), and using the rotation matrix

we obtain the following representation for the error between the estimated factors and a rotated version of the factors

Similarly, we obtain by symmetry for the loadings

These representations can be used to derive the order of the estimation error for the factors and loadings as it is done with a slightly different rotation matrix in bai2003. However, as our focus is on $\BS w_t:=\boldsymbol{\hat \xi_t}-\BS \xi_t$, we use these results to derive a representation of $\BS w_t$. With the obtained representation for $\BS w_t$, we can analyze more closely the estimation error of the second step. For this, note first that $\|1/T \sum_{t=1}^T (\BS \xi_t+\bm w_t)(\BS \xi_t+\bm w_t)^\top\|_{\max}\leq \|1/T \sum_{t=1}^T (\BS \xi_t)(\BS \xi_t)^\top\|_{\max}+2\|1/T \sum_{t=1}^T (\bm w_t)(\BS \xi_t+\bm w_t)^\top\|_{\max}+\|1/T \sum_{t=1}^T (\bm w_t)(\bm w_t)^\top\|_{\max}$. Bounds for the latter objects and the representation of $\BS w_t$ are given in the following Theorem (ref).

If $N=T^a$ and $a\leq \zeta-4$, we have $g(N,T,\zeta)\leq 1/\sqrt{NT}+1/T$ which means $g(N,T,\zeta)$ could be dropped in the above $O_p$ terms.

We focus here on the lasso itself but the above theorem is also helpful for obtaining rates for the de-sparsified/de-biased lasso in this framework. As mentioned previously, if we just plug-in the rate for $\BS w_t$ we would obtain the slower rate of $O_P(\max(1/\sqrt{T},1/\sqrt{N}))$. With the results above we can establish bounds for the estimation error of the second step, as done in the following Theorem (ref).

Recall that $q$ is the approximate sparsity parameter of Assumption (ref). Let us have a closer look on the bound $\|\hat \bm \beta^{(j)}- \bm \beta^{(j)}\|_1$. Consider $N=T^a, p=T^b$ for some $a,b>0$ and let the following moment condition hold $\zeta\geq 4(1+a+b)$ as well as $k \leq \sqrt{N}$. Then, for $j=1,\dots,N$ the bound simplifies to $$\|\hat \bm \beta^{(j)}- \bm \beta^{(j)}\|_1=O_P(k [\sqrt{\log(Np)/T}+k k_\xi/N]^{1-q}).$$ The first condition, i.e., $k [\sqrt{\log(Np)/T}]^{1-q}=o(1)$, is standard for approximately sparse models, see among others Corollary 2.4 in van2016estimation. The second condition, $k [k k_\xi/N]^{1-q}=o(1)$, is not standard for approximately sparse models and appears due to the estimation error of the first step. That means the estimation error of the first step is negligible if (ignoring log terms for simplicity) $k k_\xi/N \leq 1/\sqrt T$. Hence, in the restrictive case of $k_\xi=O(1)$ the estimation error of the first stage is negligible if $N\geq k \sqrt{T}$. Let us also mention that without the detailed expression for $\BS w_t$, we would obtain the estimation error of the first step being of order $k/\sqrt{N}$. As $k_\xi$ is upper bounded by $\sqrt{N}$, the derived error bound with detailed expression for $\BS w_t$ is in no case less tight than the one without and it is tighter when $k_\xi/\sqrt{N}=o(1)$.

Applications

Forecasting and Factor-augmented regression

Forecasting and factor-augmented regression is one of the most important uses of factor models, see among others stock2002forecasting. Let us consider a $h$-step ahead forecast. Then, in factor-augmented regression, we have the observables $\{\bm X_t,\bm Y_t,\bm W_t\}$ and the following regression model

Additionally, $\bm X_t$ possesses the factor structure $\bm X_t=\bm \varLambda \bm f_t+\bm \xi_t$ and $\bm W_t$ often consists of lagged values of $\bm Y_t$ and an intercept. To employ this regression model, $\bm f_t$ needs to be estimated and the implications of using estimated regressors are analyzed among others in bai2006confidence,gonccalves2014bootstrapping. $\bm Y_t$ can be here a scalar or vector and usually the regression model is estimated by least-squares which implies that $\bm W_t$ is considered as low-dimensional.

In many applications, $\bm Y_t$ itself is a subset of $\bm X_t$. Then, the model introduced in Section (ref) is an extension of the factor-augmented regression model. To elaborate, let $\bm Y_t$ be $N^Y$-dimensional and $\bm I_{N,Y} \in \mathds{R}^{N_Y\times N}$ be a selection matrix such that $\bm Y_t=\bm I_{N,Y}^\top \bm X_t$. Similarly, $\bm W_t=\bm I_{N,W}^\top \bm X_t$, $\bm I_{N,W} \in \mathds{R}^{N_W \times N}$. For illustrative purposes, we focus on the case of one lag. Then, the factor-augmented regression model (ref) reads as $$ \bm Y_{t+h}=\bm I_{N,Y}^\top \bm X_{t+h}=\bm I_{N,Y}^\top (\bm B \bm f_t+\bm A \bm I_{N,W} \bm I_{N,W}^\top \bm X_{t} + \bm \eta_{t+h}), $$ where $\bm \beta=\bm I_{N,Y}^\top \bm B$, $\bm \alpha=\bm I_{N,Y}^\top \bm A \bm I_{N,W}$, and $\boldsymbol{\varepsilon}_{t+h}=\bm I_{N,Y} \bm \eta_{t+h}$. In this set-up the factor-augmented regression can be understood as a feasible (and non-sparse) approximation of the regression of $\BS Y_{t+h}$ onto $\BS X_t$. With this motivation in mind, a natural and feasible extension of the previous model would be

where $\tilde \bm A$ is considered as sparse. The idea behind this extension is that instead of doing the selection by hand, i.e, choosing $\bm I_{N,W}$, the selection is done automatically by a data-driven selection procedure such as lasso. Since $\bm X_t$ decomposes into factor and idiosyncratic part, also $\bm \eta_t$ should decompose into such parts. Let $\bm \eta_{t+h}=\bm \varLambda \bm u_{t+h}+\bm v_{t+h}$ and $\bm u_{t+h}=\bm f_{t+h}-\bm D \bm f_t$, $\bm v_{t+h}=\bm \xi_{t+h}-\bm E \bm \xi_{t}$. Then, we have

which is equivalent to $$ 0=\bm I_{N,Y}^\top((\bm B+\tilde \bm A \bm \varLambda-\bm \varLambda \bm D)\bm f_t+(\tilde \bm A-\bm E) \bm \xi_t). $$ If $\bm f_t$ and $\bm \xi_t$ are uncorrelated, we have $\tilde \bm A=\bm E$ and $\bm B=\bm \varLambda \bm D-\tilde \bm A \bm \varLambda$. Thus, the sparse extended factor-augmented model (ref) reads as the following state-space model, which is a special case of the model described in Section (ref)

With this connection in mind, we see that a forecast of $\bm X_t$ is built upon forecasting $\bm f_t$ and $\bm \xi_t$. We are now going to present the forecast method for a one-step-ahead prediction. An $h$-step-ahead prediction can be done recursively. If the main interest is on an $h$-step-ahead forecast, a direct $h$-step-ahead forecast can be more accurate, see among others smeekes2018macroeconomic. A direct $h$-step-ahead forecast can be obtained by changing the regression equation from $t$ to $t+h$ as in (ref). Based on the estimation method proposed in the previous section, the approach is as follows. First, a standard linear one-step-ahead prediction is computed based on the estimated factors. Combining this prediction with the loadings gives a prediction of the common component. Second, the sparse VAR model is used to get a prediction of the idiosyncratic component. Finally, the sum of these two predicted components gives the prediction of the original process $\BS x_t$. To elaborate, consider first that the factors and idiosyncratic component are observed. Then, let $\BS f^{(1,p_f)}_{T+1}=\sum_{j=1}^{p_f} {\BS \Pi}_j^{(p_f)} {\BS f}_{T+1-j}$ be the linear one-step-ahead prediction based on ${\BS f}_{T},\dots,{\BS f}_{T-p_f}$, where $\sum_{j=1}^{p_f} {\BS \Pi}_j^{(p_f)} {\BS \Gamma}_{f}(i-j)={\BS \Gamma}_f(i),i=1,\dots,p_f$ and $ {\BS \Gamma}_f(i-j)=\mathbb{E} {\BS f}_{t+i} {\BS f}_{t-j}^\top$, see among others Section 11.4 in BrockwellDavis1991. Furthermore, since $\{ \BS \xi_t\}$ follows a VAR$(p)$ model, ${\BS \xi}_{T+1}^{(1)}=\sum_{j=1}^p {\BS A^{(j)}} {\BS \xi}_{T-j}$ is the one-step-ahead prediction for the idiosyncratic component. That means, $\BS X^{(1,p_f)}_{T+1}=\BS \Lambda \BS f^{(1,p_f)}_{T+1}+{\BS \xi}_{T+1}^{(1)}$ is the joint one-step-ahead prediction for $\BS X_{T+1}$ with the prediction error $\mbox{Var}(\BS X_{T+1}-\BS X_{T+1}^{(1,p_f)})=\BS \Lambda \mbox{Var}(\BS f^{(1,p_f)}_{T+1}-\BS f_{T+1}) \BS \Lambda^\top +\bm \varSigma_v$ and for a single variable $j$ we have $\mbox{Var}(\BS e_j^\top (\BS X_{T+1}-\BS X_{T+1}^{(1,p_f)}))=\BS \Lambda_j^\top \mbox{Var}(\BS f^{(1,p_f)}_{T+1}-\BS f_{T+1}) \BS \Lambda_j +\bm e_j^\top\bm \varSigma_v \bm e_j$. If $\{\BS f_t\}$ follows a VAR$(p_f)$ model, this simplifies to $\mbox{Var}(\BS X_{T+1}-\BS X_{T+1}^{(1,p_f)})=\BS \Lambda \bm \varSigma_u \BS \Lambda^\top +\bm \varSigma_v$.

Since the parameters are unknown and the factors and idiosyncratic component are latent, this approach is unfeasible but the results of Theorem (ref) and (ref) help to obtain a feasible approach. For this, we construct feasible counterparts of the prediction approach above. Let $\hat {\BS f}_{T+1}^{(1,p_f)}=\sum_{j=1}^{p_f} \hat {\BS \Pi}_j^{(p_f)} \hat {\BS f}_{T+1-j}$ be the linear one-step-ahead prediction based on $\hat {\BS f}_{T},\dots,\hat {\BS f}_{T-p_f}$, where $\sum_{j=1}^{p_f} \hat {\BS \Pi}_j^{(p_f)} \hat {\BS \Gamma}_{f}(i-j)=\hat {\BS \Gamma}_f(i),i=1,\dots,p_f$, and $\hat {\BS \Gamma}_f(i-j)=1/n \sum_{t=1+j}^{T-i} \hat {\BS f}_{t+i} \hat {\BS f}_{t-j}^\top$. Furthermore, let $\hat {\BS \xi}_{T+1}^{(1)}=\hat {\BS A} (\hat {\BS \xi}_{T}^\top,\dots,\hat {\BS \xi}_{T-p}^\top)^\top$ be the one-step-ahead prediction for the idiosyncratic component. Then, $\hat \bm X_{T+1}^{(1,p_f)}=\hat {\BS \Lambda} \hat {\BS f}_{T+1}^{(1,p_f)}+\hat {\BS \xi}_{T+1}^{(1)}$ is the joint and feasible one-step-ahead prediction for $\BS X_{T+1}$. Even though a high-dimensional time series system is considered, the interest is often in the prediction of some key time series. We quantify in the following Theorem (ref) the estimation error between the feasible and unfeasible approach for a single time series.

In relation to the error bound for $\|\hat \bm \beta^{(j)}- \bm \beta^{(j)}\|_1$ derived in Theorem (ref) only an additional $1/\sqrt{N}$ appears which arises due to the estimation of the factors.

Let us mention that $\bm X_t$ does not need to consist of series at the same time point. For instance, predicting inflation and GDP at time $t$ with other variables available up to time $t$, the vector $\bm X_t$ can be build as inflation$(t-1)$, GDP$(t-1)$ and other variables$(t)$. Then, using the proposed approach to obtain a prediction of $\bm X_{t+1}$ gives a prediction of inflation and GDP at time point $t$.

Bootstrap of factor models

{In factor-augmented regression, when inference for the regression coefficients of the factors is of interest, the estimation of the factors needs to be taken into account. In cases when $N$ is large in relation to $T$, the estimated factors can be treated as observed, see bai2006confidence for details. However, if $\sqrt{T}/N\to c, c>0$ some bias term appears which contains among others $\mathbb{E} \bm \xi_t \bm \xi_t=\BS \Gamma_\xi(0)$, see gonccalves2014bootstrapping. To assess this bias term, gonccalves2020bootstrapping propose a bootstrap algorithm that mimics the cross-sectional dependence structure of the idiosyncratic component. The bootstrap relies on an estimate of $\BS \Gamma_\xi$. For this, they assume sparsity of $\BS \Gamma_\xi(0)$ and estimate it via thresholding. If $\bm \xi_t$ is driven by a sparse VAR model, assuming sparsity on $\BS \Gamma_\xi(0)$ restricts (in a not traceable way) the sparsity of the slope parameter of the VAR model, see krampe2018bootstrap. To avoid this, two options exist and both rely also on an estimate of the variance matrix of $\bm v_t$, the innovations of the VAR process. First, $\BS \Gamma_\xi(0)$ can be estimated using the VAR structure, see estimator (6) in krampe2018bootstrap. With this estimated $\BS \Gamma_\xi(0)$, the bootstrap approach of gonccalves2020bootstrapping can be used. Second, one can extend the bootstrap approach of gonccalves2020bootstrapping and mimic not only the contemporaneous dependence structure of $\bm \xi_t$ but also the entire second-order structure of $\bm \xi_t$ by leveraging on the sparse VAR structure. For this, one can follow the bootstrap algorithm of krampe2018bootstrap (specifically their step 1 and step 2 in Section 3). Since this also mimics the dependence over time, it could improve finite sample performance. Furthermore, this extended bootstrap approach can also be used to obtain inference results for the loadings. The asymptotic normality of the loadings is derived in bai2003 and for $\BS \Lambda_i$ the asymptotic variance contains among others terms such as $1/T \sum_{s=1,t=1}^T \mathbb{E} \bm f_t \bm f_s^\top \bm \xi_{i,s} \bm \xi_{i,t}$. Under independence of factors and idiosyncratic component, this simplifies to $1/T \sum_{s=1,t=1}^T \BS \Gamma_f(t-s) \bm e_i^\top \BS \Gamma_{\xi}(t-s) \bm e_i$ and only second-order moments of the idiosyncratic component appear. Hence, a successful bootstrap approach for the loadings does not only need to mimic the dependence structure of the factors but also that of the autocovariance of the idiosyncratic component which our proposed extension achieves.}

Estimation of time series dependence networks

{For multivariate time series, networks are often used to display the connection structure between the individual time series. In these networks, each time series is represented by a node and an edge between two nodes is drawn if some form of connection between the two time series exists. Here, several approaches are available to define a connection. In the context of (Gaussian) graphical models, a connection is drawn based on the partial correlation structure, which translates in the time series context to the partial coherence structure, see brillinger1996remarks,dahlhaus2000graphical. Other approaches to defining a connection are based on Granger-causality granger1969investigating,hecq2021granger and forecast error variance decompositions, see diebold2014network. In the following, we elaborate on the connections based on partial coherences and graphical models. The partial coherence measures the strength of the linear relations between two time series after eliminating all indirect linear effects caused by all other time series of the system, taking into account all leads and lags relations.}

{ To elaborate, consider two components $u$ and $v$, then the partial coherence at frequency $\omega \in [0,2\pi]$ is given by

where $f_{u,v}^{-1}(\omega)$ denotes the $(u,v)$th element of the inverse of the spectral density matrix at frequency $\omega$. An edge is drawn between component $u$ and $v$ if $\sup_\omega R_{u,v}(\omega)\geq \delta$ for some $\delta \in [0,1)$. $\delta$ is a user-specified threshold determining which connections are important. Note that $\delta=0$ includes all non-zero connections. However, in the presence of a factor, it is most likely that $ \sup_\omega R_{u,v}(\omega)>0$ for all $u,v$ and it is of more interest to identify those which exceed some positive threshold. To inherit such a network from data, the spectral density needs to be estimated.} When the dimension of the time series is small, the spectral density matrix is often estimated by non-parametric approaches as lag-window estimators or smoothed periodograms, respectively, see among others brillinger2001time,koopmans1995spectral,hannan2009multiple,wu2018asymptotic. In a high-dimensional set-up, the problem of estimating the spectral density matrix or its inverse has been extensively investigated in the literature during the last decade. One approach is to combine the non-parametric lag-window estimators with regularization techniques developed for the covariance and precision matrix estimation, see among others sun2018large,fiecas2019spectral,zhang2020convergence. Such approaches work under the assumption that the spectral density matrix or its inverse is sparse. However, a direct sparsity assumption on the spectral density matrix or its inverse is contradicting the assumption of the existence of factors. That means the factors need to be taken into account in the estimation of the spectral density matrix. The procedure developed in the previous section can be used to obtain (under slightly modified assumptions) a consistent estimator of the inverse of the spectral density matrix. Since the VAR structure of the idiosyncratic component is used, we obtain a semiparametric estimator for the inverse of the spectral density matrix. Let us mention that the factors can be also taken into account by using a low-rank plus sparse approach applied to a smooth periodogram, see barigozzi2021algebraic. This estimator differs, however, from the one presented here in several aspects. First, the low-rank plus sparse approach describes in finite samples a different model than the approach used here, see also the discussion of low-rank plus sparse structures in Remark (ref). Second, they focus on consistency results regarding $\|\cdot\|_2$ whereas we present here also row- and column-wise consistency results i.e., consistency with respect to $\|\cdot\|_1$ and $\|\cdot\|_\infty$.

Let us begin with defining the spectral density matrix of the time series $\{\BS X_t\}$ given by (ref). The spectral density matrix of the factor process specified in Assumption (ref) is given by $$ \BS f_f(\omega)= \bigg[\sum_{j=0}^\infty \BS D^{(j)} \exp(-i j \omega)\bigg] \BS \Sigma_u \bigg[\sum_{j=0}^\infty \BS D^{(j)} \exp(i j \omega)\bigg]^\top, \quad \omega \in [0,2\pi] $$ and for the idiosyncratic component driven by a VAR$(p)$ we have $$ \BS f_\xi(\omega)=\bigg[\BS I_N-\sum_{j=1}^p \BS A^{(j)} \exp(-i j \omega)\bigg]^{-1} \BS \Sigma_v \bigg(\bigg[\BS I_N-\sum_{j=1}^p \BS A^{(j)} \exp(i j \omega)\bigg]^{-1}\bigg)^\top, $$ with the inverse $$ \BS f_\xi(\omega)^{-1}=\bigg[\BS I_N-\sum_{j=1}^p \BS A^{(j)} \exp(i j \omega)\bigg]^\top \BS \Sigma_v^{-1} \bigg[\BS I_N-\sum_{j=1}^p \BS A^{(j)} \exp(-i j \omega)\bigg]. $$ That is, the spectral density of the process $\{\BS X_t\}$ is given by

and its inverse using the Sherman–Morrison–Woodbury formula is given by

We estimate $\BS f_X(\omega)^{-1}$ by estimating $\BS f_f^{-1}$ and $\BS f_\xi^{-1}$ separately. Note that the factors lead to an unbounded $\|\BS f_X(\omega)\|_2$ for growing dimension but the inverse is stable, i.e., $\|\BS f_X(\omega)^{-1}\|_2$ is bounded.

As it is of fixed dimension $r$, the spectral density $\BS f_f$ or its inverse can be estimated by classical methods such as non-parametric lag-window estimators. For this, let $K$ be a kernel fulfilling Assumption 1 in wu2018asymptotic. That is, $K$ is an even and bounded function with bounded support in $(-1,1)$, continuous in $(-1,1)$, $K(0)=1, \kappa=\int_{-1}^1 K^2(u) du<1$, and $\sum_{l \in \mathds{Z}} \sup_{|s-l|<1} |K(l\omega)-K(s\omega)|=O(1)$ as $\omega \to 0$. Furthermore, let $B_T=T^b, b \in (0,1)$ be the lag-window size fulfilling Assumption 2 in wu2018asymptotic. Then, a spectral density estimator is given by

where $\hat{ \BS \Gamma}_f(h)$ is the sample autocovariance function $\hat { \BS\Gamma}_f(h)=1/T\sum_{t} \hat{\BS f}_{t+h} \hat{ \BS f}_t^T$. Based on observations $\BS f_1,\dots,\BS f_T$, let $\tilde {\BS f}_f(\omega)=\frac{1}{2\pi}\sum_{h=-T+1}^{T-1} K\left(\frac{u}{B_T}\right) \exp(-ih \omega) \tilde{ \BS \Gamma}_f(h),\; \tilde{\BS \Gamma}_{f} (h)=1/T\sum_{t} \BS f_{t+h} \BS f_t^T,$ be the (unfeasible) estimator of $\BS f_f$. Then, the results of wu2018asymptotic give that $\|\tilde {\BS f}_f(\omega)-\BS f_f(\omega)\|_{\max}=O_P(\sqrt{B_T/T})$. With this result and noting that $\{\BS f_t\}$ is a process of fixed dimension $r$, consistency of $\hat {\BS f}_f(\omega)$ follows by Lemma (ref), see Lemma (ref) in the appendix for details.

As mentioned, we use the VAR structure of the idiosyncratic component to estimate its spectral density matrix. In the previous section, we showed that the VAR parameters of the idiosyncratic component can be estimated row-wise consistently, i.e., consistency of $\boldsymbol{\frak{A}}$ for the matrix norm $\|\cdot \|_\infty$. However, the estimation of the spectral density requires additional column-wise consistency, that is consistency of $\BS A^{(j)},j=1,\dots,p,$ with respect to $\|\cdot\|_1$. Such a column-wise consistency requires additional sparsity assumptions, see also krampe2020statistical for a discussion. Furthermore, a parametric estimation of the spectral density matrix of a VAR process requires an estimate of the covariance or precision matrix of the residual process $\{\BS v_t\}$. Since our focus is on the estimation of the inverse of the spectral density matrix, we estimate the precision matrix and formulate sparsity assumption on this matrix. See Assumption (ref) for the exact definition of the additional sparsity assumptions.

As mentioned, the precision matrix of the residuals $\{\BS v_t\}$ needs to be estimated. The residuals can be estimated by $\hat {\BS v}_t=\hat {\BS \xi}_t-\sum_{j=1}^p \hat {\BS A}^{(j)} \hat {\BS \xi}_{t-j}, t=p+1,\dots,T$. Then, based on these estimated residuals, procedures like graphical lasso of friedman2008sparse or (A)CLIME of cai2011constrained,cai2016estimating,cai2016estimating2 can be used. In the proofs we consider the CLIME method and denote this estimator by $\hat {\BS \Sigma}_v^{-1,CLIME}$ but similar results can be established for the graphical lasso estimator. Then, we construct the following estimator for $\BS f_\xi^{-1}(\omega)$

where $\hat{\BS A}^{(thr,j)}=(\operatorname{THR}_{\lambda_\xi}\hat{\BS A}^{(j)})$ and $\operatorname{THR}_{\lambda_\xi}$ is a thresholding function with threshold parameter $\lambda_\xi$ fulfilling the conditions $(i)$ to $(iii)$ in Section 2 in cai2011adaptive. For instance, such a thresholding function can be the adaptive lasso thresholding function given by $\operatorname{THR}_{\lambda_\xi}^{al}(z)=z(1-|{\lambda}/z|^\nu)_+$ with $\nu\geq1$. Soft thresholding ($\nu=1$) and hard thresholding ($\nu=\infty$) are boundary cases of this function. This thresholding functions act by thresholding every element of the matrix $\hat{\BS A}^{(j)}$ and it results in a row- and column-wise consistent estimation of the VAR slope matrices. In Lemma (ref) in the appendix, we present the error bounds $\|\hat{\BS f}_\xi^{-1}(\omega)-{\BS f}_\xi^{-1}(\omega)\|_\infty$ and $\|\hat{\BS f}_\xi^{-1}(\omega)-{\BS f}_\xi^{-1}(\omega)\|_2$. Finally, replacing in (ref) all quantities with the estimators discussed above leads to our final estimator of the inverse of the spectral density matrix of $\{\BS X_t\}$. Its error bounds are given in the following Theorem (ref). We only present here explicitly the rate for a simplified case. In the general case, an explicit rate can be obtained by inserting the results of Lemma (ref) and Theorem (ref). Since it leads to a lengthy and not insightful expression, we omit it here. The rate is dominated by the estimation error of the sparse VAR and it is similar to the one in Theorem (ref). However, the rate is more affected by the sparsity parameter in the sense that its maximum growth rate is less for the spectral density than it is for prediction. Maximum growth rate refers here to the maximal rate of sparsity for which consistency can be achieved.

{Joint Selection: Number of Factors & Lag-length}

{The present context clearly requires the selection of the number of factors within the PCA step as well as the order of the VAR for the idiosyncratic component. In the literature, there is an abundance of methods available for both. Among others, the seminal work of bai2002determining introduced information criteria for a data-driven specification of the number of factors and it is perhaps the most employed method in practice. Further refinements of this method can be found in hallin2007determining,alessi2010improved. Information criteria can also be used to specify the order of a VAR. For instance, hecq2021granger propose to marginalize the (high-dimensional) VAR into a sequence of AR(p) regressions and select the lag-length via an approximated Bayesian information criterion (BIC). The consistency of the BIC has been proved in wang2009shrinkage. Under a few technical conditions on the divergence speed of the model dimension and the size of non-zero coefficients, they show how a slightly modified BIC can identify the true model consistently even when the dimension diverges. } We propose here a unified procedure able at the same time to consistently estimate the lag-length as well as the number of factors. For a given lag-length and number of factors, the penalty parameter can {also} be chosen with {an information criterion as AIC or BIC but this necessarily needs to be distinct from the joint information criteria for the number of factors and lag-length hence we briefly discuss it first. Let $\bm \xi_{t,S}^v$ be the subvector containing those columns of $\bm \xi_t^v$ belonging to the set $S$. Let further $\hat{S}$ be the active set identified by the lasso for a given $\lambda$. Then the value $\lambda^{IC}$ chosen by information criteria is found as

where $df$ represents the degrees of freedom after the penalization, i.e., the cardinality of the estimated active set. $C_T$ is the penalty specific to each criterion, where the most popular choices are: $C_T=2$, the Akaike information criterion (AIC) by akaike1974new; $C_T=\log(T)$, the Bayesian information criterion (BIC) by schwarz1978estimating.\footnote{Note: for non-Gaussian distributions, the residual sum is often used as a proxy for the likelihood.} The slight modification of the BIC proposed in wang2009shrinkage also holds for penalized estimators as the lasso, thus making it consistent asymptotically in both $N$ and $T$. } With regard to the number of factors and lag-length, as in some applications, the focus is more on forecasting a small subset of time series of the system, we present here two approaches: a global approach which gives a single lag-length and number of factors for the entire system and a local approach in which the lag-length or number of factors may differ across the time series. We present the two approaches first and then discuss their differences.

We consider that the factors are driven by a VAR model that is $\BS f_t=\sum_{j=1}^{p_f} \BS \Pi_j \BS f_{t-j}+\BS v_{t-j}$. That means we have two lag-lengths to choose: $p$ and $p_f$. The one-step ahead forecast error of model (ref) for the $i$th component is given by $$ \mbox{Var}\Big(x_{i,t}-\sum_{j=1}^{p_f} \BS \Lambda_i^\top \BS \Pi_j \BS f_{t-j}-\sum_{j=1}^p \BS e_i^\top \BS A^{(j)} \BS \xi_{t-j}\Big).$$ If we treat the factors and idiosyncratic components as known, we have to estimate for all components the parameters $\BS \Lambda \in \mathds{R}^{N\times r}, \BS \Pi_1,\dots,\BS \Pi_{p_f} \in \mathds{R}^{r\times r}, \BS A^{(1)},\dots,\BS A^{(p)} \in \mathds{R}^{N\times N}$. Note that $\BS A_1,\dots,\BS A_p$ are sparse. That means in total we have $(N+r p_f)r+ \sum_{j=1}^p \|\BS A_j\|_0$ parameters for all components. For a single component, we treat $\BS \Lambda_i \BS \Pi_j$ as $r$-dimensional vectors which gives in total for the $j$th component $r p_f+ \sum_{j=1}^p \|\bm e_i^\top\BS A_j\|_0$ parameters. {The sparsity of the idiosyncratic component has the important implication that $\sum_{j=1}^p \|\bm e_i^\top\BS A_j\|_0$ grow much slower than $Np$. To be precise, the error bounds in Theorem (ref) imply that only rates slower than $\sqrt{T}$ are reasonable. Hence, the number of parameters considered grow slower than the sample size and consequently, this fits into the framework of wang2009shrinkage and their modified BIC. Note however, that the results of wang2009shrinkage are derived under an i.i.d. set-up and also the pre-selection of the penalty parameter $\lambda_n$ is not taken into account here.} In this modified BIC set-up $C_T$ denotes a slowly diverging series which is discussed shortly. This motivates the following global information criteria

For the $i$th component we obtain the following local information criteria

In practice, the minimum is evaluated over a finite grid. That means one sets a maximal number of factors $r_{\max}$ and maximal lag-lengths $p_{\max},p_{f,\max}$. If one sets $r_{\max}=0$ or $p_{\max}=0$, this criteria can also be used to fit plain sparse VAR models or plain factor models, respectively. The series $C_T$ can be diverging very slowly and wang2009shrinkage suggest for instance, $\log(\log(T))$. We would like to consider the diverging dimension as well and follow a similar route as bai2002determining. So we set $C_T=c\frac{\log(NT/(N+T))}{\log(T)}$ with $c=1/2$. Note that for the global approach the factors are penalized by $(p_f+N)\log(NT/(N+T))/(NT)$. This also implies that this series fits into the penalization function framework of Theorem 2 in bai2002determining required to obtain a consistent estimation of the number of factors, i.e., this series converges to $0$ for $N,T\to \infty$ and diverges if scaled by $\min(N,T)$.

Some remarks on these two information criteria. First, the local approach requires for the $i$th component only an estimation of $\bm e_i^\top \hat{\bm A}^{(j)}$. If the interest is only in some time series of the system, this reduces the computational burden. Second, if the number of factors differs among the time series, the entire system cannot be written as a factor model with a maximal number of factors and a maximal number of lags. Third, the local approach takes into account that large data sets come as a -- in some sense arbitrary -- collection of series and it is most likely that some series are not driven by factors or a small lag-length is sufficient. However, the additional cross-section average in the global approach also leads to more stable results. In simulations, the local approach outperforms the global approach, see Section (ref) for further discussion.

A Combined Approach for Single Time Series

In light of the previous two paragraphs, we seek for a unified procedure able at the same time to consistently estimate the lag-length as well as the number of factors. For a given lag-length and number of factors, the penalty parameter can be chosen with the approaches discussed in the next subsection and we take this for granted here. {This means the VAR models for the different lags are considered as sparse here.} Furthermore, in some applications the focus is more on forecasting a small subset of time series of the system. To take this into account we present here two approaches. A global approach which gives a single lag-length and number of factors for the entire system and a local approach in which the lag-length or number of factors may differ across the time series. We present the two approaches first and then discuss their differences.

We consider that the factors are driven by a VAR model that is $\BS f_t=\sum_{j=1}^{p_f} \BS \Pi_j \BS f_{t-j}+\BS v_{t-j}$. That means we have two lag-lengths to choose: $p$ and $p_f$. The one-step ahead forecast error of model (ref) for the $i$th component is given by $$ \mbox{Var}\Big(x_{i,t}-\sum_{j=1}^{p_f} \BS \Lambda_i^\top \BS \Pi_j \BS f_{t-j}-\sum_{j=1}^p \BS e_i^\top \BS A^{(j)} \BS \xi_{t-j}\Big).$$ If we treat the factors and idiosyncratic components as known, we have to estimate for all components the parameters $\BS \Lambda \in \mathds{R}^{N\times r}, \BS \Pi_1,\dots,\BS \Pi_{p_f} \in \mathds{R}^{r\times r}, \BS A^{(1)},\dots,\BS A^{(p)} \in \mathds{R}^{N\times N}$. Note that $\BS A_1,\dots,\BS A_p$ are sparse. That means in total we have $(N+r p_f)r+ \sum_{j=1}^p \|\BS A_j\|_0$ parameters for all components. For a single component, we treat $\BS \Lambda_i \BS \Pi_j$ as $r$-dimensional vectors which gives in total for the $j$th component $r p_f+ \sum_{j=1}^p \|\bm e_i^\top\BS A_j\|_0$ parameters. {The sparsity of the idiosyncratic component has the important implication that $\sum_{j=1}^p \|\bm e_i^\top\BS A_j\|_0$ grow much slower than $Np$. To be precise, the error bounds in Theorem (ref) imply that only rates slower than $\sqrt{T}$ are reasonable. Hence, the number of parameters considered grow slower than the sample size and consequently, this fits into the framework of wang2009shrinkage and their modified BIC. Note however, that the results of wang2009shrinkage are derived under an i.i.d. set-up and also the pre-selection of the penalty parameter $\lambda_n$ is not taken into account here.} In this modified BIC set-up $C_T$ denotes a slowly diverging series which is discussed shortly. This motivates the following global information criteria

For the $i$th component we obtain the following local information criteria

In practice, the minimum is evaluated over a finite grid. That means one sets a maximal number of factors $r_{\max}$ and maximal lag-lengths $p_{\max},p_{f,\max}$. If one sets $r_{\max}=0$ or $p_{\max}=0$, this criteria can also be used to fit plain sparse VAR models or plain factor models, respectively. The series $C_T$ can be diverging very slowly and wang2009shrinkage suggest for instance, $\log(\log(T))$. We would like to consider the diverging dimension as well and follow a similar route as bai2002determining. So we set $C_T=c\frac{\log(NT/(N+T))}{\log(T)}$ with $c=1/2$. Note that for the global approach the factors are penalized by $(p_f+N)\log(NT/(N+T))/(NT)$. This also implies that this series fits into the penalization function framework of Theorem 2 in bai2002determining required to obtain a consistent estimation of the number of factors, i.e., this series converges to $0$ for $N,T\to \infty$ and diverges if scaled by $\min(N,T)$.

Some remarks to these two information criteria. First, the local approach requires for the $i$th component only an estimation of $\bm e_i^\top \hat{\bm A}^{(j)}$. If the interest is only in some time series of the system, this reduces the computational burden. Second, if the number of factors differs among the time series, the entire system cannot be written as a factor model with maximal number of factors and maximal number of lags. Third, the local approach takes into account that large data sets come as a -- in some sense arbitrary -- collection of series and it is most likely that some series are not driven by factors or a small lag-length is sufficient. However, the additional cross-section average in the global approach also leads to more stable results. In simulations as well as for the real data set considered, the local approach outperforms the global approach, see Section (ref) for further discussion.

Penalty Tuning

Another important aspect to consider within the framework of $\ell_1$-regularizations is the choice of the tuning parameter $\lambda$. The latter should be set in order to balance between the fit of the model and its complexity, thus trading off bias with variance. Whenever the tuning parameter is large in its magnitude, the consequence is strong variable selection, i.e., many potentially relevant variables might be set to zero by the regularization technique (e.g lasso), thus implying a larger estimation bias. In parallel, when instead $\lambda\approx 0$, no variable selection is performed and thus regularization techniques such as lasso converge in the limit to the standard OLS estimator. Among the most popular techniques to tune $\lambda$ is cross-validation (CV). While CV have seen a surge of applications in statistics in the last decades, it can suffer from some shortcomings. It is often computationally demanding, especially in high dimensions, given it has to recursively train and validate on batches of the sample. Also, it needs to be adapted for different data, e.g., $K$-fold cross validation needs to be used for time series.

Alternatively, a fast and reliable way of tuning $\lambda$ is by minimizing an information criterion (IC). Let $\bm \xi_{t,S}^v$ be the subvector containing those columns of $\bm \xi_t^v$ belonging to the set $S$. Let further $\hat{S}$ be the active set identified by the lasso for a given $\lambda$. Then the value $\lambda^{IC}$ chosen by information criteria is found as

where $df$ represents the degrees of freedom after the penalization, i.e., the cardinality of the estimated active set. $C_T$ is the penalty specific to each criterion, where the most popular choices are: $C_T=2$, the Akaike information criterion (AIC) by akaike1974new; $C_T=\log(T)$, the Bayesian information criterion (BIC) by schwarz1978estimating.\footnote{Note: for non-Gaussian distributions the residual sum is often used as a proxy for the likelihood.} As for the case of the lag-length selection in Remark (ref), the slight modification of the BIC proposed in wang2009shrinkage also holds for penalized estimators as the lasso, thus making it consistent asymptotically in both $N$ and $T$.

\fi

Numerical Results

All results presented in this section are based on implementations in R R. We compute the data generating processes (DGPs) at random and consider the following model class: $\BS x_t=\BS \Lambda \BS f_t+\BS \xi_t, \BS f_t=\sum_{j=1}^{p_f} \BS \Pi^{(j)} \BS f_{t-j}+\BS u_t, \BS \xi_t=\sum_{j=1}^p \BS A^{(j)} \BS \xi_{t-j}+\BS v_t$. The innovations $\{\BS u_t\}, \{\BS v_t\}$ are generated as Gaussian processes and $\BS \Sigma_u=\mbox{Var}(\BS u_t)$ is generated as a positive definite matrix with eigenvalues in the range $1$ to $10$, using the implementation of the package clusterGeneration Clustergen. If not denoted otherwise, sparsity of a matrix is obtained by setting entries -- beginning with the absolute smallest values -- to zero such that the specified amount of sparsity is obtained. The entries of $\bm A^{(j)}$ are generated randomly using a $t$ distribution with 3 degrees of freedom. After sparsifying, the matrices are rescaled to fit the eigenvalue conditions of $0.8$. In real data, it is often observed that for a component of a multivariate time series the history of the component itself is quite an important predictor. That means that the diagonals of $\bm A^{(j)},j=1,\dots,p$ are (at least for one $j$) often non-sparse. To take this into account we put more weight onto the diagonal of $\bm A^{(1)}$ by adding $0.4 \BS I_N$ before sparsifying the randomly generated matrix. This results in a much more dominant diagonal and the diagonal of $\bm A^{(1)}$ is more dominant the smaller $p$ is.

Furthermore, we consider the following specifications:

Note that if the sparsity parameter is of similar size as the dimension, we have no sparsity. Also $r=0$ gives a pure sparse and $p=0$ a pure factor case. Dropping unnecessary combinations, e.g., varying the sparsity for $p=0$, we end up in $2352$ different set-ups for the DGP. We run each set-up $100$ times which results in $235200$ different DGPs. To evaluate the performance, we consider the average one-step ahead prediction error of the first ten time series. To compute the one-step ahead prediction error a test set of $10000$ time points is used. That is, the one-step ahead prediction error of component $i$ is given by $MSFE_{\boldsymbol{{x_i}}}= \Big[ \frac{1}{10000} \sum_{t=1}^{10000} (\hat{{{x}}}_{i,T+t}^{(1)}-{x}_{i,T+t}^{(1)})^2\Big].$ These are then averaged over the $10$ components as well as over the $100$ DGPs of each set-up.

We consider the following models to predict:

{$\boldsymbol{F_{BN}AR}$} use the information criteria of bai2002determining to determine the number of factors. Let us mention that in preliminary simulations we also considered the method of alessi2010improved as well as the global information criteria of Section (ref) to determine the number of factors. The obtained results are almost identical to the ones with the information criteria of bai2002determining. So for this specification we focus the presentation on the criteria of bai2002determining only. Furthermore, for {$\boldsymbol{L_{sel}}$}, {$\boldsymbol{FL_{sel}}$} we also considered the global information criteria of Section (ref) but do not present its result here. The findings here can be summarized as follows. The local information criteria outperforms the global criteria and the differences are larger for small sample sizes and dimensions. In the following, we present the MSFE-results in relation to the MSFE of {$\boldsymbol{FL_{sel}}$}. That means values larger than $1$ indicate a performance worse than {$\boldsymbol{FL_{sel}}$} and values smaller than $1$ vice versa. The overall performance is summarized in Table (ref).

The relative performance overall $2352$ different DGP set-ups is displayed in Figure (ref)-(ref). Each dot represents the relative MSFE for one DGP set-up averaged over the runs. The set-ups are sorted by sample size$(T)$, lag-length of the idiosyncratic part$(p)$, sparsity$(k)$, lag-length of the factors$(p_f)$, number of factors$(r)$, and dimension $(N)$. The obtained groups for $T,p,$ and $k$ are highlighted by vertical bars and the specific parameter values are given at the bottom of the figure. This sorting is chosen because these specification parameters matter the most in the sense that the results can differ substantially among different specification of the parameter values.

Let us discuss the three Figures (ref) to (ref) starting with the performance relation of {$\boldsymbol{AR_{BIC}}$} and {$\boldsymbol{FL_{sel}}$}. If the sample size is small ($T=100$), {$\boldsymbol{FL_{sel}}$} outperform {$\boldsymbol{AR_{BIC}}$} only in the case of $p=3$. In all other cases it performs equally good or even worse. This behavior changes for the larger sample size settings. Here, {$\boldsymbol{FL_{sel}}$} performs equally well in the cases of no sparsity or no dependence and clearly outperforms {$\boldsymbol{AR_{BIC}}$} in the other cases with a smaller MSFE of 40% or more. \\ In cases in which factors are present, {\emph{$\boldsymbol{FL_{sel}}$}} is outperforming {\emph{$\boldsymbol{L_{sel}}$}}. The outperformance do not differ much for different sample sizes or lag length and the MSFE is around 10% smaller for {\emph{$\boldsymbol{FL_{sel}}$}}.\\ If the sample size is small and the idiosyncratic component is mainly driven by a diagonal VAR (note the construction of slope matrices) {\emph{$\boldsymbol{F_{BN}AR}$}} outperforms {\emph{$\boldsymbol{FL_{sel}}$}}. It is the other way around for all other cases, i.e., for a less diagonal dominant VAR model and also for larger sample sizes. Then, except for the case where the idiosyncratic component has no dependency or no sparsity, {\emph{$\boldsymbol{FL_{sel}}$}} strongly outperforms {\emph{$\boldsymbol{F_{BN}^{ARfilt}}$}} with a 10% to 40% smaller MSFE.

To conclude, for the smaller sample size ($T=100$) the additional modelling of the idiosyncratic parts of {$\boldsymbol{FL_{sel}}$} does not always pay off but for the larger sample size ($T=200$) {$\boldsymbol{FL_{sel}}$} does perform best among all competitors and if there is dependence in the idiosyncratic component the gain can be quite substantial.

Conclusion

We blend the dense dimensionality reduction of factor models with the one of sparsity-inducing high-dimensional VARs. We propose a factor model whose factors and relative loadings are estimated via standard principal components while its idiosyncratic components are assumed to follow a high-dimensional sparse VAR model and are thus estimated via $\ell_1-$norm regularization techniques such as the adaptive lasso. We derive error bounds of this estimation procedure and show in which situations the lasso suffers from the estimation of the idiosyncratic components. {We discuss the implications of our model to forecasting, factor augmented regression, bootstrapping factor models and semi-parametric estimation of the inverse of the spectral density matrix.} To choose the number of factors and the lag-length of the VAR, we propose a unified procedure able to simultaneously estimate both. In simulations, we compare the performance of our proposed method with several workhorse forecasting models in the literature and find that the advantage of the procedure proposed can be substantial for moderate to large sample sizes.

\bigbreak

\bigbreak {\bf Acknowledgments.} We thank the participants of the workshop “Dimensionality Reduction and Inference in High-Dimensional Time Series" at Maastricht University for very helpful comments. The research of the first author was supported by the Research Center (SFB) 884 “Political Economy of Reforms”(Project B6), funded by the German Research Foundation (DFG). Furthermore, the first author acknowledges support by the state of Baden-W{\"u}rttemberg through bwHPC.