Optimal Portfolio Using Factor Graphical Lasso ⁰⁰footnotetext: The authors would like to thank the editor Fabio Trojani and three anonymous referees for helpful and constructive comments on the paper.

Estimating the inverse covariance matrix, or precision matrix, of excess stock returns is crucial for constructing weights of financial assets in a portfolio and estimating the out-of-sample Sharpe Ratio. In high-dimensional setting, when the number of assets, $p$, is greater than or equal to the sample size, $T$, using an estimator of covariance matrix for obtaining portfolio weights leads to unstable investment allocations. This is known as the Markowitz’ curse: a higher number of assets increases correlation between the investments, which calls for a more diversified portfolio, and yet unstable corner solutions for weights become more likely. The reason behind this curse is the need to invert a high-dimensional covariance matrix to obtain the optimal weights from the quadratic optimization problem: when $p\geq T$, the condition number of the covariance matrix (i.e., the absolute value of the ratio between maximal and minimal eigenvalues of the covariance matrix) is high. Hence, the inverted covariance matrix yields an unstable estimator of the precision matrix. To circumvent this issue one can estimate precision matrix directly, rather than inverting an estimated covariance matrix.

Graphical models were shown to provide consistent estimates of the precision matrix (Friedman et al., (2008); Meinshausen and Bühlmann, (2006); Cai et al., (2011)). Goto and Xu, (2015) estimated a sparse precision matrix for portfolio hedging using graphical models. They found out that their portfolio achieves significant out-of-sample risk reduction and higher return, as compared to the portfolios based on equal weights, shrunk covariance matrix, industry factor models, and no-short-sale constraints. Awoye, (2016) used Graphical Lasso (Friedman et al., (2008)) to estimate a sparse covariance matrix for the Markowitz mean-variance portfolio problem and reduce the realized portfolio risk. Millington and Niranjan, (2017) conducted an empirical study that applies Graphical Lasso for the estimation of covariance for the portfolio allocation. Their empirical findings suggest that portfolios using Graphical Lasso enjoy lower risk and higher returns compared to those using empirical covariance matrix. Millington and Niranjan, (2017) also construct a financial network using the estimated precision matrix to explore the relationship between the companies and show how the constructed network helps to make investment decisions. Callot et al., (2021) use the nodewise-regression method of Meinshausen and Bühlmann, (2006) to establish consistency of the estimated covariance matrix, weights and risk of high-dimensional financial portfolio. Their empirical application demonstrates that the precision matrix estimator based on the nodewise-regression outperforms the principal orthogonal complement thresholding estimator (POET) (Fan et al., (2013)) and linear shrinkage (Ledoit and Wolf, (2004)). Cai et al., (2020) use constrained $\ell_{1}$-minimization for inverse matrix estimation (Clime) of the precision matrix (Cai et al., (2011)) to develop a consistent estimator of the minimum variance for high-dimensional global minimum-variance portfolio. It is important to note that all the aforementioned methods impose some sparsity assumption on the precision matrix of excess returns.

An alternative strategy to handle high-dimensional setting uses factor models to acknowledge common variation in the stock prices, which was documented in many empirical studies (see Campbell et al., (1997) among many others). A common approach decomposes covariance matrix of excess returns into low-rank and sparse parts, the latter is further regularized since, after the common factors are accounted for, the remaining covariance matrix of the idiosyncratic components is still high-dimensional (Fan et al., (2013, 2011, 2018)). This stream of literature, however, focuses on the estimation of a covariance matrix. The accuracy of precision matrices obtained from inverting the factor-based covariance matrix was investigated by Ait-Sahalia and Xiu, (2017), but they did not study a high-dimensional case. Factor models are generally treated as competitors to graphical models: as an example, Callot et al., (2021) find evidence of superior performance of nodewise-regression estimator of precision matrix over a factor-based estimator POET (Fan et al., (2013)) in terms of the out-of-sample Sharpe Ratio and risk of financial portfolio. The root cause why factor models and graphical models are treated separately is the sparsity assumption on the precision matrix made in the latter. Specifically, as pointed out in Koike, (2020), when asset returns have common factors, the precision matrix cannot be sparse because all pairs of assets are partially correlated conditional on other assets through the common factors. One attempt to integrate factor modeling and high-dimensional precision estimation was made by Fan et al., (2018) (Section 5.2): the authors referred to such class of models as “conditional graphical models”. However, this was not the main focus of their paper which concentrated on covariance estimation through elliptical factor models. As Fan et al., (2018) pointed out, “though substantial amount of efforts have been made to understand the graphical model, little has been done for estimating conditional graphical model, which is more general and realistic”. Concretely, to the best of our knowledge there are no studies that examine theoretical and empirical performance of graphical models integrated with the factor structure in the context of portfolio allocation.

