Sparse Approximate Factor Estimation for High-Dimensional Covariance Matrices

We propose a novel estimation approach for the covariance matrix based on the $l_1$-regularized approximate factor model. Our sparse approximate factor (SAF) covariance estimator allows for the existence of weak factors and hence relaxes the pervasiveness assumption generally adopted for the standard approximate factor model. We prove consistency of the covariance matrix estimator under the Frobenius norm as well as the consistency of the factor loadings and the factors. Our Monte Carlo simulations reveal that the SAF covariance estimator has superior properties in finite samples for low and high dimensions and different designs of the covariance matrix. Moreover, in an out-of-sample portfolio forecasting application the estimator uniformly outperforms alternative portfolio strategies based on alternative covariance estimation approaches and modeling strategies including the $1/N$-strategy.


Introduction
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 → ∞, 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 highdimensional 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 Fan, Liao, and Liu (2016) 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 1regularization 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 Bickel andLevina (2008a, 2008b), Cai and Liu (2011) and Cai and Zhou (2012). 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 Fan, Fan, and Lv (2008) 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 Fan, Liao, and Mincheva (2011) and Fan, Liao, and Mincheva (2013), who propose a covariance estimator based on an approximate factor model representation. While Fan, Liao, and Mincheva (2011) shrink the entries of the covariance matrix of the idiosyncratic errors to zero using the adaptive thresholding technique by Cai and Liu (2011), the approach proposed in Fan, Liao, and Mincheva (2013) 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 O(N)).
The recent paper by Fan, Liu, and Wang (2018) 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. Fan, Liao, and Mincheva (2013)). 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 2 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 3 discusses the theoretical setup and provides the convergence results. In Section 4, we present Monte-Carlo evidence on the finite sample properties of our new covariance estimator, while in Section 5 we show the performance of our approach when applied to the empirical portfolio allocation problem. Section 6 summarizes the main findings and gives an outlook on future research.
Throughout the paper we will use the following notation: π max (A A A) and π min (A A A) are the maximum and minimum eigenvalue of a matrix A A A. Further, A A A , A A A F and A A A 1 denote the spectral, Frobenius and the l 1 -norm of A A A, respectively. They are defined as

