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.
79,094 characters · 17 sections · 3 citation commands
Sparse Approximate Factor Estimation for High-Dimensional Covariance Matrices
\def\spacingset#1{ {#1}} \spacingset{1}
\footnotetext[1]{Financial support by the Graduate School of Decision Sciences (GSDS), the German Science Foundation (DFG) and the German Academic Exchange Service (DAAD) is gratefully acknowledged. For helpful comments on an earlier draft of the paper we would like to thank Lyudmila Grigoryeva and Karim Abadir. The usual disclaimer applies.}
\footnotetext[1]{Department of Economics, Universit\"atsstra\ss e 1, D-78457 Konstanz, Germany. Phone: +49-7531-88-2657, email: [email removed].} \footnotetext[2]{Department of Economics, Universit\"atsstra\ss e 1, D-78457 Konstanz, Germany.}
{\it Keywords:} Approximate Factor model, weak factors, $l_{1}$-regularization, high dimensional covariance matrix, portfolio allocation \\ {\it JEL classification: } C38, C55, G11, G17
\spacingset{1.37} {5pt plus 0.3ex}
The estimation of high-dimensional covariance matrices and their inverses (precision matrices) has recently received a great attention. In economics and finance, it is central for portfolio allocation, risk measurement, asset pricing and graphical network analysis. The list of important applications from other areas of research includes, for example, the analysis of climate data, gene classification and image classification. What appears to be a trivial estimation problem for a large sample size $T$ and a low dimensional vector of covariates, turns out to be demanding, if $N$ is of the same order of magnitude or even larger than $T$. In these cases, the sample covariance matrix becomes nearly singular and estimates the population covariance matrix poorly. Moreover, assumptions of standard asymptotic theory with $T \to \infty $, holding $N $ fixed, turns out to be inappropriate and have to be replaced by assumptions allowing for both, $T$ and $N$, approaching infinity.
In recent years numerous studies proposed alternative estimation approaches for high-dimensional covariance matrices, which differ in the way of bounding the dimensionality problem. Two major approaches are factor models imposing a lower dimensional factor structure for the underlying multivariate process and regularization strategies for the parameters of the covariance matrix or its eigenvalues (see \citemain{Fan/Liao/Liu2016} for a recent survey on the estimation of large covariances and precision matrices). In this paper, we present an effective novel approach to the estimation of high-dimensional covariances, which profits from both branches of the literature. Our sparse approximate factor (SAF) approach to the estimation of high-dimensional covariance matrices is based on $l_1$-regularization of the factor loadings and thereby is able to account for weak factors and shrinks elements in the covariance matrix towards zero.
Approaches to obtain consistent estimators by imposing a sparse structure on the covariance matrix directly include\nocitemain{Bickel2008}\nocitemain{BickelLevina2008} Bickel2008 (Bickel2008, BickelLevina2008), \citemain{Cai/Liu2011} and \citemain{Cai/Zhou2012}. These thresholding approaches are shrinking small elements in the covariance matrix exactly to zero. While this may be a reasonable strategy, e.g. for genetic data, this assumption may not be appropriate for economic or financial data, where variables are driven by common underlying factors. Such a feature may be more appropriately captured by covariance matrices based on factor representations.
In the literature on factor based covariance estimation \citemain{Fan2008a} consider the case of a strict factor representation with observed factors. This approach requires knowledge of additional observable variables (e.g. the Fama-French factors in the asset pricing framework), which may be an additional source of misspecification. Moreover, strict factor model representations impose the overly strong assumption of strictly uncorrelated idiosyncratic errors. This assumption was relaxed in \citemain{Fan2011a} and \citemain{FanLiaoMincheva2013}, who propose a covariance estimator based on an approximate factor model representation. While \citemain{Fan2011a} shrink the entries of the covariance matrix of the idiosyncratic errors to zero using the adaptive thresholding technique by \citemain{Cai/Liu2011}, the approach proposed in \citemain{FanLiaoMincheva2013} rests on the more general principal orthogonal complement thresholding method (POET) to allow for sparsity in the covariance matrix of the idiosyncratic errors.
Our SAF covariance matrix estimator extends the existing framework on factor based approaches by imposing sparsity on both, the factor loadings and the covariance matrix of the idiosyncratic errors. Unlike imposing sparsity for the covariance matrix directly by thresholding or $l_1$-norm regularization, the $l_1$-regularization of the factor loadings does not necessarily imply zero entities of the covariance matrix, but simply reduces the dimensionality problem in the estimation of the factor driven part of the covariance matrix. Moreover, the sparsity in the matrix of factor loadings allows for weak factors, which only affect a subset of the observed variables. Thus the SAF-approach relaxes the identifying assumption on the pervasiveness of the factors in the standard framework. This further implies that the eigenvalues of the covariance matrix corresponding to the common component are allowed to diverge at a slower rate than commonly considered (i.e. slower than $\mathcal{O}(N)$).
The recent paper by \citemain{FanLiuWang2018} claims that the relaxation of the pervasiveness assumption in the approximate factor model framework is the next major concern which should be addressed in future research. Hence, in this paper we focus exactly on this issue and build a bridge between the standard factor model and a relaxed pervasiveness assumption.
The weaker conditions on the eigenvalues allow us to derive the consistency for the SAF covariance matrix estimator under the average Frobenius norm under rather mild regularity conditions. To our knowledge this convergence result is new. Because of the fast diverging eigenvalues for estimators based on the approximate factor model, convergence has only be shown under the weaker weighted quadratic norm but not for the more general Frobenius norm (see, e.g. \citemain{FanLiaoMincheva2013}). As a byproduct of our proof for the SAF covariance matrix estimator, we also prove the consistency for the estimators of the sparse factor loadings, the factors and the covariance matrix of the idiosyncratic errors.
The favorable asymptotic properties of the SAF covariance matrix estimator are well supported by our Monte Carlo study based on different dimensions and alternative designs of the population covariance matrix. More precisely, the SAF covariance matrix estimator yields the lowest difference in the Frobenius norm to the true underlying covariance matrix compared to several competing estimation strategies.
Finally, in an empirical study on the portfolio allocation problem, we show that the SAF covariance matrix estimator is a superior choice to construct the weights of the Global Minimum Variance Portfolio (GMVP) for low and large dimensional portfolios. Based on returns data from the S&P 500 the estimator uniformly outperforms portfolio strategies based on alternative covariance estimation approaches and modeling strategies including the $1/N$-strategy in terms of different popular out-of-sample portfolio performance measures.
The rest of the paper is organized as follows. In Section (ref) we introduce the approximate factor model approach and show how sparsity can be obtained with respect to the factor loadings matrix by $l_{1}$-regularization. Section (ref) discusses the theoretical setup and provides the convergence results. In Section (ref), we present Monte-Carlo evidence on the finite sample properties of our new covariance estimator, while in Section (ref) we show the performance of our approach when applied to the empirical portfolio allocation problem. Section (ref) summarizes the main findings and gives an outlook on future research.
Throughout the paper we will use the following notation: $\pi_{\max}(\operatorname{\pmb{A}})$ and $\pi_{\min}(\operatorname{\pmb{A}})$ are the maximum and minimum eigenvalue of a matrix $\operatorname{\pmb{A}}$. Further, $\left\lVert \operatorname{\pmb{A}}\right\lVert$, $\left\lVert \operatorname{\pmb{A}}\right\lVert_F$ and $\left\lVert \operatorname{\pmb{A}}\right\lVert_1$ denote the spectral, Frobenius and the $l_1$-norm of $\operatorname{\pmb{A}}$, respectively. They are defined as $\left\lVert \operatorname{\pmb{A}}\right\lVert = \sqrt{\pi_{\max}(\operatorname{\pmb{A}}'\operatorname{\pmb{A}})}$, $\left\lVert \operatorname{\pmb{A}}\right\lVert_F = \sqrt{\textup{\text{tr}}\left(\operatorname{\pmb{A}}'\operatorname{\pmb{A}}\right)}$ and $ \left\lVert \operatorname{\pmb{A}}\right\lVert_1= \max_{j} \sum_i |a_{ij} |$. For some constant $c > 0$ and a non-random sequence $b_N$, we use the notation $b_N = \mathcal{O}(N)$, if $N^{-1} b_N \to c$, for $N \to \infty$. Moreover, $b_N = o(N)$, if $N^{-1} b_N \to 0$, for $N \to \infty$. Similarly, for a random sequence $d_N$, we say $d_N = \mathcal{O}_p(N)$, if $N^{-1} d_N \overset{p}{\to} c$, for $N \to \infty$ and $d_N = o_p(N)$, if $N^{-1} d_N \overset{p}{\to} 0$, for $N \to \infty$, where $\overset{p}{\to}$ denotes convergence in probability.
The following analysis is based on the approximate factor model (AFM) proposed by \citemain{Chamberlain1983} to obtain a lower dimensional representation of a possibly high-dimensional covariance matrix. Let $x_{it}$ be the $i$-th observable variable at time $t$ for $i = 1, \dots, N$ and $t = 1, \dots, T$, such that $N$ and $T$ denote the sample size in the cross-section and in the time dimension, respectively. The AFM is given by:
where $\operatorname{\pmb{\lambda}}_i$ is a $(r \times 1)$-dimensional vector of factor loadings for variable $i$ and $\operatorname{\pmb{f}}_t$ is a $(r \times 1)$-dimensional vector of latent factors at time $t$, where $r$ denotes the number of factors common to all variables in the model. Typically, we assume that $r$ is much smaller than the number of variables $N$. Finally, the idiosyncratic component $u_{it}$ accounts for variable-specific shocks, which are not captured by the common component $\operatorname{\pmb{\lambda}}_i' \operatorname{\pmb{f}}_t$. The AFM allows for weak serial and cross-sectional correlations among the idiosyncratic components with a dense covariance matrix of the idiosyncratic error term vector, $\operatorname{\pmb{\Sigma}}_u = \operatorname*{Cov}{\left[(u_{1t}, u_{2t}, \ldots u_{Nt})'\right]}$. In matrix notation, (ref) can be written as:
where $\operatorname{\pmb{X}}$ denotes a $(N \times T)$ matrix containing $T$ observations for $N$ weakly stationary time series. It is assumed that the time series are demeaned and standardized. $\operatorname{\pmb{F}} = (\operatorname{\pmb{f}}_1, \dots, \operatorname{\pmb{f}}_T)'$ is referred to as a $(T \times r)$-dimensional matrix of unobserved factors, $\operatorname{\pmb{\Lambda}} = (\operatorname{\pmb{\lambda}}_1, \dots, \operatorname{\pmb{\lambda}}_N)'$ is a $N \times r$ matrix of corresponding factor loadings and $\operatorname{\pmb{u}}$ is a $(N \times T)$-dimensional matrix of idiosyncratic shocks.
There are several estimation approaches for a factor model as given by (ref). The principal component analysis (PCA)\footnote{See e.g., \citemain{Bai2002} for a detailed treatment of the PCA in approximate factor models.} and the quasi-maximum likelihood estimation (QMLE) under normality (see i.e. \citemain{Bai2016a}) are the two most popular ones. In the following, we pursue estimating the factor model by QMLE. This allows us to introduce sparsity in the factor loadings by penalizing the likelihood function. Moreover, contrary to PCA, all model parameters including the covariance matrix $\operatorname{\pmb{\Sigma}}_u$ can be estimated jointly, while PCA-based second stage estimates of $\operatorname{\pmb{\Sigma}}_{u}$ require consistent estimation of $\operatorname{\pmb{\Lambda}}$ and $\operatorname{\pmb{F}}$ in the first stage. This, however, may be problematic for the case of a relatively small $N$, because $\operatorname{\pmb{F}}$ can no longer be estimated consistently (\citemain{Bai2016}).
The negative quasi log-likelihood function for the data in the AFM is defined as:
where $\operatorname{\pmb{S}}_{x} = \frac{1}{T} \sum_{t = 1}^{T} \operatorname{\pmb{x}}_t\operatorname{\pmb{x}}_t'$ denotes the sample covariance matrix based on the observed data. $\operatorname{\pmb{\Sigma}}_{F}$ is the low dimensional covariance matrix of the factors. Within the framework of an AFM, the estimation of a full $\operatorname{\pmb{\Sigma}}_u$ is cumbersome, as the number of parameters to estimate is $\frac{N(N+1)}{2}$ which may exceed the sample size $T$. In order to overcome this problem, we treat $\operatorname{\pmb{\Sigma}}_{u}$ as a diagonal matrix in the first step and define $\operatorname{\pmb{\Phi}}_u = \text{diag}\left(\operatorname{\pmb{\Sigma}}_{u}\right)$ denoting a diagonal matrix that contains only the elements of the main diagonal of $\operatorname{\pmb{\Sigma}}_{u}$. Furthermore, we restrict the covariance matrix of the factors to $\operatorname{\pmb{\Sigma}}_{F} = \operatorname{\pmb{I}}_r$.
Imposing these restrictions has the advantage that the estimation of the covariance matrix of the factors becomes redundant. Hence, our objective function reduces to:
As the true covariance matrix of $\operatorname{\pmb{u}}_t$ allows for correlations of general form, but the previous objective function incorporates the error term structure of a strict factor model, (ref) may be seen as a quasi-likelihood. \citemain{Bai2016a} show that the QML estimator based on (ref) yields consistent parameter estimates. Hence, the consistency of $\operatorname{\pmb{\Phi}}_{u}$ is not affected by the general form of cross-section and serial correlations in $\operatorname{\pmb{u}}_t$.
The factors $\operatorname{\pmb{f}}_t$ can be estimated by generalized least squares (GLS):
where the estimates $\hat{\operatorname{\pmb{\Lambda}}}$ and $\hat{\operatorname{\pmb{\Phi}}}_{u}$ are the ones obtained from the optimization of the objective function in (ref).
The sparse approximate factor (SAF) model allows for sparsity in the factor loadings matrix $\operatorname{\pmb{\Lambda}}$ by shrinking single elements of $\operatorname{\pmb{\Lambda}}$ to zero. This is obtained by the $l_1$-norm penalized MLE of (ref) based on the following optimization problem:
where $\mu \geq 0$ denotes a regularization parameter. Note that the number of factors $r$ is predetermined and assumed to be fixed. Sparsity is obtained by shrinking some elements of $\operatorname{\pmb{\Lambda}}$ to zero, such that not all $r$ factors load on each $x_{it}$. Hence, this framework allows for weak factors (see, e.g. \citemain{Onatski2012}) that affect only a subset of the $N$ time series.
It is well known that the factors and factor loadings in AFM model in (ref) are only identified up to an arbitrary non-singular rotation matrix $\operatorname{\pmb{P}}$. This follows from the fact that $\operatorname{\pmb{X}} = \operatorname{\pmb{\Lambda}} \operatorname{\pmb{P}} \operatorname{\pmb{P}}^{-1'}\operatorname{\pmb{F}}' + \operatorname{\pmb{u}} = \operatorname{\pmb{\Lambda}}^{*} \operatorname{\pmb{F}}^{*'} + \operatorname{\pmb{u}}$, with $\operatorname{\pmb{\Lambda}}^{*} = \operatorname{\pmb{\Lambda}} \operatorname{\pmb{P}}$ and $\operatorname{\pmb{F}}^{*'} = \operatorname{\pmb{P}}^{-1'}\operatorname{\pmb{F}}'$. \\ In contrast to the standard AFM model, which needs additional restrictions to identify $\operatorname{\pmb{P}}$, our SAF model with embedded $l_1$-norm penalty function ensures the identification of the factors and factor loadings up to a unitary generalized permutation matrix $\operatorname{\pmb{P}}$.\footnote{A short demonstration of the fact that $\operatorname{\pmb{P}}$ can only be a unitary generalized permutation matrix for the $l_1$-norm is given in Section (ref) in the Supplement, as well as in \citemain{HornJohnson2012}.}
Hence, by fixing the ordering of columns, e.g. by sorting the columns of the factor loadings matrix according to their respective sparsity, and assuming that the SAF estimator $\hat{\operatorname{\pmb{\Lambda}}}$ has identical column signs as the true factor loadings $\operatorname{\pmb{\Lambda}}_0$, as part of the identification conditions, the SAF model is fully identified. However, it should be noted that the identification of the SAF model only holds if the $l_1$-norm penalty on $\operatorname{\pmb{\Lambda}}$ enters the penalized optimization problem (ref), i.e. for $\mu > 0$. For $\mu = 0$, we are in the standard ML setting for the AFM and solely for this case we identify the model, following \citemain{Lawley1971}, by imposing the identification restriction that $\operatorname{\pmb{\Lambda}}'\operatorname{\pmb{\Phi}}_{u}^{-1}\operatorname{\pmb{\Lambda}}$ is diagonal, with distinct diagonal entries that are arranged in a decreasing order.
In contrast to the weak factor assumption introduced in the following, the pervasiveness assumption conventionally made for standard approximate factor models (e.g. \citemain{Bai2002}, \citemain{stock2002forecasting}), implies that the $r$ largest eigenvalues of $\operatorname{\pmb{\Lambda}}' \operatorname{\pmb{\Lambda}}$ diverge at the rate $\mathcal{O}(N)$. Intuitively, this means that all factors are strong and the entire set of time series is affected. Consequently, the sparsity in the factor loadings matrix introduced in Assumption (ref) below considerably relaxes the conventional pervasiveness assumption.
\footnotetext{The lower limit 1/2 for $\beta$ is necessary to consistently estimate the factors. See \lemref{lem_est_factor} in Section (ref) in the Supplement.} Assumption (ref) implies that the $r$ largest eigenvalues of $\operatorname{\pmb{\Lambda}}' \operatorname{\pmb{\Lambda}}$ diverge with the rate $\mathcal{O}\left( N^{\beta}\right)$, which can be much slower than in the standard AFM. Furthermore, the parameter $\beta$ can take on different values for each of the eigenvalues of $\operatorname{\pmb{\Lambda}}' \operatorname{\pmb{\Lambda}}$. Hence, the eigenvalues can diverge at different rates. On the other hand, the special case of $\beta = 1$, implies the standard AFM framework with strong factors (i.e. \citemain{FanLiaoMincheva2013}, \citemain{Bai2016}). Hence, our sparse approximate factor model offers a convenient generalization of the standard one. Furthermore, Assumption (ref) has a direct implication on the sparsity of $\operatorname{\pmb{\Lambda}}$. In fact, this can be deduced by upper bounding the spectral norm of $\operatorname{\pmb{\Lambda}}$ according to the following expressions:
This result shows that imposing the weak factor assumption limits the amount of affected time series across all factors and hence requires a non-negligible amount of zero elements in each column of the factor loadings matrix. Nevertheless, the number of zero factor loadings can be arbitrarily small as $\beta$ increases. Note, that the lower bound of equation (ref) restricts the number of zero elements in each column of $\operatorname{\pmb{\Lambda}}$, so that we can disentangle the common component from the idiosyncratic one.
The pervasiveness assumption imposed by the standard AFM, further implies a clear separation of the eigenvalues of the data covariance matrix into two groups, corresponding to the diverging eigenvalues of the common component and the bounded eigenvalues of the covariance matrix of the idiosyncratic errors. These characteristics can be observed in Figure (ref), where both panels illustrate the eigenvalue structure of datasets, that are simulated only based on strong factors for $T=450$ and different $N$.
The panels differ solely in the number of factors included, where the left panel includes one strong factor and the right panel depicts the case of four strong factors. Both graphs reveal a clear partition in their respective eigenvalue structures, into sets of eigenvalues that diverge with the sample size $N$ corresponding to the number of included strong factors and sets of bounded eigenvalues associated to the idiosyncratic components.
However, such a clear separation in the eigenvalue structure of the covariance matrix cannot typically be found in real datasets. An example offers a dataset that contains the monthly asset returns of stocks constituents of the S&P 500 stock index available for the entire period of 450 months,\footnote{The same dataset is also used in our empirical application and is described in more detail in Section (ref).} whose eigenvalue distribution is illustrated in Figure (ref). The graph shows a clear distinction between the first eigenvalue and the remaining eigenvalues. However, the remaining eigenvalues diverge at a slower rate and a clear separation between the common and idiosyncratic component as implied by the standard AFM is impossible. Hence, the weak factor framework that allows for a slower divergence rate in the eigenvalues of the common component is more realistic for modeling the eigenvalue structure of real datasets. Furthermore, the weak factor assumption supports the well-documented empirical evidence that the eigenvalues of the sample covariance matrix of asset returns diverge at different rates (see, e.g. \citemain{Ross1976} and \citemain{Trzcinka1986}). Figure (ref) depicts the eigenvalue structure of a dataset, which is generated by one strong factor and three weak factors. This model with weak factors nicely mimics the decaying eigenvalue structure we observe for the S&P 500 asset returns.
In order to relax the imposed diagonality assumption on $\operatorname{\pmb{\Sigma}}_{u}$ in the first step of our estimation, we re-estimate the covariance matrix of the idiosyncratic error term by means of the principal orthogonal complement thresholding (POET) estimator by \citemain{FanLiaoMincheva2013}. The POET estimator is based on soft-thresholding the off-diagonal elements of the sample covariance matrix of the residuals obtained from the estimation of an approximate factor model. Hence, it introduces sparsity in the idiosyncratic covariance matrix and offers a solution to the non-invertibility problem, generated using the sample covariance estimator, especially in high dimensional settings, where $N$ is close or even larger than $T$. More specifically, the estimated idiosyncratic error covariance matrix $\hat{\operatorname{\pmb{\Sigma}}}_{u}^{\tau}$ based on the POET method is defined as:
where $\hat{\sigma}_{u,ij}$ is the $ij$-th element of the sample covariance matrix \\ $\operatorname{\pmb{S}}_u = \frac{1}{T} \sum_{t = 1}^{T} (\operatorname{\pmb{x}}_{t} - \hat{\operatorname{\pmb{\Lambda}}}\hat{\operatorname{\pmb{f}}}_t)(\operatorname{\pmb{x}}_{t} - \hat{\operatorname{\pmb{\Lambda}}}\hat{\operatorname{\pmb{f}}}_t)'$ of the estimated factor model residuals, $\tau = \frac{1}{\sqrt{N}}+\sqrt{\frac{\log(N)}{T}}$ is a threshold\footnote{The threshold $\tau$ is based on the convergence rate of the idiosyncratic error covariance estimator specified in \lemref{lem_idio}. in Section (ref) in the Supplement.} and $\mathcal{S}(\cdot)$ denotes the soft-thresholding operator defined as:
In contrast to \citemain{FanLiaoMincheva2013}, who use the residuals of a static factor model based on the PCA estimator, our estimates are based on the residuals obtained from our sparse factor model.
Thus, we follow the approach of \citemain{FanLiaoMincheva2013} and use a two-step procedure, where at first step we identify the common part, however, unlike in the PCA framework we allow for the weak factors; and in the second step, we model the general covariance structure for the idiosyncratic component. By using a two-step procedure, we control for the sparsity patterns in $\operatorname{\pmb{\Lambda}}$ and $\operatorname{\pmb{\Sigma}}_u$ separately and hence, this ensures that the sparsity in the loadings matrix is not distorted by the sparsity in the idiosyncratic error covariance matrix.
Moreover, the joint estimation of two high-dimensional matrices with embedded $l_1$-norms, would become computationally burdensome and lead to considerable numerical instabilities. By separating the joint estimation into our two-step procedure we obtain a numerical stable optimization method that is computationally time-efficient.
The estimator of the data covariance matrix based on the approximate factor model is obtained according to $ \operatorname{\pmb{\Sigma}} = {\operatorname*{Cov}}\left[ {\operatorname{\pmb{X}}} \right] = \operatorname{\pmb{\Lambda}} \operatorname{\pmb{\Sigma}}_{F} \operatorname{\pmb{\Lambda}}' + \operatorname{\pmb{\Sigma}}_{u}$. Hereby, we first estimate the factors $\operatorname{\pmb{f}}_t$ and the factor loadings $\operatorname{\pmb{\Lambda}}$ according to our sparse factor model introduced in Section (ref). Consistent estimates of $\operatorname{\pmb{\Lambda}}$ and $\operatorname{\pmb{f}}_t$ are obtained by MLE and GLS as given by (ref) and (ref), respectively. This yields the estimates of the common and idiosyncratic components of the AFM defined in (ref). The latter one is used as input to estimate $\operatorname{\pmb{\Sigma}}_u$ by the POET estimator introduced in Section (ref). Hence, our SAF covariance matrix estimator is given by:
where $\operatorname{\pmb{S}}_{\hat{F}}$ denotes the sample estimator for the covariance matrix of the estimated factors, which is positive definite because the number of observations exceeds the number of factors. Further, using the convergence rate of the idiosyncratic error covariance matrix for the threshold $\tau$ also guarantees that $\hat{\operatorname{\pmb{\Sigma}}}_{u}^{\tau}$ is positive definite with probability tending to one according to \citemain{Bickel2008}. Hence, the covariance matrix estimator $\hat{\operatorname{\pmb{\Sigma}}}_\text{SAF}$ is positive definite by construction.
The implementations issues, the choice of the number of factors and the selection of the tuning parameter $\mu$ are described in Section (ref) in the Supplement.
In order to establish the consistency of the factor loadings matrix $\operatorname{\pmb{\Lambda}}$ and the data covariance matrix $\operatorname{\pmb{\Sigma}}$ estimators, we adapt the following standard assumptions:
The assumptions in (ref) impose regularity conditions on the data generating process and are identical to those imposed by \citemain{Bai2016}. Condition (ref) imposes strict stationarity for $\operatorname{\pmb{u}}_t$ and $\operatorname{\pmb{f}}_t$ and requires that both terms are not correlated. Condition (ref) requires exponential-type tails, which allows to use the large deviation theory for $\frac{1}{T} \sum_{t = 1}^{T} u_{it} u_{jt} - \sigma_{u, ij}$ and $\frac{1}{T} \sum_{t = 1}^{T} f_{jt} u_{it}$. In order to allow for weak serial dependence, we impose a strong mixing condition specified in Condition (ref). Further, Condition (ref) implies bounded eigenvalues of the idiosyncratic error covariance matrix, which is a common identifying assumption in the factor model framework.
Assumptions (ref) imposes sparsity conditions on $\operatorname{\pmb{\Lambda}}$ and $\operatorname{\pmb{\Sigma}}_{u}$, where condition (ref) defines the quantity $L_N$ that reflects the number of non-zero elements in the factor loadings matrix $\operatorname{\pmb{\Lambda}}$. As the number of factors $r$ is assumed to be fixed, (ref) restricts the number of non-zero elements in each column of $\operatorname{\pmb{\Lambda}}$ to be upper bounded by $N$. At the same time, this assumption allows for a sparse factor loadings matrix with less than $N$ non-zero elements. Condition (ref) specifies $S_N$ that quantifies the maximum number of non-zero elements in each row of $\operatorname{\pmb{\Sigma}}_{u}$, following the definition of \citemain{Bickel2008}. Furthermore, it restricts the number of zero elements in each row of $\Sigma_u$. Hence, it requires that $\Sigma_u$ is not too dense.
The proof of Theorem (ref) is given in the Sections (ref) and (ref) in the Supplement. Under the given regularity conditions this theorem establishes the average consistency in the Frobenius norm of the estimators for the factor loadings matrix and idiosyncratic error covariance matrix based on our sparse factor model. More specifically, $\operatorname{\pmb{\Lambda}}$ and $\operatorname{\pmb{\Phi}}$ can be estimated consistently, regardless of the diagonality restriction on $\operatorname{\pmb{\Sigma}}_{u}$ in the first step of our estimation procedure. Consequently, the factors $\operatorname{\pmb{f}}_t$ estimated based on GLS are as well consistent. The lower limit $1/2$ on $\beta$ is a necessary condition to achieve consistency. Intuitively this means that the factors should not be too weak such that there is still a clear distinction between the common and idiosyncratic component. Furthermore, the second step estimator of $\operatorname{\pmb{\Sigma}}_u$ can be consistently estimated under the spectral norm.
Finally, in this section we take a closer look on the asymptotic properties of the SAF covariance matrix estimator, given in Section (ref). The following theorem gives the convergence rates of the covariance matrix estimator and of its inverse under different matrix norms.
The proof of Theorem (ref) is given in Section (ref) in the Supplement. Similar as for Theorem (ref), we assume that the regularization parameter $\mu = o(1)$ and $\log(N) = o(T)$. Equation (ref) in Theorem (ref) shows that the covariance matrix estimator based on the sparse factor model in equation (ref) is consistent if we consider the weighted quadratic norm for the entire set of possible values for $\beta$.
Generally, convergence under the average Frobenius norm is hard to achieve because of the too fast diverging eigenvalues of the common component (see \citemain{FanLiaoMincheva2013}). However, according to equation (ref) our SAF covariance matrix estimator is consistent, if $\mu = o\left(N^{-\beta/2}\right)$ and $1/2 \leq \beta \lessapprox 9/10$. Hence, the relaxation of the pervasiveness assumption in the standard approximate factor model to allow for weak factors leads to convergence of the covariance estimator under the average Frobenius norm. The upper bound for $\beta$ follows from the expression $\frac{N^{\beta} \log N^{\beta}}{N}$ in Equation (ref) of Theorem (ref).\footnote{A closed form solution for the upper bound of $\beta$ is not feasible, hence we numerically approximate the maximum value of $\beta$ in the neighbourhood of one such that the expression $\frac{N^{\beta} \log N^{\beta}}{N}$ converges to zero.} Further, Equation (ref) of Theorem (ref) shows that the inverse of $\operatorname{\pmb{\Sigma}}_{\normalfont{\text{SAF}}}$ is consistently estimated under the average Frobenius norm.
In the following, we present Monte Carlo evidence on the finite sample properties of our new covariance estimator. In particular, we focus on the accuracy of the covariance matrix estimates depending on the dimensionality as well as on the strength of correlations in the true covariance matrix to be estimated. The simulation results for the SAF estimator are compared to the ones obtained from eight competing estimators that are popular in the literature.
For our Monte Carlo experiments we use three different designs of the true covariance matrix $\operatorname{\pmb{\Sigma}}$. In the first case, we consider the uniform covariance matrix design used in \citemain{AbadirDistasoZikes2014}, which takes the following form:
where $\mathcal{U}_{(0,1)}$ denotes a standard uniform random variable, and we set $\eta \in \{0.025, 0.05, 0.075\}$. In this setting, $\eta$ controls for the correlations among the variables, where an increase in $\eta$ amplifies the strength of the dependencies among the covariates.
For the second design, we use the sparse covariance matrix suggested by \citemain{BienTibshiraniothers2011}, which contains zero entries for the off-diagonals with a certain probability. More specifically, the $ij$-th element of the covariance matrix $\sigma_{ij} = \sigma_{ji}$ is assigned to be non-zero with probability $p$, where $p \in \{0.05, 0.075, 0.1\}$. Similar as in the uniform design, the diagonal elements are set to 1. The non-zero off-diagonal elements are independently drawn from the uniform distribution $\mathcal{U}_{(0,0.2)}$.
Finally, the last design we consider is based on a generalized spiked covariance model as in \citemain{BaiYao2012}. More precisely, we use the following definition:
where $r_1 - r_4$ correspond to four spiked eigenvalues and $\operatorname{\pmb{\Sigma}}_u$ is a covariance matrix based on the uniform design in equation (ref). As this covariance matrix design complies with the approximate factor model framework, estimation approaches that are based on a factor model specification might benefit from this setting. More precisely, the first part of equation (ref) is in accordance with the eigenvalue distribution of the common component in an AFM with four factors, whereas the second part in (ref) corresponds to the covariance matrix of the idiosyncratic component and allows for weak correlations among the errors. In the simulation, we consider the following specification for the spiked eigenvalues: $r_1 = r_2 = N, r_3 = N^{0.8}, r_4 = N^{0.5}$. This design is in line with the weak factor framework, where the first two factors are strong and the last two correspond to weak factors.
For all three covariance matrix designs, we draw a time independent random data series $\operatorname{\pmb{X}}$ from a multivariate normal distribution with zero population mean.\footnote{The same Monte Carlo experiments are carried out based on data from a multivariate t-distribution with five degrees of freedom. The results are rather similar to the multivariate normal setting and can be obtained upon request.} The time dimension $T$ is set to 60, which relates to a dataset with 5 years of monthly data. The number of replications is 1000. Further, we consider several dimensions for $\operatorname{\pmb{X}}$ and set $N \in \{30, 50, 100, 200\}$. As goodness of fit criterion for the difference between the true and the estimated covariance matrix, we use the Frobenius norm.
Table (ref) gives an overview of the methods for the covariance matrix estimation that are compared in our Monte Carlo experiments. A more detailed description of the alternative strategies is provided in Section (ref) in the Supplement.
Since the observed factor models SIM and FF3F require additional, observable economic factors, they will only be considered in our empirical application to portfolio choice.
Table (ref) below contains the Monte Carlo results for the uniform design of the true covariance matrix, Table (ref) gives the results based on the sparse covariance matrix design, while Table (ref) shows the results for the covariance matrix design with spiked eigenvalues. Interestingly, we find a very similar and clear picture. In terms of the goodness of fit, our sparse approximate factor model approach provides the smallest Frobenius norm, i.e. the SAF fits the true covariance matrix best. These results hold for all of the three rather different designs, all dimensions and degrees of correlation between the variables. Note that the advantage of the SAF model in accurately estimating the true covariance matrix is even more pronounced when $N$ increases, especially for the two high dimensional settings with $N= 100, 200$ and $T = 60$. Concerning the alternative approaches, ST, which is rather similar to our approach, performs second best in most of the scenarios.
However, for small samples ($N = 30, 50$) it is outperformed by LW-NL, for the uniform and sparse covariance matrix designs. Furthermore, for the uniform covariance matrix design for high dimensions and very strong dependencies ($N= 100,200, \eta = 0.075$), ADZ performs slightly better than ST. It is also interesting to note that direct $l_1$-norm penalization of the covariance matrix as suggested by \citemain{ BienTibshiraniothers2011} does not do nearly as well as our approach, which profits from sparsity in the factor loadings matrix and thresholding of the covariance matrix of the idiosyncratic component. Moreover, the results for the POET estimator by \citemain{FanLiaoMincheva2013} that allows only for sparsity in the idiosyncratic error covariance matrix indicate that allowing for sparsity in the factor loadings matrix leads to a considerable improvement in the estimation accuracy.
\FloatBarrier
Empirical portfolio models, particularly when applied to large asset spaces, suffer from a high degree of instability. The estimation of $N $ mean and $N(N + 1)/2$ variance-covariance parameters yields extremely noisy estimates of portfolio weights with large standard errors. It is well-documented that these estimated portfolios show poor out-of-sample performance, extreme short positions and no diversification (e.g. \citemain{Jobson/Korkie1980} and \citemain{Michaud1989}). In order to mitigate these shortcomings and to improve portfolio estimates against extreme estimation noise, a range of alternative strategies have been proposed including the shrinkage estimation of the covariance matrix of asset returns (\citemain{Ledoit2003}; \citemain{LedoitWolf2018} and \citemain{Kourtis2012}).
In the following, we investigate to what extent the SAF model can be used to obtain robustified estimates of high-dimensional covariance matrices of asset returns as input for empirical portfolio models. In an out-of-sample portfolio forecasting experiment, we compare the performance of the global minimum variance portfolio (GMVP) strategy based on a covariance matrix estimated by our sparse factor model to popular alternative portfolio strategies with regularized covariance estimators. As in many other studies, we restrict our analysis to the GMVP, because its vector of portfolio weights, $ \operatorname{\pmb{\omega}} = \frac{\operatorname{\pmb{\Sigma}}^{-1} \mathbf{1}_N}{\mathbf{1}_N' \operatorname{\pmb{\Sigma}}^{-1} \mathbf{1}_N}$, is solely a function of the covariance matrix of the asset returns. Thus, for estimating the GMVP the mean vector of asset returns is redundant and its empirical performance only depends on the quality of the covariance matrix estimator.
In a first step, we theoretically analyze the properties of the GMVP weights based on the SAF estimator. The results are summarized in the following proposition:
The proof is given in Section (ref) in the Supplement. \propref{prop_eig_cov} shows that allowing for sparsity in the factor loadings matrix leads to shrinking the eigenvalues of the precision matrix towards the ones of an identity matrix. Hence, the portfolio weights based on our SAF model are shrunken towards those of the $1/N$ portfolio. This result makes intuitively sense as it is reasonable to invest in the equally weighted portfolio in the case of great estimation instabilities regarding the covariance matrix.
The dataset comprises the monthly excess returns of stocks of the S&P 500 index, that were constituents of the index in December, 2016. The excess returns are obtained by subtracting the corresponding one-month Treasury bill rate from the asset returns. We consider the time period from January, 1980 until December, 2016, which yields $T = 443$ monthly returns for each of the 205 available stocks.\footnote{The return data are retained from Thompson Reuters Datastream.} In order to check the performance of our estimator with respect to the dimensionality of the asset space, we consider the following portfolio sizes: $N \in \{30, 50, 100, 200\}$. Out of the 205 stocks, we select at random individual subsets from the overall number of assets and work with the selected assets for the entire forecasting experiment.
Since by construction, a theoretical portfolio built on a subset of assets from a larger portfolio cannot outperform the larger one, an observed inferiority of the larger empirical portfolio can only be the consequence of higher estimation noise due to the larger dimensionality, which overcompensates for the ex-ante theoretical superiority. Therefore, this selection strategy provides us with insights into the impact of estimation noise on the performance of empirical portfolios.
In order to estimate the portfolio weights for each strategy, we apply a rolling window approach with $h=60$ months, corresponding to 5 years of historic data. Thus, at time $t$ we use the last 60 months from $t-59$ until $t$ for our estimation. Using the estimated portfolio weights, we compute the out-of-sample portfolio return $\hat{r}^p_{t+1}(s) = \hat{\operatorname{\pmb{\omega}}}(s)' \pmb{r}_{t+1}$ for the period $t+1$ for the 12 different estimation strategies $s=1, \ldots, 12$. All portfolios are rebalanced on a monthly basis. This generates a series of $T-h$ out-of-sample portfolio returns. The results are then used to estimate the mean $\mu(s)$ and variance $\sigma^2(s)$ of the portfolio returns for each strategy by their empirical counterparts:
We repeat this procedure 100 times to avoid that the out-of-sample results depend on the initially randomly selected stocks. Hence, all results reported below are average outcomes across the 100 forecasting experiments.
The criteria for the performance evaluation are the out-of-sample standard deviation (SD), the average return (AV), the Certainty Equivalent (CE) and the Sharpe ratio (SR).
Table (ref) contains the annualized results of our comparative study on the out-of-sample performance of different portfolio estimation approaches. The results represent average outcomes across the 100 different forecasting experiments for each of the four performance measures. Our sparse approximate factor model (SAF) yields the lowest out-of-sample portfolio standard deviation for all portfolio dimensions, i.e. it is performing best for the performance criterion the GMVP-strategy is designed for.
In theory, the GMVP-strategy may not necessarily outperform the $1/N$-strategy in terms of the remaining three performance criteria, since it completely disregards optimization with respect to the expected portfolio return. Nevertheless, our SAF model also outperforms the $1/N$-strategy and the other estimation approaches in terms of AV, CE and SR, which depend on the expected return. In the portfolio forecasting experiment, our regularization method does best for the expected out-of-sample portfolio return.
It is of utmost importance to note that the superiority of our approach does not only hold for different performance measures, but also for all portfolio dimensions. The SAF model performs best for low, but also for high dimensional portfolios, for which the sample size is much smaller than the portfolio dimension, i.e. $T \ll N$. This indicates, at least for this specific application, that the selection of the penalty parameter is reasonable.
As mentioned earlier, increasing the portfolio dimension does not necessarily improve the out-of-sample performance of an empirical portfolio as the theoretical gains maybe overcompensated by the increase in estimation noise due to the increase in the number of parameters to be estimated. It is not too surprising that this phenomenon is most dramatically pronounced for the plug-in estimator of the GMVP, but we also find it to some extent for the DFM. Moreover, for the SIM and FF3F, we do not find a strict monotonicity between portfolio dimension and portfolio performance, while the performance of our SAF model strictly increases with $N$.
While in the portfolio forecasting experiment for any performance measure and any portfolio dimension our sparse factor model shows the best performance, there is no clear further ranking regarding the other approaches. FF3F is performing second best in terms of the minimization of portfolio risk for all portfolio dimensions, but it is outperformed by other estimation approaches when performance measures other than the portfolio risk are considered.
Our comparative study also confirms the findings of \citemain{DeMiguel2009a} that the $1/N$ portfolio is a strong competitor for many alternative portfolio strategies. For low dimensions ($N = 30$ and $N = 50$), we can see that, apart from our estimator only the single factor model generates a higher average SR compared to the equally weighted portfolio, although it is very close to it. In terms of the portfolio risk, only our method and FF3F reveal performance superior to the $1/N$ portfolio for low dimensions of the asset space. The picture slightly changes, when higher asset dimensions ($N > 50$) are considered. For higher dimensions, the method by \citemain{AbadirDistasoZikes2014} is a serious competitor to the $1/N$ portfolio. This mirrors our finding from the simulation study in Section (ref), where the ADZ estimator performs comparatively well in high dimensional settings with strong linear dependencies.
Table (ref) in Appendix (ref) provides additional insights into the quality of the weight estimates. The summary statistics indicate that the outstanding performance of the SAF model results from effectively stabilizing the estimated portfolio weights by avoiding extreme positions (moderate minima and maxima in the weight estimates) and by the low standard deviations. Furthermore, the results show that the weights of our SAF estimator shrink towards the weights of the equally weighted portfolio as $N$ increases. This is in line with the theoretical results in Proposition (ref). The relative good performance of SIM and FF3F result from very low variation in the portfolio weights, which come for the SIM with $N=200$ close to the constant weights of the equally weighted portfolio.
In order to check the robustness of our findings, which are based on data from January 1980 until December 2016, we also consider forecasts based on subperiods. We restrict our attention to the standard deviation of the out-of-sample portfolio returns and consider how a gradual increase of the evaluation sample affects the performance of the competing estimators. The results are illustrated in Figure (ref) in Appendix (ref), where the portfolio standard deviation at time $t$ incorporates the out-of-sample portfolio returns until $t$ (e.g. the out-of-sample portfolio standard deviation in January 2005 incorporates the out-of-sample portfolio returns from January 1985 until January 2005). Special attention is given to the periods before and after the financial crisis in 2007. The graphs indicate that the SAF estimator also provides for different subperiods the lowest portfolio standard deviation compared to FF3F and LW-NL. Note, that the difference is more pronounced when the recent financial crisis period is included. Hence, in comparison to our SAF model both, FF3F and LW-NL, fail to pick up the changing risk during the crisis and, as a result, they provide more volatile portfolio estimates.
In this paper, we propose a novel approach for the estimation of high-dimensional covariance matrices based on a sparse approximate factor model. The estimator allows for sparsity in the factor loadings matrix by shrinking single elements of the factor loadings matrix to zero. Hence, this setting reduces the number of parameters to be estimated and therefore leads to a reduction in estimation noise. Furthermore, the sparse factor model framework allows for weak factors, which only affect a subset of the available time series. Thus, our framework offers a convenient generalization to the pervasiveness assumption in the standard approximate factor model that solely leads to strong factors.
We prove average consistency under the Frobenius norm for the factor loadings matrix estimator and consistency in the spectral norm for the idiosyncratic component covariance matrix estimator based on our sparse approximate factor model. The factors estimated using the GLS method are also shown to be consistent. Furthermore, we derive average consistency for our factor model based covariance matrix estimator under the Frobenius norm for a particular rate of divergence for the eigenvalues of the covariance matrix corresponding to the common component. To the best of our knowledge, this result has not been shown in the existing literature because of the fast diverging eigenvalues. Additionally, we provide consistency results of our covariance matrix estimator under the weighted quadratic norm.
In our Monte Carlo study, we analyze the finite sample properties of our covariance matrix estimator for different simulation designs for the true underlying covariance matrix. The results show that our estimator offers the lowest difference in Frobenius norm to the true covariance matrix compared to the competing estimators. Further, the benefit of the covariance matrix estimator based on our sparse factor model is even more pronounced if the dimensionality of the problem increases.
In an out-of-sample portfolio forecasting experiment, we compare the performance of the global minimum variance portfolio based on the covariance matrix estimator of our sparse approximate factor model to alternative estimation approaches frequently used in the literature. The forecasting results reveal that our estimator yields the lowest average out-of-sample portfolio standard deviation across different portfolio dimensions. At the same time, it generates the highest Certainty Equivalent and Sharpe Ratio compared to all considered portfolio strategies. The performance gains of our SAF model are especially pronounced during the recent financial crisis. Hence, our estimator has a stabilizing impact on the portfolio weights, especially during highly volatile periods.
The results of our out-of-sample portfolio forecasting study show a substantial reduction of the portfolio standard deviation of the dynamic factor model compared to the standard approximate factor model, especially for small asset dimensions. Hence, it would be interesting to analyze if a possible extension of our SAF model by considering dynamic factors, would as well lead to a more efficient estimation of the covariance matrix. We leave this for future research.