In this paper we fill this gap and develop a new conditional precision matrix estimator for the excess returns under the approximate factor model that combines the benefits of graphical models and factor structure. We call our algorithm the Factor Graphical Lasso (FGL). We use a factor model to remove the co-movements induced by the factors, and then we apply the Weighted Graphical Lasso for the estimation of the precision matrix of the idiosyncratic terms. We prove consistency of FGL in the spectral and $\ell_{1}$ matrix norms. In addition, we prove consistency of the estimated portfolio weights and risk exposure for three formulations of the optimal portfolio allocation.

Our empirical application uses daily and monthly data for the constituents of the S&P500: we demonstrate that FGL outperforms equal-weighted portfolio, index portfolio, portfolios based on other estimators of precision matrix (Clime, Cai et al., (2011)) and covariance matrix, including POET (Fan et al., (2013)) and the shrinkage estimators adjusted to allow for the factor structure (Ledoit and Wolf, (2004), Ledoit and Wolf, (2017)), in terms of the out-of-sample Sharpe Ratio. Furthermore, we find strong empirical evidence that relaxing the constraint that portfolio weights sum up to one leads to a large increase in the out-of-sample Sharpe Ratio, which, to the best of our knowledge, has not been previously well-studied in the empirical finance literature.

From the theoretical perspective, our paper makes several important contributions to the existing literature on graphical models and factor models. First, to the best of out knowledge, there are no equivalent theoretical results that establish consistency of the portfolio weights and risk exposure in a high-dimensional setting without assuming sparsity on the covariance or precision matrix of stock returns. Second, we extend the theoretical results of POET (Fan et al., (2013)) to allow the number of factors to grow with the number of assets. Concretely, we establish uniform consistency for the factors and factor loadings estimated using PCA. Third, we are not aware of any other papers that provide convergence results for estimating a high-dimensional precision matrix using the Weighted Graphical Lasso under the approximate factor model with unobserved factors. Furthermore, all theoretical results established in this paper hold for a wide range of distributions: Sub-Gaussian family (including Gaussian) and elliptical family. Our simulations demonstrate that FGL is robust to very heavy-tailed distributions, which makes our method suitable for the financial applications. Finally, we demonstrate that in contrast to POET, the success of the proposed method does not heavily depend on the factor pervasiveness assumption: FGL is robust to the scenarios when the gap between the diverging and bounded eigenvalues decreases.

This paper is organized as follows: Section 2 reviews the basics of the Markowitz mean-variance portfolio theory. Section 3 provides a brief summary of the graphical models and introduces the Factor Graphical Lasso. Section 4 contains theoretical results and Section 5 validates these results using simulations. Section 6 provides empirical application. Section 7 concludes.

Suppose we observe $p$ assets (indexed by $i$) over $T$ period of time (indexed by $t$). Let $\widetilde{{\mathbf{r}}}_{t}=(\widetilde{r}_{1t},\widetilde{r}_{2t},\ldots,\widetilde{r}_{pt})^{\prime}\sim\mathcal{D}({\mathbf{m}},{\bm{\Sigma}})$ be a $p\times 1$ vector of excess returns drawn from a distribution $\mathcal{D}$, where ${\mathbf{m}}$ and ${\bm{\Sigma}}$ are the unconditional mean and covariance matrix of the returns. The goal of the Markowitz theory is to choose asset weights in a portfolio optimally. We will study two optimization problems: the well-known Markowitz weight-constrained (MWC) optimization problem, and the Markowitz risk-constrained (MRC) optimization that relaxes the constraint on portfolio weights.