The Approximate Factor Model
The following analysis is based on the approximate factor model (AFM) proposed by Chamberlain and Rothschild (1983) 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, . . . , N and t = 1, . . . , 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 λ λ λ i is a (r × 1)-dimensional vector of factor loadings for variable i and f f f t is a (r × 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 variablespecific shocks, which are not captured by the common component λ λ λ ′ i f f 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, Σ Σ Σ u = Cov (u 1t , u 2t , . . . u N t ) ′ .
In matrix notation, (1) can be written as: where X X X denotes a (N × T ) matrix containing T observations for N weakly stationary time series. It is assumed that the time series are demeaned and standardized.
is referred to as a (T × r)-dimensional matrix of unobserved factors, Λ Λ Λ = (λ λ λ 1 , . . . , λ λ λ N ) ′ is a N × r matrix of corresponding factor loadings and u u u is a (N × T )-dimensional matrix of idiosyncratic shocks.
There are several estimation approaches for a factor model as given by (2). The principal component analysis (PCA) a and the quasi-maximum likelihood estimation (QMLE) under normality (see i.e. Bai and Li (2016)) 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 Σ Σ Σ u can be estimated jointly, while PCA-based second stage estimates of Σ Σ Σ u require consistent estimation of Λ Λ Λ and F F F in the first stage.
This, however, may be problematic for the case of a relatively small N, because F F F can no longer be estimated consistently (Bai and Liao (2016)).
a See e.g., Bai and Ng (2002) for a detailed treatment of the PCA in approximate factor models.
The negative quasi log-likelihood function for the data in the AFM is defined as: where S S S x = 1 T T t=1 x x x t x x x ′ t denotes the sample covariance matrix based on the observed data. Σ Σ Σ F is the low dimensional covariance matrix of the factors. Within the framework of an AFM, the estimation of a full Σ Σ Σ u is cumbersome, as the number of parameters to estimate is N (N +1) 2 which may exceed the sample size T . In order to overcome this problem, we treat Σ Σ Σ u as a diagonal matrix in the first step and define Φ Φ Φ u = diag (Σ Σ Σ u ) denoting a diagonal matrix that contains only the elements of the main diagonal of Σ Σ Σ u . Furthermore, we restrict the covariance matrix of the factors to Σ Σ Σ F = I I 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 u u u t allows for correlations of general form, but the previous objective function incorporates the error term structure of a strict factor model, (4) may be seen as a quasi-likelihood. Bai and Li (2016) show that the QML estimator based on (4) yields consistent parameter estimates. Hence, the consistency of Φ Φ Φ u is not affected by the general form of cross-section and serial correlations in u u u t .
The factors f f f t can be estimated by generalized least squares (GLS): where the estimatesΛ Λ Λ andΦ Φ Φ u are the ones obtained from the optimization of the objective function in (4).

The Sparse Approximate Factor Model
The sparse approximate factor (SAF) model allows for sparsity in the factor loadings matrix Λ Λ Λ by shrinking single elements of Λ Λ Λ to zero. This is obtained by the l 1 -norm penalized MLE of (4) based on the following optimization problem: where µ ≥ 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 Λ Λ Λ to zero, such that not all r factors load on each x it . Hence, this framework allows for weak factors (see, e.g. Onatski (2012)) that affect only a subset of the N time series.
It is well known that the factors and factor loadings in AFM model in (2) are only identified up to an arbitrary non-singular rotation matrix P P P . This follows from the fact that X X X = Λ Λ Λ P P P P P P In contrast to the standard AFM model, which needs additional restrictions to identify P P 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 P P P . b 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 Λ Λ Λ has identical column signs as the true factor loadings Λ Λ Λ 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 Λ Λ Λ enters the penalized optimization problem (6), i.e. for µ > 0. For µ = 0, we are in the standard ML setting for the AFM and solely for this case we identify the model, following Lawley and Maxwell b A short demonstration of the fact that P P P can only be a unitary generalized permutation matrix for the l 1 -norm is given in Section S.1.4 in the Supplement, as well as in Horn and Johnson (2012).
In contrast to the weak factor assumption introduced in the following, the pervasiveness assumption conventionally made for standard approximate factor models (e.g. Bai and Ng (2002), Stock and Watson (2002)), implies that the r largest eigenvalues of Λ Λ Λ ′ Λ Λ Λ diverge at the rate 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 2.1 below considerably relaxes the conventional pervasiveness assumption. There exists a constant c > 0 such that, for all N, Assumption 2.1 implies that the r largest eigenvalues of Λ Λ Λ ′ Λ Λ Λ diverge with the rate O N β , which can be much slower than in the standard AFM. Furthermore, the parameter β can take on different values for each of the eigenvalues of Λ Λ Λ ′ Λ Λ Λ. Hence, the eigenvalues can diverge at different rates. On the other hand, the special case of β = 1, implies the standard AFM framework with strong factors (i.e. Fan, Liao, and Mincheva (2013), Bai and Liao (2016)). Hence, our sparse approximate factor model offers a convenient generalization of the standard one. Furthermore, Assumption 2.1 has a direct implication on the sparsity of Λ Λ Λ. In fact, this can be deduced by upper bounding the spectral norm of Λ Λ Λ 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 β increases. Note, that the lower bound of equation (7) restricts the number of zero elements in each column of Λ Λ Λ, 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 1, where both panels illustrate the eigenvalue structure of datasets, that are simulated only based on strong factors for T = 450 and different N. 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, d whose eigenvalue distribution is illustrated in Figure 2. 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 d The same dataset is also used in our empirical application and is described in more detail in Section 5.
evidence that the eigenvalues of the sample covariance matrix of asset returns diverge at different rates (see, e.g. Ross (1976) and Trzcinka (1986)). Figure 3 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.

Estimation of the idiosyncratic error covariance matrix Σ Σ Σ u
In order to relax the imposed diagonality assumption on Σ Σ Σ 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 Fan, Liao, and Mincheva (2013). 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Σ Σ Σ τ u based on the POET method is defined as: whereσ u,ij is the ij-th element of the sample covariance matrix T is a threshold e and S(·) denotes the soft-thresholding operator defined as: In contrast to Fan, Liao, and Mincheva (2013), 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 Fan, Liao, and Mincheva (2013) 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 Λ Λ Λ and Σ Σ Σ 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 1norms, 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.

SAF covariance matrix estimation
The estimator of the data covariance matrix based on the approximate factor model is obtained according to Σ Σ Σ = Cov [X X X] = Λ Λ Λ Σ Σ Σ F Λ Λ Λ ′ + Σ Σ Σ u . Hereby, we first estimate the factors f f f t and the factor loadings Λ Λ Λ according to our sparse factor model introduced in Section 2.2. Consistent estimates of Λ Λ Λ and f f f t are obtained by MLE and GLS as given by (4) and (5), respectively. This yields the estimates of the common and idiosyncratic components e The threshold τ is based on the convergence rate of the idiosyncratic error covariance estimator specified in Lemma S.1.10. in Section S.1.2 in the Supplement.
of the AFM defined in (1). The latter one is used as input to estimate Σ Σ Σ u by the POET estimator introduced in Section 2.3. Hence, our SAF covariance matrix estimator is given where S S SF 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 τ also guarantees thatΣ Σ Σ τ u is positive definite with probability tending to one according to Bickel and Levina (2008a). Hence, the covariance matrix estimatorΣ Σ Σ SAF is positive definite by construction.
The implementations issues, the choice of the number of factors and the selection of the tuning parameter µ are described in Section S.2 in the Supplement.

Large Sample Properties
In order to establish the consistency of the factor loadings matrix Λ Λ Λ and the data covariance matrix Σ Σ Σ estimators, we adapt the following standard assumptions: Assumption 3.1 (Data generating process).
The assumptions in 3.1 impose regularity conditions on the data generating process and are identical to those imposed by Bai and Liao (2016). Condition (i) imposes strict stationarity for u u u t and f f f t and requires that both terms are not correlated. Condition (ii) requires exponential-type tails, which allows to use the large deviation theory for 1 T T t=1 u it u jt −σ u,ij and 1 T T t=1 f jt u it . In order to allow for weak serial dependence, we impose a strong mixing condition specified in Condition (iii). Further, Condition (iv) implies bounded eigenvalues of the idiosyncratic error covariance matrix, which is a common identifying assumption in the factor model framework.
where 1l {·} defines an indicator function that is equal to one if the boolean argument in braces is true, and µ denotes the regularization parameter.
Assumptions 3.2 imposes sparsity conditions on Λ Λ Λ and Σ Σ Σ u , where condition (i) defines the quantity L N that reflects the number of non-zero elements in the factor loadings matrix Λ Λ Λ. As the number of factors r is assumed to be fixed, (i) restricts the number of non-zero elements in each column of Λ Λ Λ 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 (ii) specifies S N that quantifies the maximum number of non-zero elements in each row Under Assumptions 2.1, 3.1 and 3.2 the sparse factor model in (6) satisfies the following properties, as T and N → ∞ and for 1/2 ≤ β ≤ 1: Hence, for log(N) = o(T ) and the regularization parameter µ = o(1), we have: For the covariance matrix estimator of the idiosyncratic errors in the second step, specified in Section 2.3, we get: Hence, for S 2 N d T = o(1) and S N µ = o(1), this yields: The proof of Theorem 3.1 is given in the Sections S.1.1 and S.1.2 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, Λ Λ Λ and Φ Φ Φ can be estimated consistently, regardless of the diagonality restriction on Σ Σ Σ u in the first step of our estimation procedure. Consequently, the factors f f f t estimated based on GLS are as well consistent. The lower limit 1/2 on β 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 Σ Σ Σ u can be consistently estimated under the spectral norm.

Consistency of the Covariance Matrix Estimator
Finally, in this section we take a closer look on the asymptotic properties of the SAF covariance matrix estimator, given in Section 2.4. The following theorem gives the convergence rates of the covariance matrix estimator and of its inverse under different matrix norms.

Theorem 3.2 (Convergence Rates for the Covariance Matrix Estimator).
Under Assumptions 2.1, 3.1 and 3.2, the covariance matrix estimator based on the SAF model in equation (9) satisfies the following properties, as T , N → ∞ and 1/2 ≤ β ≤ 1: denotes the weighted quadratic norm introduced by Fan, Fan, and Lv (2008).
The proof of Theorem 3.2 is given in Section S.1.3 in the Supplement. Similar as for Theorem 3.1, we assume that the regularization parameter µ = o(1) and log(N) = o(T ).
Equation (10) in Theorem 3.2 shows that the covariance matrix estimator based on the sparse factor model in equation (9) is consistent if we consider the weighted quadratic norm for the entire set of possible values for β.
Generally, convergence under the average Frobenius norm is hard to achieve because of the too fast diverging eigenvalues of the common component (see Fan, Liao, and Mincheva (2013)). However, according to equation (11) our SAF covariance matrix estimator is consistent, if µ = o N −β/2 and 1/2 ≤ β 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 β follows from the expression N β log N β N in Equation (11)

Monte Carlo Evidence
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.

Monte Carlo Designs
For our Monte Carlo experiments we use three different designs of the true covariance matrix Σ Σ Σ. In the first case, we consider the uniform covariance matrix design used in f A closed form solution for the upper bound of β is not feasible, hence we numerically approximate the maximum value of β in the neighbourhood of one such that the expression N β log N β N converges to zero. Abadir, Distaso, andŽikeš (2014), which takes the following form: where U (0,1) denotes a standard uniform random variable, and we set η ∈ {0.025, 0.05, 0.075}.
In this setting, η controls for the correlations among the variables, where an increase in η amplifies the strength of the dependencies among the covariates.
For the second design, we use the sparse covariance matrix suggested by Bien and Tibshirani (2011), which contains zero entries for the off-diagonals with a certain probability.
More specifically, the ij-th element of the covariance matrix σ ij = σ ji is assigned to be nonzero with probability p, where p ∈ {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 U (0,0.2) .
Finally, the last design we consider is based on a generalized spiked covariance model as in Bai and Yao (2012). More precisely, we use the following definition: where r 1 − r 4 correspond to four spiked eigenvalues and Σ Σ Σ u is a covariance matrix based on the uniform design in equation (13). 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 (14) is in accordance with the eigenvalue distribution of the common component in an AFM with four factors, whereas the second part in (14) 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 X X X from a multivariate normal distribution with zero population mean. g 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 X X X and set N ∈ {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 1 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 S.3 in the Supplement.

Alternative covariance estimation strategies
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 3 gives the results based on the sparse covariance matrix design, while Table 4 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 g 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 sparse covariance matrix estimator by Bien and Tibshirani (2011) 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.

Simulation results
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, η = 0.075), ADZ performs slightly better than ST. It is also interesting to note that direct l 1 -norm penal- ization of the covariance matrix as suggested by Bien and Tibshirani (2011) 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 Fan, Liao, and Mincheva (2013) 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.

An Application to Portfolio Choice
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. Jobson andKorkie (1980) andMichaud (1989)). 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 (Ledoit and Wolf (2003); Ledoit and Wolf (2018) and Kourtis, Dotsis, and Markellos (2012)).
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, 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: Proposition 5.1. Based on the general definition of the covariance matrix of an approxi-mate factor model given in Section 2.4, we obtain: The proof is given in Section S.1.5 in the Supplement. Proposition 5.1 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.

Data and Design of the Forecasting Experiment
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 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 returnr p t+1 (s) =ω ω ω(s) ′ r r r t+1 for the period t + 1 for the 12 different estimation strategies s = 1, . . . , 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 µ(s) and variance σ 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). 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.

Out-of-Sample Portfolio Performance
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 ≪ 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 DeMiguel, Garlappi, and Uppal (2009) 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 Abadir, Distaso, andŽikeš (2014) is a serious competitor to the 1/N portfolio. This mirrors our finding from the simulation study in Section 4, where the ADZ estimator performs comparatively well in high dimensional settings with strong linear dependencies. Table 6 in Appendix A provides additional insights into the quality of the weight es-

Conclusions
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 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.

Supplement
The supplementary material provides the omitted proofs and additional related results.   Define the penalized log-likelihood Therefore, we can see that equation (S.1) can be written as Define the set, Further, Φ Φ Φ u = diag (Σ Σ Σ u ) and denotes a covariance matrix that contains only the elements of the main diagonal of Σ Σ Σ u .
We impose the following sparsity assumptions on Λ Λ Λ and Σ Σ Σ u : where 1l {·} defines an indicator function that is equal to one if the boolean argument in braces is true, and µ denotes the regularization parameter.
We introduce a lemma that will be necessary for the forthcoming derivations.

Proof. See Lemmas A.3 and B.1 in
Proof. The unpenalized log-likelihood can be decomposed in a similar fashion as in Lemma A.2. in Bai and Liao (2016).
The first term in equation (S.3) can be written as: Hence, we have Now, we consider the second term 1 N tr S S S x Λ Λ Λ Λ Λ Λ ′ + Σ Σ Σ u −1 . Hereby, S S S x is defined as: By the matrix inversion formula we have: Hence, we get: Subsequently, we look at the terms A 1 − A 4 , respectively.
Since π max (Σ Σ Σ u ) and π −1 min Λ Λ Λ ′ Λ Λ Λ are bounded from above uniformly in Ψ δ , we can derive the following expressions similarly as in Bai and Liao (2016): By applying the matrix inversion formula we have, where the second term can be bounded using (S.6) and (S.7), by the following: Hence, In the following, we define s i (A A A) as the i-th singular value of a (m × n) matrix A A A. Furthermore, s max (A A A) denotes the largest singular value of A A A. Using Lemma S.1.1., Fact 9.14.3 and Fact 9.14.23 in Bernstein (2009) and the fact that we have: By the matrix inversion formula, we have for A 4 the following: From equations (S.6) and (S.7), we see that the second term on the right hand side is uniformly of a smaller order than the first term. The first term of A 4 is bounded by: Hence, we can bound the unpenalized log-likelihood function by: By the definition of Q 3 (Λ Λ Λ, Σ Σ Σ u ) we have: Proof. If we consider equation (S.2) at the true parameter values, we get: Hence, by (S.1) and (S.8), we have: Therefore, by Lemma S.1.2. we have: Proof. By equation (S.9) and the definition of Q 1 (Λ Λ Λ,Σ Σ Σ u ) and Q 2 (Λ Λ Λ,Σ Σ Σ u ), we get: where B 1 and B 2 are defined as: By equation (S.10), we can see that where Φ Φ Φ u = diag (Σ Σ Σ u ) and denotes a covariance matrix that contains only the elements of the main diagonal of Σ Σ Σ u . Using the same argument as in the proof of Lemma B.1. in Bai and Liao (2016), we get: Hence, we can conclude the proof by the following derivation: In the following, we establish the consistency of the factor loadings estimator. Initially, we bound the first part of B 2 defined in equation (S.10).
The consistency result forΛ Λ Λ is summarized in the following lemma.
Proof. If we consider equation (S.10), Lemma S.1.4. and Lemma S.1.5., we have Proof. By the definition of the factor estimator in equation (5) we have: As L N = O N β , the first term on the right-hand side can be bounded by: Now, we are going to bound the second term on the right-hand side of (S.12). For this we first analyze the termΛ Λ Λ Using Lemma S.1.6., the first term can be bounded by: The second term can be bounded using Lemma S.1.4. according to: Hence, using (S.13), (S.14) and (S.15) yields: we have by Lemma S.1.6. and Lemma S.1.7.: Proof. Based on Lemma A.3.(iii) by Fan, Liao, and Mincheva (2011) we have: where the authors show that the first term on the right-hand side is O p