The first optimization problem searches for asset weights such that the portfolio achieves a desired expected rate of return with minimum risk, under the restriction that all weights sum up to one. This can be formulated as the following quadratic optimization problem:

where ${\mathbf{w}}$ is a $p\times 1$ vector of asset weights in the portfolio, ${\bm{\iota}}_{p}$ is a $p\times 1$ vector of ones, and $\mu$ is a desired expected rate of portfolio return. Let ${\bm{\Theta}}\equiv{\bm{\Sigma}}^{-1}$ be the precision matrix.

If ${\mathbf{m}}^{\prime}{\mathbf{w}}>\mu$, then the solution to (2.1) yields the global minimum-variance (GMV) portfolio weights ${\mathbf{w}}_{GMV}$:

The MRC problem maximizes Sharpe Ratio (SR) subject to either target risk or target return constraints, but portfolio weights are not required to sum up to one:

When $\mu=\sigma\sqrt{{\mathbf{m}}^{\prime}{\bm{\Theta}}{\mathbf{m}}}$, the solution to either of the constraints is given by

Equation (2.4) tells us that once an investor specifies the desired return, $\mu$, and maximum risk-tolerance level, $\sigma$, the MRC weight maximizes the Sharpe Ratio of the portfolio.

Therefore, we have three alternative portfolio allocations commonly used in the existing literature: GMV in (2.2), MWC in (2.3) and MRC in (2.5). It is clear that all formulations require an estimate of the precision matrix ${\bm{\Theta}}$.

In this section we introduce a framework for estimating precision matrix for the aforementioned financial portfolios which accounts for the fact that the returns follow approximate factor structure. We examine how to solve the Markowitz mean-variance portfolio allocation problems using factor structure in the returns. We also develop Factor Graphical Lasso Algorithm that uses the estimated common factors to obtain a sparse precision matrix of the idiosyncratic component. The resulting estimator is used to obtain the precision of the asset returns necessary to form portfolio weights.

The arbitrage pricing theory (APT), developed by Ross, (1976), postulates that the expected returns on securities should be related to their covariance with the common components or factors. The goal of the APT is to model the tendency of asset returns to move together via factor decomposition. Assume that the return generating process ($\widetilde{{\mathbf{r}}}_{t}$) follows a $K$-factor model:

where ${\mathbf{f}}_{t}=(f_{1t},\ldots,f_{Kt})^{\prime}$ are the factors, ${\mathbf{B}}$ is a $p\times K$ matrix of factor loadings, and ${\bm{\varepsilon}}_{t}$ is the idiosyncratic component that cannot be explained by the common factors. Without loss of generality, we assume throughout the paper that unconditional means of factors and idiosyncratic component are zero. Factors in (3.1) can be either observable, such as in Fama and French, (1993, 2015), or can be estimated using statistical factor models. Unobservable factors and loadings are usually estimated by the principal component analysis (PCA), as studied in Bai, (2003), Bai and Ng, (2002), Connor and Korajczyk, (1988), and Stock and Watson, (2002).

In this paper our main interest lies in establishing asymptotic properties of the estimators of precision matrix, portfolio weights and risk-exposure for the high-dimensional case. We assume that the number of common factors, $K=K_{p,T}\rightarrow\infty$ as $p\rightarrow\infty$, or $T\rightarrow\infty$, or both $p,T\rightarrow\infty$, but we require that $\max\{K/p,K/T\}\rightarrow 0$ as $p,T\rightarrow\infty$.

Our setup is similar to the one studied in Fan et al., (2013): we consider a spiked covariance model when the first $K$ principal eigenvalues of ${\bm{\Sigma}}$ are growing with $p$, while the remaining $p-K$ eigenvalues are bounded.

Rewrite equation (3.1) in matrix form:

where ${\bm{\iota}}_{T}$ is a $T\times 1$ vector of ones. We further demean the returns using the sample mean, $\widehat{{\mathbf{m}}}$, to obtain ${\mathbf{R}}\equiv\widetilde{{\mathbf{R}}}-\widehat{{\mathbf{m}}}{\bm{\iota}}^{\prime}_{T}$. We assume that $\left\lVert\widehat{{\mathbf{m}}}-{\mathbf{m}}\right\rVert_{\text{max}}=\mathcal{O}_{P}(\sqrt{\log p/T})$, which was proven to hold in Chang et al., (2018) (see their Lemma 1).

Let ${\bm{\Sigma}}_{\varepsilon}=T^{-1}{\mathbf{E}}{\mathbf{E}}^{\prime}$ and ${\bm{\Sigma}}_{f}=T^{-1}{\mathbf{F}}{\mathbf{F}}^{\prime}$ be covariance matrices of the idiosyncratic components and factors, and let ${\bm{\Theta}}_{\varepsilon}={\bm{\Sigma}}_{\varepsilon}^{-1}$ and ${\bm{\Theta}}_{f}={\bm{\Sigma}}_{f}^{-1}$ be their inverses. The factors and loadings in (3.2) are estimated by solving the following minimization problem: $(\widehat{{\mathbf{B}}},\widehat{{\mathbf{F}}})=\arg\!\min_{{\mathbf{B}},{\mathbf{F}}}\left\lVert{\mathbf{R}}-{\mathbf{B}}{\mathbf{F}}\right\rVert^{2}_{F}$ s.t. $\frac{1}{T}{\mathbf{F}}{\mathbf{F}}^{\prime}={\mathbf{I}}_{K},\ {\mathbf{B}}^{\prime}{\mathbf{B}}\ \text{is diagonal}$. The constraints are needed to identify the factors (Fan et al., (2018)). It was shown (Stock and Watson, (2002)) that $\widehat{{\mathbf{F}}}=\sqrt{T}\text{eig}_{K}({\mathbf{R}}^{\prime}{\mathbf{R}})$ and $\widehat{{\mathbf{B}}}=T^{-1}{\mathbf{R}}\widehat{{\mathbf{F}}}^{\prime}$. Given $\widehat{{\mathbf{F}}},\widehat{{\mathbf{B}}}$, define $\widehat{{\mathbf{E}}}={\mathbf{R}}-\widehat{{\mathbf{B}}}\widehat{{\mathbf{F}}}$. Given a sample of the estimated residuals $\{\widehat{{\bm{\varepsilon}}}_{t}={\mathbf{r}}_{t}-\widehat{{\mathbf{B}}}\widehat{{\mathbf{f}}_{t}}\}_{t=1}^{T}$ and the estimated factors $\{\widehat{{\mathbf{f}}}_{t}\}_{t=1}^{T}$, let $\widehat{{\bm{\Sigma}}}_{\varepsilon}=(1/T)\sum_{t=1}^{T}\widehat{{\bm{\varepsilon}}}_{t}\widehat{{\bm{\varepsilon}}}_{t}^{\prime}$ and $\widehat{{\bm{\Sigma}}}_{f}=(1/T)\sum_{t=1}^{T}\widehat{{\mathbf{f}}}_{t}\widehat{{\mathbf{f}}}_{t}^{\prime}$ be the sample counterparts of the covariance matrices. Since our interest is in constructing portfolio weights, our goal is to estimate a precision matrix of the excess returns ${\bm{\Theta}}$.

We impose a sparsity assumption on the precision matrix of the idiosyncratic errors, ${\bm{\Theta}}_{\varepsilon}$, which is obtained using the estimated residuals after removing the co-movements induced by the factors (see Barigozzi et al., (2018); Brownlees et al., (2018); Koike, (2020)).