S.1.2 Rate of convergence for the idiosyncratic error covariance matrix estimator
In what follows, we are going to determine the convergence rate of the idiosyncratic error covariance matrix estimator based on soft-thresholding.

S.1.3 Convergence Rates for the Covariance Matrix Estimator
Proof: Theorem 3.2 (Convergence Rates for the Covariance Matrix Estimator) whereΣ Σ Σ τ u corresponds to the POET estimator of Fan, Liao, and Mincheva (2013). Similar as in Fan, Liao, and Mincheva (2013), we consider the weighted quadratic norm introduced by Fan, Fan, and Lv (2008) and which is defined as: Proof. The weighted quadratic norm of the difference between the estimated covariance matrixΣ Σ Σ SAF and the true one Σ Σ Σ can be expressed as: If we consider C C C =Λ Λ Λ − Λ Λ Λ 0 we can introduce the following definitions: Using the previous definitions, we can rewrite the first term in (S.19) in the following form Hence, equation (S.19) can be expressed as: (S.20) Now we analyze each term in (S.20) separately: Hence, by Lemma S.1.10. we get: Under the Frobenius norm we have: Lemma S.1.12.
Proof. A similar argument as in Lemma S.1.11 leads to: where the second term can be bounded by Furthermore, the first term in (S.21) has the same upper bound. Hence, again by using Lemma S.1.10 we get:

Inverse of the covariance matrix
Define,Ĝ Proof.
(ii) By the definition ofĜ G G and G G G 0 , we have: Hence, the previous quantitiy can be decomposed according to: If we bound all three terms on the right hand side of equation (S.22), we get: Lemma S.1.14.
Proof. Using the Sherman-Morrison-Woodbury inverse formula, we get In the following, we bound each of the six terms, separately. Further, Hence, also by Lemma S.1.13. (i) Similarly, L 5 = O p (L 4 ). Finally, Adding up the terms L 1 − L 6 gives

S.1.4 Isometry group for the l p -norm
In this section, we provide a short demonstration that only a unitary generalized permutation matrix can be an isometry for the l p -norm.
Define the l p norm of a (n × n)-dimensional matrix A A A as: Consider the isometry P P P , which is (n × n)-dimensional for the l p norm, with 1 ≤ p < 2 < q ≤ ∞ and the (n × n)-dimensional isometry P P P * for the l q norm. For each standard basis vector e j , we have: | P P P e j | p = | e j | p = 1.
Hence, n i=1 |p ij | p = 1, for each j = 1, . . . , n. Further, we have |p ij | ≤ 1, for all i, j = 1, . . . , n and n i=1 n j=1 |p ij | p = n. Similar, considering P P P * and the l q norm, yields n i=1 n j=1 |p ij | q = n, as well. However, |p ij | q ≤ |p ij | p holds with equality if and only if p ij = 0, 1 or −1. Hence, each column of P P P contains only one non-zero element and it has unit modulus. Furthermore, the non-singularity of P P P ensures that each row of P P P has no more than one non-zero element.
S.1.5 Proof of Proposition 5.1 Proof. As presented in Section 2.4, the general equation of the covariance matrix estimator based on an approximate factor model is given by: Correspondingly, the precision matrix is given by the inverse of the two matrices on the right-hand side of (S.23): Further, by the definition of the factor loadings matrix Λ Λ Λ, the first part on the right-hand side of (S.23) can be expressed as: . . .
where C C C denotes the upper and lower diagonal block of the matrix Λ Λ Λ Λ Λ Λ ′ .
Hence, the sum of the eigenvalues of Λ Λ Λ Λ Λ Λ ′ is calculated as: From equation (S.25), we can clearly see that sparsity or zeros in the factor loadings matrix corresponds to shrinking the sum of the eigenvalues of the covariance of the common component.
In the next step, we want to analyze the global minimum variance portfolio weights based on the estimate of the covariance matrix of our SAF model.
Without loss of generality we assume that the idiosyncratic error covariance matrix is an identity matrix, which corresponds to a high penalization of the off-diagonal elements based on the POET method. Hence, the precision matrix in (S.24) simplifies to: In the following, we have a look at the eigenvalues of the precision matrix of our SAF estimator based on equation (S.26): From equations (S.25) and (S.27), we can see that the possible sparsity in Λ Λ Λ allowed by our SAF model shrinks the precision matrix based on the SAF model towards an identity matrix.
As the GMVP weights directly depend on an estimate of the precision matrix this implies a shrinkage of the SAF portfolio weights towards the weights of the 1/N-portfolio.

S.2 Implementation Issues
For the implementation of the SAF model, we use a two-step estimation procedure that treats Σ Σ Σ u in the first step as a diagonal matrix, denoted as Φ Φ Φ u and re-estimates the idiosyncratic error covariance matrix in a second step by the method introduced in Section 2.3.
Theorem 3.1 shows that this two-step procedure yields consistent estimates for Λ Λ Λ and Σ Σ Σ u .

S.2.1 Majorized Log-Likelihood Function
The numerical minimization of the objective function (6) is cumbersome as it is not globally convex. This problem arises because the first term in (6) log det Λ Λ Λ Λ Λ Λ ′ + Φ Φ Φ u is concave in Λ Λ Λ and Φ Φ Φ u , whereas the second term tr S S S x Λ Λ Λ Λ Λ Λ ′ + Φ Φ Φ u −1 is convex. For our implementation we employ the majorize-minimize EM algorithm introduced by Bien and Tibshirani (2011). The idea of this optimization approach is to approximate the numerically unstable concave part log det Λ Λ Λ Λ Λ Λ ′ + Φ Φ Φ u by its tangent plane, which corresponds to the following expression: where the subscript m denotes the m-th step in an iterative procedure outlined in Section S.2.2. Replacing the concave part in (4) by the convex expression in (S.28), yields the following majorized log-likelihood function: Augmenting the majorized log-likelihood by the l 1 -penalty term, leads to the following optimization problem for our SAF model: As all three components in (S.30) are convex, the optimization problem simplifies considerably compared to the original problem in equation (6).

S.2.2 Projection Gradient Algorithm
In order to minimize (S.30) efficiently, we apply the fast projected gradient algorithm proposed by Bien and Tibshirani (2011). More specifically, we approximate the majorized log-likelihoodL m (Λ Λ Λ) in (S.29) by the following expression: where t is the depth of projection i and which corresponds to the first derivative ofL(Λ Λ Λ) with respect to Λ Λ Λ. Hence, our final optimization problem corresponds to: i We set t = 0.01 for all our applications below.
The minimization of the objective function (S.32) can be carried out by computing its gradient with respect to λ ik and setting it to zero, which yields: where ν ik denotes the subgradient of |λ ik |. Solving (S.33) for a specificλ ik , gives: where S denotes the soft-thresholding operator defined in equation (8). Equation (S.34) can be used to update the estimated factor loadingsλ ik,m+1 given the estimate from the previous stepλ ik,m .
In order to obtain an update for the estimate of the covariance matrix of the idiosyncratic errors Φ Φ Φ u , we use the EM algorithm suggested by Bai and Li (2012): Our iterative estimation procedure for the SAF model can be briefly summarized as given below.

Iterative Algorithm
Step 1: Obtain an initial consistent estimate for the factor loadings matrix Λ Λ Λ and for the diagonal idiosyncratic error covariance matrix Φ Φ Φ u , i.e. by using unpenalized MLE and set m = 1.
Step 6: Re-estimate the covariance matrix of the idiosyncratic errors based on the procedure introduced in Section 2.3.
For the high dimensional case of N > T , the sample covariance matrix S S S x is not of full rank and hence leads to inconsistent parameter estimates. To overcome this problem, we adopt the solution proposed by Bien and Tibshirani (2011), who suggest augmenting the diagonal elements of S S S x by an arbitrarily small ε > 0, when S S S x is not of full rank. This augmentation stabilizes S S S x and yields a non-degenerate solution for our sparse factor model.

S.2.3 Selecting the number of factors
In order to select the number of latent factors r, we follow Onatski (2010). To the best of our knowledge, Onatski's method is the only one, which does not explicitly require that all factors are strong. Therefore, it is suitable for our setting, which allows as well for weak factors. The method uses the difference in subsequent eigenvalues and chooses the largest r such that: {r ≤ r max : πr((X X X ′ X X X)/T ) − πr +1 ((X X X ′ X X X)/T ) > ξ}, where ξ is a fixed positive constant, r max is an upper bound for the possible number of factors and πr((X X X ′ X X X)/T ) denotes ther-th largest eigenvalue of the covariance matrix of X X X. For the choice of ξ, the empirical distribution of the eigenvalues of the data sample covariance matrix is taken into account. j However, the estimation of the number of factors based on the empirical distribution of the eigenvalues of the sample covariance matrix still requires a clear separation of the eigenvalues of the common and idiosyncratic component.
Therefore, its selection accuracy depends on the degree of differentiability between the two components. Nevertheless, even if the selection method of Onatski (2010) overestimates the true number of factors, the sparsity assumption in our setting would allow us to disentangle the informative factors from those that are too weak. Thus, compared to the standard approximate factor model we avoid including redundant factors that amplify the misspecification error. Moreover, to further support the above argument, we refer to Yu and Samworth (2013), who show that in the weak factor setting the true number of factors is not asymptotically overestimated.