Let us elaborate on three reasons justifying the assumption of sparsity on the precision matrix of residuals. First, from the technical viewpoint, this assumption is widely used in high-dimensional settings when $p>T$. Second, a more intuitive rationale for the sparsity assumption on ${\bm{\Theta}}_{\varepsilon}$ stems from its implication for the structure of corresponding optimal portfolios. Let $r_{t}^{\text{portf}}\equiv\widetilde{{\mathbf{r}}}_{t}^{\prime}{\mathbf{w}}_{t}$ be the optimal portfolio. Plugging in the definition of $\widetilde{{\mathbf{r}}}_{t}$ from (3.1), we get $r_{t}^{\text{portf}}=({\mathbf{m}}+{\bm{\varepsilon}}_{t})^{\prime}{\mathbf{w}}_{t}+{\mathbf{f}}_{t}^{\prime}{\mathbf{B}}{\mathbf{w}}_{t}$. Hence, after hedging factor risk, we can isolate the excess return component only loading on non-factor risk. In this context, since ${\mathbf{w}}_{t}$ is a function of ${\bm{\Theta}}_{\varepsilon}$, imposing sparsity on ${\bm{\Theta}}_{\varepsilon}$ translates into reducing the contribution of more volatile non-factor risk on the optimal portfolio and thus leading to less sensitive (more robust) investment strategies.

Third, another rationale comes from relatively high “concentration” of S&P 500 Composite Index: as evidenced from SP Global Index methodology and financial data on S&P 500 constituents by weight, 15 large companies (top 3%) comprise 30% of the total index weights (starting from Apple that has the highest weight of nearly 7%). As the number of firms, $p$, increases, one reasonable assumption is that the number of large firms increases at a rate slower than $p$ (Chudik et al., (2011); Gabaix, (2011)). This suggests that one could divide the firms into dominant ones and followers. After the effect of common factors is accounted for, dominant firms still have significant idiosyncratic movements that influence other firms and must be taken into account when constructing a portfolio. When it comes to fringe firms (or market followers), idiosyncratic movements are smaller in magnitude and might be less relevant for portfolio allocation purposes. Hence, the network of the idiosyncratic returns is sparse and the sparsity increases with $p$. By imposing sparsity, we only keep relatively large partial correlations among idiosyncratic components: as illustrated in Supplemental Appendix D.2, in our empirical application the estimated number of zeroes in off-diagonal elements of ${\bm{\Theta}}_{\varepsilon}$ varies over time from 74.5%-98.8%.

A common way to induce sparsity is by utilizing Lasso-type penalty. This strategy is used in the Graphical Lasso (GL) together with the objective function based on the Bregman divergence for estimating inverse covariance. The discussion of GL is presented in Supplemental Appendix A. We now elaborate on the Bregman divergence class which unifies many commonly used loss functions, including the quasi-likelihood function. Let ${\mathbf{W}}_{\varepsilon}$ be an estimate of ${\bm{\Sigma}}_{\varepsilon}$. Ravikumar et al., (2011) showed that Bregman divergence of the form $\text{trace}({\mathbf{W}}_{\varepsilon}{\bm{\Theta}}_{\varepsilon})-\log\det({\bm{\Theta}}_{\varepsilon})$, known as the log-determinant Bregman function, is suitable to be used as a measure of the quality of constructed sparse approximations of signals such as precision matrices. As pointed out by Ravikumar et al., (2011), in principle one could use other Bregman divergences including the von Neumann Entropy or the Frobenius divergence which would lead to alternative forms of divergence minimizations for estimating precision matrix. We proceed with the log-determinant Bregman function since (i) it ensures positive definite estimator of precision matrix; (ii) the population optimization problem involves only the population covariance and not its inverse; (iii) the log-determinant divergence gives rise to the likelihood function in the multivariate Gaussian case. At the same time, despite its resemblance with the Gaussian log-likelihood, Bregman divergence was shown to be applicable for non-Gaussian distributions (Ravikumar et al., (2011)). Let $\widehat{{\mathbf{D}}}_{\varepsilon}^{2}\equiv\textup{diag}({\mathbf{W}}_{\varepsilon})$. To sparsify entries of precision matrix of the idiosyncratic errors ${\bm{\Theta}}_{\varepsilon}$, we use the following penalized Bregman divergence with the Weighted Graphical Lasso penalty:

The subscript $\lambda$ in $\widehat{{\bm{\Theta}}}_{\varepsilon,\lambda}$ means that the solution of the optimization problem in (3.3) will depend upon the choice of the tuning parameter which is discussed below. Section 4 establishes sparsity requirements that guarantee convergence of (3.3). In order to simplify notation, we will omit the subscript $\lambda$.

The objective function in (3.3) extends the family of linear shrinkage estimators of the first moment to linear shrinkage estimators of the inverse of the second moments. Instead of restricting the number of regressors for estimating conditional mean, equation (3.3) restricts the number of edges in a graph by shrinking some off-diagonal entries of precision matrix to zero. Note that shrinkage occurs adaptively with respect to partial covariances.

Let us discuss the choice of the tuning parameter $\lambda$ in (3.3). Let $\widehat{{\bm{\Theta}}}_{\varepsilon,\lambda}$ be the solution to (3.3) for a fixed $\lambda$. Following Koike, (2020), we minimize the following Bayesian Information Criterion (BIC) using grid search:

The grid $\mathcal{G}\equiv\{\lambda_{1},\ldots,\lambda_{M}\}$ is constructed as follows: the maximum value in the grid, $\lambda_{M}$, is set to be the smallest value for which all the off-diagonal entries of $\widehat{{\bm{\Theta}}}_{\varepsilon,\lambda_{M}}$ are zero. The smallest value of the grid, $\lambda_{1}\in\mathcal{G}$, is determined as $\lambda_{1}\equiv\vartheta\lambda_{M}$ for a constant $0<\vartheta<1$. The remaining grid values $\lambda_{1},\ldots,\lambda_{M}$ are constructed in the ascending order from $\lambda_{1}$ to $\lambda_{M}$ on the log scale:

We use $\vartheta=\omega_{3T}$ which is defined in Theorem 2 of the next section and $M=10$ in the simulations and the empirical exercise.

Having estimated factors, factor loadings and precision matrix of the idiosyncratic components, we combine them using the Sherman-Morrison-Woodbury formula to estimate the precision matrix of excess returns:

To solve (3.3) we use the procedure based on the GL. However, the original algorithm developed by Friedman et al., (2008) is not suitable under the factor structure. Our procedure called Factor Graphical Lasso (FGL), which is summarized in Procedure 1, augments the standard GL: it starts with estimating factors, loadings (low-rank part) and error terms (sparse part), then it proceeds by recovering sparse precision matrix of the errors using GL, and, finally, low-rank and sparse components are combined through Shermann-Morrison-Woodbury formula in (3.5).

The estimator produced by GL in general and FGL in particular is guaranteed to be positive definite. We have verified it in the simulations (Section 5) and the empirical application (Section 6). In Section 4, consistency properties of estimators are established for the factors and loadings (Theorem 1), the precision matrix of ${\bm{\varepsilon}}$ (Theorem 2), the precision matrix ${\bm{\Theta}}$ (Theorem 3), portfolio weights (Theorem 4), and the portfolio risk exposure (Theorem 5). We can use $\widehat{{\bm{\Theta}}}$ obtained from (3.5) using Step 4 of Procedure 1 to estimate portfolio weights in (2.2), (2.3) and (2.5):

In this section we first provide a brief review of the terminology used in the literature on graphical models and the approaches to estimate a precision matrix. After that we establish consistency of the Factor Graphical Lasso in Procedure 1. We also study consistency of the estimators of weights in (2.2), (2.3) and (2.5) and the implications on the out-of sample Sharpe Ratio. Throughout the main text we assume that errors and factors have exponential-type tails (1(c)). Supplemental Appendix B.10 proves that the conclusions of all theorems studied in Section 4 continue to hold when this assumption is relaxed.

The review of the Gaussian graphical models is based on Hastie et al., (2001) and Bishop, (2006). A graph consists of a set of vertices (nodes) and a set of edges (arcs) that join some pairs of the vertices. In graphical models, each vertex represents a random variable, and the graph visualizes the joint distribution of the entire set of random variables. The edges in a graph are parameterized by potentials (values) that encode the strength of the conditional dependence between the random variables at the corresponding vertices. Sparse graphs have a relatively small number of edges. Among the main challenges in working with the graphical models are choosing the structure of the graph (model selection) and estimation of the edge parameters from the data.