S.2.4 Choosing the tuning parameter
As for any penalized estimation approach, the selection of the tuning parameter µ is crucial, as it controls the degree of sparsity in the factor loadings matrix and it affects the performance of our estimator. In our case we select µ based on a type of Bayesian information j We refer to Onatski (2010) for the detailed description of the determination of ξ.
criterion, according to: where κ µ denotes the number of non-zero elements in the factor loadings matrixΛ Λ Λ for a given value of µ and L Λ Λ Λ, S S SF ,Σ Σ Σ τ u is the value of the log-likelihood function in equation (3), evaluated at the estimates of the factors, the factor loadings and the covariance matrix of the idiosyncratic errors. The penalty term in (S.35) has the property of converging to zero as both N and T approach infinity. Hence, the penalization vanishes as the sample size increases and a smaller value for µ is selected. The characteristics of our information criterion are therefore convenient with respect to the asymptotic properties we require for the regularization parameter µ. In fact, we need µ = o(1) in order to achieve estimation consistency, as elaborated in Section 3. The representation of the penalty term is based on the convergence rate of the factor loadings estimator in Lemma S.1.6.
To select the optimal µ, we estimate the criterion in (S.35) for a grid of different values for µ and choose the one that minimizes our information criterion. For the grid of the shrinkage parameter we consider the interval µ = (0, µ max ), where µ max denotes the highest value for the shrinkage parameter such that all imposed model restrictions are still fulfilled.

S.3 Competing Approaches
In this section, we summarize the estimation methods that are used in the simulation study and in the empirical application. k k We also included in our extended comparative study the approaches by Frahm and Memmel (2010) and Pollak (2011), which are based on direct shrinkage of the portfolio weights. The performance of these two models was clearly inferior, so that we refrained from giving the results here. However, they can be obtained from the authors upon request.
• Equally Weighted Portfolio (1/N) The equally weighted or 1/N portfolio strategy comprises identical portfolio weights of size 1/N, for each of the risky assets. By ignoring any type of optimizing portfolio strategy it often serves as a benchmark case to be outperformed in empirical performance comparisons.
As the weights have not to be estimated, the 1/N-strategy is free from any estimation risk. Moreover, the 1/N portfolio weights can be considered as the outcome for portfolio weights under extreme l 2 -penalization. DeMiguel, Garlappi, and Uppal (2009) find that the mean-variance portfolio and most of its extensions cannot significantly outperform the 1/N portfolio. l

• Sample (GMVP)
As the extreme alternative to the 1/N-strategy, we consider the plug-in estimator of the GMVP based on the sample covariance matrix of the asset returns. The plug-in estimator is free from any type of regularization. The plug-in approach yields unbiased estimates of the true weights (Okhrin and Schmid (2006)), but the weight estimates are extremely unstable when the asset space is large relative to the time series dimension. For some of our empirical designs with N = 100, 200, the asset dimension exceeds the sample size, T = 60.
For these cases the plug-in estimator is infeasible, because the sample covariance matrix is singular.

A. Factor Models
Factor models with latent factors l Kazak and Pohlmeier (2018) show, however, that conventional portfolio performance tests suffer from very low power, so that the rejection of null hypothesis of equal performance of a given data-based strategy and the 1/N -strategy is very unlikely.

Fan, Liao, and Mincheva (2013) (POET)
In our comparative study we include the POET estimator by Fan, Liao, and Mincheva (2013) that is based on the standard approximate factor model with a dense factor loadings matrix and a sparse idiosyncratic error covariance matrix. Similar to SAF, we use the number of factors selected by Onatski (2010).
2. Doz, Giannone, and Reichlin (2011) (DFM) To allow for some dynamics in the latent factors, we consider also a dynamic factor model originally proposed by Geweke (1977). Specifically, the dynamic factor model is represented by the following equation: where B B B i (L) = b i1 + b i2 L + · · · + b ip L p and L corresponds to the lag operator such that, In this setup f f f t = f 1t , f 2t , . . . , f qt ′ is a (q × 1)-dimensional vector of dynamic factors following a VAR process and b ij , j = 1, . . . , p denote the corresponding q-dimensional factor loadings. In order to estimate the dynamic factor model in (S.36), we use the two step procedure of Doz, Giannone, and Reichlin (2011). The estimation requires that the number of dynamic factors is given ex-ante. We use the consistent method by Bai and Ng (2007) to determine q.

Factor models with observable factors
In addition, we consider two factor models that have been frequently used in the empirical finance literature. Contrary to the approximate factor models, the factors in these models are not latent but observable time series variables. In this respect, these type of models incorporate more information than approaches, which solely use the information on the return process itself to estimate the covariance matrix of returns. However, the inclusion of additional time series information may give rise to an additional source of misspecification, if the factor specification fails to describe the true data generating process properly.

The Single Index Model (SIM)
The single index model by Sharpe (1963) is based on a single observable factor, f 1t , representing the excess market return: In our study, we use as a proxy for the market return, the value-weighted returns of all Center for Research in Security Prices (CRSP) firms incorporated in the US and listed on the AMEX, NASDAQ, or the NYSE. The one-month treasury bill rate serves as the risk free rate to construct the excess market returns. The estimator for the covariance matrix of the single index model is given by: whereσ f 1 denotes the sample variance of the market excess returns.β 1 represents the OLS estimates of the factor loadings andD D D is a diagonal matrix of the OLS residual variances of regression model (S.37) assuming that the observed factor picks up the cross-correlations of the returns completely.