Let $A\in\mathcal{S}_{p}$. Define the following set for $j=1,\ldots,p$:

where $d_{j}(A)$ is the number of edges adjacent to the vertex $j$ (i.e., the degree of vertex $j$), and $d(A)$ measures the maximum vertex degree. Define $S(A)\equiv\bigcup_{j=1}^{p}D_{j}(A)$ to be the overall off-diagonal sparsity pattern, and $s(A)\equiv\sum_{j=1}^{p}d_{j}(A)$ is the overall number of edges contained in the graph. Note that $\text{card}(S(A))\leq s(A)$: when $s(A)=p(p-1)/2$ this would give a fully connected graph.

In order to validate our theoretical results, we perform several simulation studies which are divided into four parts. The first set of results computes the empirical convergence rates and compares them with the theoretical expressions derived in Theorems 3-5. The second set of results compares the performance of the FGL with several alternative models for estimating covariance and precision matrix. To highlight the benefit of using the information about factor structure as opposed to standard graphical models, we include Graphical Lasso by Friedman et al., (2008) (GL) that does not account for the factor structure. To explore the benefits of using FGL for error quantification in (3.5), we consider several alternative estimators of covariance/precision matrix of the idiosyncratic component in (3.5): (1) linear shrinkage estimator of covariance developed by Ledoit and Wolf, (2004) further referred to as Factor LW or FLW; (2) nonlinear shrinkage estimator of covariance by Ledoit and Wolf, (2017) (Factor NLW or FNLW); (3) POET (Fan et al., (2013)); (4) constrained $\ell_{1}$-minimization for inverse matrix estimator, Clime (Cai et al., (2011)) (Factor Clime or FClime). Furthermore, we discovered that in certain setups the estimator of covariance produced by POET is not positive definite. In such cases we use the matrix symmetrization procedure as in Fan et al., (2018) and then use eigenvalue cleaning as in Callot et al., (2017) and Hautsch et al., (2012). This estimator is referred to as Projected POET; it coincides with POET when the covariance estimator produced by the latter is positive definite. The third set of results examines the performance of FGL and Robust FGL (described in Supplemental Appendix B.10) when the dependent variable follows elliptical distribution. The fourth set of results explores the sensitivity of portfolios constructed using different covariance and precision estimators of interest when the pervasiveness assumption 1 is relaxed, that is, when the gap between the diverging and bounded eigenvalues decreases. All exercises in this section use 100 Monte Carlo simulations.

We consider the following setup: let $p=T^{\delta}$, $\delta=0.85$, $K=2(\log T)^{0.5}$ and $T=[2^{h}],\ \text{for}\ h=7,7.5,8,\ldots,9.5$. A sparse precision matrix of the idiosyncratic components is constructed as follows: we first generate the adjacency matrix using a random graph structure. Define a $p\times p$ adjacency matrix ${\mathbf{A}}_{\varepsilon}$ which is used to represent the structure of the graph:

Let $a_{\varepsilon,ij}$ denote the $i,j$-th element of the adjacency matrix ${\mathbf{A}}_{\varepsilon}$. We set $a_{\varepsilon,ij}=a_{\varepsilon,ji}=1,\ \text{for}\ i\neq j$ with probability $q$, and $0$ otherwise. Such structure results in $s_{T}=p(p-1)q/2$ edges in the graph. To control sparsity, we set $q=1/(pT^{0.8})$, which makes $s_{T}=\mathcal{O}(T^{0.05})$. The adjacency matrix has all diagonal elements equal to zero. Hence, to obtain a positive definite precision matrix we apply the procedure described in Zhao et al., (2012): using their notation, ${\bm{\Theta}}_{\varepsilon}={\mathbf{A}}_{\varepsilon}\cdot v+{\mathbf{I}}(\left\lvert\tau\right\rvert+0.1+u)$, where $u>0$ is a positive number added to the diagonal of the precision matrix to control the magnitude of partial correlations, $v$ controls the magnitude of partial correlations with $u$, and $\tau$ is the smallest eigenvalue of ${\mathbf{A}}_{\varepsilon}\cdot v$. In our simulations we use $u=0.1$ and $v=0.3$.