Fama and French 3-Factor Model (FF3F)
The Fama and French 3-factor model offers an extension to the single index model by Sharpe (1963) and is defined as: x x x t = β β β 1 f 1t + β β β 2 f 2t + β β β 3 f 3t + ε ε ε t . (S.38) The first factor f 1 is identical to the one of the one-factor model in (S.37). The second factor f 2t often denoted by the acronym SMB is composed as the average returns on the three small portfolios minus the average returns on the three big portfolios. In particular, it defines a zero-cost portfolio that is long in stocks with a small market capitalization and short in stocks with a large market capitalization m . The third factor f 3t , denoted as HML, comprises a zero-cost portfolio that is long in stocks with a high book-to-market value and short in low book-to-market stocks n . In matrix notation (S.38) is given by: X X X = β β β F F F ′ + ε ε ε , (S.39) where F F F = [f f f 1 , f f f 2 , f f f 3 ] with dimension T × 3 and β β β = [β β β 1 , β β β 2 , β β β 3 ] with dimension N × 3.
The estimator for the covariance matrix for the 3-factor model by Fama and French (1993) Σ Σ Σ FF is equal to the following equation: whereΣ Σ Σ F denotes the covariance matrix of the three factors andD D D FF represents a diagonal matrix that contains the variances of the OLS residuals covariance matrix on its main diagonal.
m It is important to note that securities with a long position in a portfolio are expected to rise in value and on the other hand securities with short positions in a portfolio are expected to decline in value. n A detailed definition of the factors can be found on the website of Kenneth R. French. See http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data library.html 1. Ledoit and Wolf (2003) (LW) The LW approach shrinks the sample covariance matrix S S S x towards the covariance matrix of a single index model that is well-conditioned. This yields the following definition: where α * ∈ (0, 1) is a constant, which corresponds to the shrinkage intensity. Ledoit and Wolf (2003) propose the following estimator to be used in practiceα * = 1 T τ −ρ γ , where τ denotes the error on the sample covariance matrix, ρ measures the covariance between the estimation errors ofΣ Σ Σ SIM and S S S x and γ accounts for the misspecification of the shrinkage targetΣ Σ Σ SIM .

Kourtis, Dotsis, and Markellos (2012) (KDM)
The estimation method by Kourtis, Dotsis, and Markellos (2012) directly shrinks the inverse of the sample-based covariance matrix S S S x towards the identity matrix I N and the inverse of the covariance matrix resulting from a single index model by Sharpe (1963) The authors show that the resulting weights are a three-fund strategy, i.e. a linear combination of the sample-based weightsω ω ω, the equally weighted portfolio weightsω ω ω 1/N and those of the single index model modelω ω ω SIM . In order to select the optimal shrinkage coefficients in (S.40), the authors suggest minimizing the out-of-sample portfolio variance using cross-validation. It is important to note that this portfolio strategy is also applicable for the case when N > T . In order to obtain reliable results for the inverse of S S S x the authors use the Moore-Penrose pseudo-inverse.

Abadir, Distaso, andŽikeš (2014) (ADZ)
The design-free estimator for the covariance matrix by Abadir, Distaso, andŽikeš (2014) aims to improve the estimation of the eigenvaluesP P P of S S S x , that is a possible source of ill-conditioning. The authors consider the following spectral decomposition of S S S x : In order to obtain an improved estimator for P P P , X X X is split into two subsamples X X X = X X X 1 N ×n , X X X 2 N ×(T −n) . Calculating the sample covariance matrix for the first n observations yields: where M M M n = I I I n − 1 n 1 n 1 ′ n is the de-meaning matrix of dimension n and 1 n denotes a n × 1 vector of ones. The spectral decomposition of S S S 1 provides the matrix of eigenvectorsΓ Γ Γ 1 and the diagonal matrix of eigenvaluesP P P 1 .
In the second step, an improved estimator for P P P is computed from the remaining or- Another method that aims to improve on the estimation of the eigenvalues of S S S x is provided by Ledoit and Wolf (2018). The covariance estimator is given by: (S.45) whereΓ Γ Γ are the sample eigenvectors of S S S x and the eigenvalues in the diagonal matrixD D D are estimated in a non-linear fashion as in Theorem 6.2. in Ledoit and Wolf (2018).

C. Sparse Covariance Estimators
The following estimators are explicitly designed to provide sparse covariance matrices.
Hence, these models are appropriate for empirical settings that are reflected by our second simulation design.

Bien and Tibshirani (2011) (BT)
The authors propose a penalized maximum likelihood estimator based on a lasso penalty in order to allow for sparsity in the covariance matrix and to reduce the effective number of parameters. More specifically, the following objective function is optimized: where α N is a regularization parameter selected based on 5-fold cross-validation. The ij-th element of the selection matrix H is defined as h ij = 1l {i = j} and enables an equal penalization of the off-diagonal elements and leaves the diagonal elements unaffected.
Furthermore, Bien and Tibshirani (2011) show that the estimated sparse covariance matrix is positive definite.