Factors are assumed to have the following structure:

where $m_{i}\sim\mathcal{N}(1,1)$ independently for each $i=1,\ldots,p$, ${\bm{\varepsilon}}_{t}$ is a $p\times 1$ random vector of idiosyncratic errors following $\mathcal{N}(\bm{0},{\bm{\Sigma}}_{\varepsilon})$, with sparse ${\bm{\Theta}}_{\varepsilon}$ that has a random graph structure described above, ${\mathbf{f}}_{t}$ is a $K\times 1$ vector of factors, $\phi_{f}$ is an autoregressive parameter in the factors which is a scalar for simplicity, ${\mathbf{B}}$ is a $p\times K$ matrix of factor loadings, ${\bm{\zeta}}_{t}$ is a $K\times 1$ random vector with each component independently following $\mathcal{N}(0,\sigma^{2}_{\zeta})$. To create ${\mathbf{B}}$ in (5.3) we take the first $K$ rows of an upper triangular matrix from a Cholesky decomposition of the $p\times p$ Toeplitz matrix parameterized by $\rho$. For the first set of results we set $\rho=0.2$, $\phi_{f}=0.2$ and $\sigma^{2}_{\zeta}=1$. The specification in (5.3) leads to the low-rank plus sparse decomposition of the covariance matrix of stock returns ${\mathbf{r}}_{t}$.

As a first exercise, we compare the empirical and theoretical convergence rates of the precision matrix, portfolio weights and exposure. A detailed description of the procedure and the simulation results is provided in Supplemental Appendix C.1. We confirm that the empirical rates and theoretical rates from Theorems 3-5 are matched.

As a second exercise, we compare the performance of FGL with the alternative models listed at the beginning of this section. We consider two cases: Case 1 is the same as for the first set of simulations ($p<T$): $p=T^{\delta}$, $\delta=0.85$, $K=2(\log T)^{0.5}$, $s_{T}=\mathcal{O}(T^{0.05})$. Case 2 captures the cases when $p>T$ with $p=3\cdot T^{\delta}$, $\delta=0.85$, all else equal. The results for Case 2 are reported in Figure 1-3, and Case 1 is located in Supplemental Appendix C.2. FGL demonstrates superior performance for estimating precision matrix and portfolio weights in both cases, exhibiting consistency for both Case 1 and Case 2 settings. Also, FGL outperforms GL for estimating portfolio exposure and consistently estimates the latter, however, depending on the case under consideration some alternative models produce lower averaged error.

As a third exercise, we examine the performance of FGL and Robust FGL (described in Supplemental Appendix B.10) when the dependent variable follows elliptical distributions. A detailed description of the data generating process (DGP) and simulation results are provided in Supplemental Appendix C.3. We find that the performance of FGL for estimating the precision matrix is comparable with that of Robust FGL: this suggests that our FGL algorithm is robust to heavy-tailed distributions even without additional modifications.

As a final exercise, we explore the sensitivity of portfolios constructed using different covariance and precision estimators of interest when the pervasiveness assumption 1 is relaxed. A detailed description of the data generating process (DGP) and simulation results are provided in Supplemental Appendix C.4. We verify that FGL exhibits robust performance when the gap between the diverging and bounded eigenvalues decreases. In contrast, POET and Projected POET are most sensitive to relaxing pervasiveness assumption which is consistent with our empirical findings and also with the simulation results by Onatski, (2013).

In this section we examine the performance of the Factor Graphical Lasso for constructing a financial portfolio using daily data. The description and empirical results for monthly data can be found in Supplemental Appendix D. We first describe the data and the estimation methodology, then we list four metrics commonly reported in the finance literature, and, finally, we present the results.