Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
107,438 characters · 32 sections · 69 citation commands
A Simple Method for Predicting Covariance Matrices of Financial Returns
\makeabstracttitle
\chapter{Introduction}
We consider cross-sections, {\it e.g.}, a vector time series of $n$ financial returns, denoted $r_t\in {\mbox{\bf R}}^n$, $t=1,2,\ldots$, where $(r_t)_i$ is the return of asset $i$ from $t-1$ to $t$. We focus on the case where the mean $\mathop{\bf E{}} r_t$ is small enough that the second moment $\mathop{\bf E{}} r_tr_t^T\in {\mbox{\bf R}}^{n \times n}$ is a good approximation of the covariance matrix $\mathop{\bf cov{}}(r_t) = \mathop{\bf E{}} r_tr_t^T - (\mathop{\bf E{}} r_t)(\mathop{\bf E{}} r_t)^T$, where $\mathop{\bf E{}}$ denotes expectation. This is the case for most daily, weekly, or monthly stock, bond, and futures returns, factor returns, and index returns. We start by focussing on the case where the number of assets $n$ is modest, say, on the order 10--100 or so; in chapter (ref) we explain how to extend the method to much larger universes using ideas such as factor models.
We model the returns $r_t$ as independent random variables with zero mean and covariance $\Sigma_t \in {\mbox{\bf S}}_{++}^n$ (the set of symmetric positive definite matrices). We focus on the problem of predicting or estimating $\Sigma_t$, based on knowledge of $r_1, \ldots, r_{t-1}$. The prediction is denoted as $\hat{\Sigma}_t \in {\mbox{\bf S}}_{++}^n$. The predicted volatilities of assets are given by \[ \hat\sigma_t = \mathop{\bf diag} (\hat{\Sigma}_t)^{1/2}\in {\mbox{\bf R}}^n, \] where $\mathop{\bf diag}$ with a matrix argument is the vector of diagonal entries of the matrix, and the squareroot of a vector above is elementwise. We denote the predicted correlations as \[ \hat R_t = \mathop{\bf diag}(\hat\sigma_t)^{-1} \hat{\Sigma}_t \mathop{\bf diag}(\hat\sigma_t)^{-1}, \] where $\mathop{\bf diag}$ with a vector argument is the diagonal matrix with entries from the vector argument.
Covariance estimation comes up in several areas of finance, including Markowitz portfolio construction markowitz_1952,grinold2000_portfolio, risk management mcneil2015quantitative, and asset pricing sharpe_1964. Much attention has been devoted to this problem, and a Nobel Memorial Prize in Economic Sciences was awarded for work directly related to volatility estimation engle_1982.
While it is well known that the tails of financial returns are poorly modeled by a Gaussian distribution, our focus here is on the bulk of the distribution, where the Gaussian assumption is reasonable. For future use, we note that the log-likelihood of an observed return $r_t$, under the Gaussian distribution $r_t \sim \mathcal N(0,\hat{\Sigma}_t)$, is
The Gaussian log-likelihood is closely related to a popular metric for evaluating covariance predictors in econometrics, called the (Gaussian) quasi-likelihood (QLIKE) patton2011volatility, patton2009evaluating, laurent2013loss. QLIKE is the negative log-likelihood, under the Gaussian assumption, up to an additive constant and a positive scale factor. Roughly speaking, we seek covariance predictors that achieve large values of log-likelihood, or small values of QLIKE, on realized returns. We will describe evaluation of covariance predictors in detail in chapter (ref).
This monograph makes three contributions. First, we propose a new method for predicting the time-varying covariance matrix of a vector of financial returns, building on a specific covariance estimator suggested by Engle in 2002. Our method is a relatively simple extension that requires very little tuning and is readily interpretable. It relies on solving a small convex optimization problem, which can be carried out very quickly and reliably boyd2004convex. Our method performs as well as much more complex methods, as measured by several metrics.
Our second contribution is to propose a new method for evaluating a covariance predictor, by considering the regret of the log-likelihood over some time period such as a quarter. This approach allows us to evaluate how quickly a covariance estimator reacts to changes in market conditions.
Our third contribution is an extensive empirical study of covariance predictors. We compare our new method to other popular predictors, including rolling window, exponentially weighted moving average (EWMA), and generalized autoregressive conditional heteroscedastic (GARCH) type methods. We find that our method performs slightly better than other predictors. However, even the simplest predictors perform well for practical problems like portfolio optimization.
Everything needed to reproduce our results, together with an open source implementation of our proposed covariance predictor, is available online at
In chapter (ref) we describe some common predictors, including the one that our method builds on. We introduce our proposed covariance predictor in chapter (ref). In chapter (ref) we discuss methods for validating covariance predictors that measure both overall performance and reactivity to market changes. We describe the data we use in our first empirical studies in chapter (ref), and give the results in chapter (ref).
In the next chapters we discuss some extensions of and variations on our method, including realized covariance prediction (chapter (ref)), handling large universes via factor models (chapter (ref)), obtaining smooth covariance estimates (chapter (ref)), and using our covariance model to generate simulated returns (chapter (ref)).
\chapter{Some common covariance predictors} In this chapter we review some common covariance predictors, ranging from simple to complex, with the goal of giving context and fixing our notation. To simplify some formulas, we take $r_\tau=0$ for $\tau \leq 0$.
The rolling window predictor with window length or memory $M$ is the average of the last $M \geq n$ outer products, \[ \hat{\Sigma}_t = \alpha_t\sum_{\tau=t-M}^{t-1} r_\tau r_\tau^T, \quad t=2,3, \ldots, \] where $\alpha_t = 1/\min\{t-1,M\}$ is the normalization constant. The rolling window predictor can be evaluated via the recursion \[ \hat{\Sigma}_{t+1} = \frac{\alpha_{t+1}}{\alpha_{t}}\hat{\Sigma}_{t} + \alpha_{t+1}(r_t r_t^T-r_{t-M}r_{t-M}^T), \quad t=1,2, \ldots, \] with initialization $\hat{\Sigma}_1=0$.
For $t<n$, the rolling window covariance estimate is not full rank. To handle this, as well as to improve the quality of the prediction, we can add regularization or shrinkage, for example by adding a positive multiple of $\mathop{\bf diag}(\hat{\Sigma}_t)$ to our estimate ledoit110_honey_cov, ledoit2003improved, or approximating the predicted covariance matrix by a diagonal plus low rank matrix, as described in chapter (ref).
The exponentially weighted moving average (EWMA) estimator, with forgetting factor $\beta \in(0,1)$, is
where \[ \alpha_t=\left( \sum_{\tau=1}^{t-1}\beta^{t-1-\tau} \right)^{-1} = \frac{1-\beta}{1-\beta^{t-1}} \] is the normalization constant. The forgetting factor $\beta$ is usually expressed in terms of the half-life $H=-\log 2 / \log\beta$, for which $\beta^H =1/2$. The half-life $H$ is the number of periods when the exponential weight has decreased by a factor of two. For example, for a half-life of one year, the current observed return has twice the impact on our covariance prediction as the return observed one year ago. The EWMA predictor is widely used in practice; for example RiskMetrics suggests the forgetting factor $\beta=0.94$, which corresponds to a half-life of around 11 days menchero2011barra, longerstaey1996riskmetrics.
The EWMA covariance predictor can be computed recursively as \[ \hat{\Sigma}_{t+1} = \frac{\beta-\beta^{t}}{1-\beta^t} \hat{\Sigma}_t + \frac{1-\beta}{1-\beta^t}r_tr_t^T, \quad t=1,2, \ldots, \] with initialization $\hat{\Sigma}_1=0$. Like the rolling window predictor, the EWMA predictor is singular for $t<n$, which can be handled using the same regularization methods described above.
\paragraph{GARCH.} The generalized autoregressive conditional heteroscedastic (GARCH) predictor decomposes the return of a single asset as \[ r_t = \mu + \epsilon_t, \] where $\mu$ is the mean return and $\epsilon_t$ is the innovation, and models the innovation as \[ \epsilon_t = \sigma_t z_t, \qquad \sigma_t^2 = \omega + \sum_{\tau=1}^q a_{\tau} \epsilon_{t-\tau}^2 + \sum_{\tau=1}^p b_{\tau} \sigma_{t-\tau}^2, \] where $\sigma_t$ is the asset volatility, $z_t$ are independent $\mathcal N(0,1)$, and $q$ and $p$ (often both set to one in practice) determine the GARCH order bollerslev_1986. (Recall that we assume zero mean.) The model parameters are $\omega$, $a_1, \ldots, a_q$, and $b_1, \ldots, b_p$. Estimating the model parameters requires solving a nonconvex optimization problem cov_barrat_2022.
With $p=0$ we recover the autoregressive conditional heteroscedastic (ARCH) predictor, introduced in the seminal paper by engle_1982. This paper set the foundation for a wide variety of popular volatility and correlation predictors and earned him the 2003 Nobel Memorial Prize in Economic Sciences.
\paragraph{MGARCH.} There are several ways of extending the GARCH predictor to a multivariate or vector setting. The most popular is the dynamic conditional correlation (DCC) predictor DCC, which is a two-step approach described below.
Many other MGARCH predictors have been proposed. The most straightforward generalization from the univariate to multivariate predictors is the VEC predictor, where the covariance matrix is vectorized and each element is modeled as a GARCH process with dependencies on all other elements bollerslev_engle_1988. However, this extension requires estimating $n(n+1)(n(n+1)+1)/2 \approx n^4/2$ parameters, which can be impractical even for modest values of $n$.
Following the VEC extension of GARCH, multivariate GARCH (MGARCH) predictors have been proposed in two lines of development silvennoinen2009multivariate. The first line involves models that impose restrictions on the parameters of the VEC predictor, including DVEC bollerslev_1986, BEKK engle_kroner_1995, FF-MGARCH FF_MGARCH, O-GARCH OGARCH, and GO-GARCH GOGARCH, to name some. However, these predictors have been shown to be hard to fit and can yield inconsistent estimates brooks_2003. (These inconsistencies may not have much practical impact.) For detailed reviews of MGARCH predictors we refer the reader to silvennoinen2009multivariate, garch_survey
The second line of extensions of GARCH to vector time series models conditional covariances through separate estimates of conditional variances and correlations DCC, engle2001theoretical. In CCC Bollerslev introduced the constant conditional correlation predictor (CCC) where the individual asset volatilities are modeled as separate GARCH processes, while the correlation matrix is assumed constant and equal to the unconditional correlation matrix. This predictor was later extended to the dynamic conditional correlation (DCC) predictor where the correlation matrix is allowed to change over time DCC. The DCC model has the form \[ \Sigma_t = D_tR_tD_t, \] where $D_t$ is the diagonal matrix of standard deviations, {\it i.e.}, $(D_t)_{ii}=(\Sigma_t)_{ii}^{1/2}$, and $R_t$ is the correlation matrix associated with $\Sigma_t$.
DCC GARCH models the diagonal elements of $D_t$ as separate univariate GARCH processes as described above. The correlation matrix $R_t$ is then modeled as a constrained multivariate GARCH (MGARCH) process, {\it e.g.}, as
where $\bar{Q}$ is the unconditional correlation matrix, $a$ and $b$ are the MGARCH parameters, and $\tilde{r}_t$ are the volatility adjusted returns defined as \[ \tilde{r}_t = D_t^{-1} r_t. \] The parameters can be estimated in two steps via (quasi) maximum likelihood, but requires solving non-convex optimization problems DCC. This predictor has become a popular choice amongst MGARCH predictors due to its interpretability. Variants of the DCC predictor are widely used in finance, where it is also often used in combination with EWMA estimates. Conditional correlation predictors are easier to estimate than other multivariate GARCH predictors, and their parameters are more interpretable.
\paragraph{Iterated covariance estimation.} DCC, which separately estimates the volatilities and correlations, is closely related to the idea of iterated covariance predictors cov_barrat_2022. Iterated covariance predictors estimate the covariance matrix in multiple iterations. In a two-step iteration we first form a first covariance estimate $\hat{\Sigma}_t^{(1)}$ of the returns $r_t$, at each time $t$, and form the whitened returns \[ \tilde{r}_t = \left(\hat{\Sigma}_t^{(1)}\right)^{-1/2} r_t. \] In the second iteration we form the covariance estimate $\hat{\Sigma}_t^{(2)}$ of the whitened returns $\tilde{r}_t$. The final covariance estimate (of the returns $r_t$) is then formed as \[ \hat{\Sigma}_t = \left(\hat{\Sigma}_t^{(1)}\right)^{1/2} \hat{\Sigma}_t^{(2)} \left(\hat{\Sigma}_t^{(1)}\right)^{1/2}. \] This procedure can be iterated further, and has been shown empirically to improve the quality of the covariance estimate; see cov_barrat_2022 for details. In DCC, $\hat \Sigma^{(1)}$ is diagonal and models the volatilities; $\hat \Sigma^{(2)}$ is a correlation matrix.
Iterated EWMA (IEWMA) was proposed by DCC and is analogous to DCC GARCH but with EWMA estimates of the volatilities and correlations instead of GARCH. Engle proposed IEWMA as an efficient alternative to the DCC GARCH predictor, although he did not refer to it as IEWMA; we use this term to emphasize its connection to iterated whitening, as proposed in cov_barrat_2022. Specifically, IEWMA can be viewed as an iterated whitener, where we first use a diagonal whitener (which estimates the volatilities) and then a full matrix whitener (which estimates the correlations). This is analogous to the two-step iterated covariance predictor where $\hat{\Sigma}_t^{(1)}$ is the diagonal matrix of squared volatility estimates and $\hat{\Sigma}_t^{(2)}$ estimates the correlation matrix of the volatility adjusted returns.
First we form an estimate of the volatilities $\hat\sigma_t=\mathop{\bf diag}(\hat{\Sigma}_t)^{1/2}$ using EWMA predictors for each asset. We denote the half-life of these volatility estimates as $H^\text{vol}$. We then form the marginally standardized returns as
where $\hat D_t = \mathop{\bf diag}(\hat \sigma_t)$. These vectors should have entries with standard deviation near one. It is common practice to winsorize the standardized returns; a good rule of thumb is to clip $\tilde r_t$ at $\pm 4.2$, which corresponds to clipping $r_t$ at $\pm 4.2 \hat \sigma_t$.
Then we form a EWMA estimate of the covariance of $\tilde r_t$, which we denote as $\tilde R_t$, using half-life $H^\text{cor}$ for this EWMA estimate. (We use the superscript `cor' since the diagonal entries of $\tilde R_t$ should be near one, so $\tilde R_t$ is close to a correlation matrix.) From $\tilde R_t$ we form its associated correlation matrix $\hat R_t$, {\it i.e.}, we scale $\tilde R_t$ on the left and right by a diagonal matrix with entries $(\tilde R_t)_{ii}^{-1/2}$. Since the diagonal entries of $\tilde R_t$ should be near one, $\tilde R_t$ and $\hat R_t$ are not too different.
Our IEWMA covariance predictor is \[ \hat{\Sigma}_t = \hat D_t \hat R_t \hat D_t, \quad t=2,3, \ldots. \] This is the covariance predictor proposed in DCC; replacing $\hat R_t$ with $\tilde R_t$ we obtain the iterated whitener proposed by Barratt and Boyd in cov_barrat_2022. As mentioned above, they are typically quite close.
It is common to choose the volatility half-life $H^\text{vol}$ to be smaller than the correlation half-life $H^\text{cor}$. The intuition here is that we can average over fewer past samples when we predict the $n$ volatilities $\hat{\sigma}_t$, but need more past samples to reliably estimate the $n(n-1)/2$ off-diagonal entries of $\hat R_t$. Empirical studies on real return data confirm that choosing a faster volatility half-life than correlation half-life yields better estimates.
\chapter{Combined multiple iterated EWMAs} In this chapter we introduce a novel covariance predictor, which we call combined multiple iterated EWMAs, for which we use the acronym CM-IEWMA. The CM-IEWMA predictor is constructed from a modest number of IEWMA predictors, with different pairs of half-lives, which are combined using dynamically varying weights that are based on recent performance.
The CM-IEWMA predictor is motivated by the idea that different pairs of half-lives may work better for different market conditions. For example, short half-lives perform better in volatile markets, while long half-lives perform better for calm markets where conditions are changing slowly.
We first describe the idea in a general setting. We start with $K$ different covariance predictors, denoted $\hat \Sigma^{(k)}_t$, $k=1, \ldots, K$. These could be any of the predictors described above, or predictors of the same type with different parameter values, {\it e.g.}, half-lives (for EWMA) or pairs of half-lives (for IEWMA). In some contexts these different predictors are referred to as a set of $K$ experts hastie2009elements, jordan1994hierarchical.
We denote the Cholesky factorizations of the associated precision matrices $(\hat{\Sigma}_t^{(k)})^{-1}$ as $\hat L_t^{(k)}$, {\it i.e.}, \[ \left(\hat{\Sigma}_t^{(k)}\right)^{-1} = \hat L_t^{(k)} ( \hat L_t^{(k)} )^T, \quad k=1, \ldots, K, \] where $\hat L_t^{(k)}$ are lower triangular with positive diagonal entries. We will combine these Cholesky factors with nonnegative weights $\pi_1,\ldots, \pi_K$ that sum to one, to obtain
From this we recover the weighted combined predictor
We will see below why we combine the Cholesky factors of the precision matrices, and not the covariance or precision matrices themselves.
The log-likelihood (ref) can be expressed in terms of the Cholesky factor of the precision matrix $\hat L_t$ as \[ l_t (\hat{\Sigma}_t) = -(n/2)\log(2\pi) + \sum_{i=1}^n \log \hat{L}_{t,ii} - (1/2)\|\hat{L}^T_tr_t\|^2_2, \] where $\| \cdot \|_2$ denotes the Euclidean norm. This is a concave function of the weights $\pi\in {\mbox{\bf R}}_+^K$ boyd2004convex.
We choose the weights at time $t$ as the solution of the convex optimization problem
with variables $\pi_1, \ldots, \pi_K$, where $N$ is the look-back, $\mathbf 1$ denotes the vector with entries one, and $\geq$ between vectors means entrywise. In words: we choose the (mixture) weights in each period so as to maximize the average log-likelihood of the combined prediction over the trailing $N$ periods. The problem (ref) is convex, and can be solved very quickly and reliably by many methods boyd2004convex. The covariance predictor is then recovered using (ref) and (ref).
The look-back $N$ is a parameter that can be adjusted to give good performance. Numerical experiments suggest that the predictor is not very sensitive to the choice of $N$, and that a choice $N=10$ seems to work well for asset universes up to a few hundred assets.
We mention several extensions of the weight problem (ref). First, we can add one prediction which is diagonal, using any estimates of the volatilities (including constant). This gives us shrinkage, automatically chosen. We can also add a constraint or objective term that encourages the weights to vary smoothly over time, as discussed more in chapter (ref).
The CM-IEWMA predictor is a special case of the dynamically weighted prediction combiner described above, where the $K$ predictions are each IEWMA, with different pairs of half-lives $H^\text{vol}$ and $H^\text{cor}$.
\chapter{Evaluating covariance predictors} There are several ways of evaluating a covariance predictor, often divided into two categories, direct and indirect patton2009evaluating, ANDERSEN2006777. Direct methods use a proxy for the true covariance matrix to evaluate the predictor, while indirect methods use the covariance predictor on tasks of interest, such as portfolio construction or portfolio tracking.
Popular direct methods are the Mincer-Zarnowitz (MZ) regression and its variants, based on statistical tests of the regression coefficients of a predicted variable on an observed variable (or in the case of variance and covariance, a proxy for the observed variable) mincer1969evaluation, theil1961economic. Direct methods also include the comparison between different predictors in terms of some loss function. Common loss functions are the mean squared error (MSE) and quasi-likelihood (QLIKE) patton2011volatility, patton2009evaluating. To select good models, the model confidence set (MCS) is usually used hansen2011model, or the Ledoit--Wolf test ledoit2008robust to compare Sharpe ratios.
Indirect methods use applications to rank covariance predictors, and include the minimum variance and mean-variance portfolios, as well as portfolio tracking tasks.
The difference in performance between various predictors can also be evaluated using statistical tests. For a more detailed discussion of both direct and indirect methods, we refer the reader to patton2009evaluating.
In this chapter we discuss several evaluation metrics for covariance predictors. The first three metrics are direct, and include the mean squared error and two metrics based on a statistical measure, the log-likelihood under a Gaussian distribution. The remaining metrics judge a covariance predictor by the performance of a portfolio using a method that depends on a covariance matrix. We are mainly interested in illustrating how simple methods can perform just as well as or better than more complex ones, rather than finding optimal predictors in a statistical sense. Therefore we look at the absolute performance of covariance predictors on these metrics.
The mean squared error (MSE) is a common metric for evaluating a covariance predictor $\hat{\Sigma}_t$, defined as \[ \frac{1}{T} \sum_{t=1}^T \|r_tr_t^T - \hat{\Sigma}_t\|_F^2, \] {\it i.e.}, the average squared Frobenius norm of the difference between the realized (rank one) covariance matrix $r_tr_t^T$ and the covariance predictor $\hat{\Sigma}_t$. Lower values of MSE are better. One variation on the MSE error assumes that $\hat{\Sigma}_t$ is constant over some number of time periods and replaces the rank one realized covariance $r_tr_t^T$ with an average of the rank one terms over the periods, {\it i.e.}, the realized empirical covariance.
A natural way of judging a covariance predictor is via its average log-likelihood on realized returns, \[ \frac{1}{2T}\sum_{t=1}^T\Big( -n\log(2\pi) - \log\det \hat{\Sigma}_t - r_t^T \hat{\Sigma}_t^{-1}r_t\Big), \] with larger values being better. This metric can be used to compare different predictors.
To understand the performance of a covariance predictor over time and changing market conditions, we can examine the average log-likelihood over periods such as quarters, and look at the distribution of quarterly average log-likelihood values. We are particularly interested in poor, {\it i.e.}, low values.
Recall that the best constant predictor, in terms of the log-likelihood, is the empirical sample covariance \[ \Sigma^\text{emp} = \frac{1}{T} \sum_{t=1}^T r_t r_t^T , \] with value \[ \frac{1}{2}\big(-n(\log(2\pi)+1) - \log\det \Sigma^\text{emp}\big). \] For any other constant $\Sigma \in {\mbox{\bf S}}_{++}^n$, the log-likelihood is lower than the log-likelihood of $\Sigma^\text{emp}$. We define the average log-likelihood regret as the average log-likelihood of the (constant) empirical covariance, minus the average log-likelihood of the covariance predictor. The regret is a measure of how much the covariance predictor $\hat{\Sigma}_t$, $t=1, \ldots, T$, underperforms the best possible constant covariance predictor ({\it i.e.}, the sample covariance matrix). The term regret comes from the field of online optimization; see, {\it e.g.}, zinkevich2003_regret, Mokhtari2016_regret, hazan2007_regret, hazan2016introduction.
We want our covariance predictor to have small regret. The regret is typically positive, but it can be negative, {\it i.e.}, our time-varying covariance can have higher log-likelihood than the best constant one. The regret is not any more useful than the log-likelihood when comparing predictors over one time interval, since it simply adds a constant and switches the sign. But it is interesting when we compute the regret over multiple periods, like months or quarters. The regret over multiple quarters removes the effect of the log-likelihood of the empirical covariance varying due to changing market conditions, and allows us to assess how well the covariance predictor adapts.
We can also judge the performance of a covariance predictor by the investment performance of portfolio construction methods that depend on the estimated covariance matrix. As with log-likelihood or log-likelihood regret, we can examine the portfolio performance in periods such as quarters, to see how evenly the performance is spread over time.
One obvious metric of interest is how close the ex-ante and realized portfolio volatilities are. The metrics described above, MSE, log-likelihood, and log-likelihood regret, are agnostic to the portfolio; with specific real portfolios we can see how well our covariance predictors predict portfolio volatility.
We will assess a covariance predictor using five simple portfolio construction methods. The first is an equally weighted (or $1/n$) portfolio, which does not by itself depend on the covariance, but does when we adjust it with cash to achieve a given ex-ante risk. The second, third, and fourth portfolios depend only on the covariance matrix. They are minimum variance, risk parity, and maximum diversification portfolios. For an in depth discussion of these portfolios, see braga2015risk. The last portfolio we consider is a mean-variance portfolio, using a very simple mean estimator.
For each portfolio we look at four metrics: realized return, volatility, Sharpe ratio, and maximum drawdown. The returns, volatilities, and Sharpe ratios are reported in annualized values. The Sharpe ratio is defined as the ratio of the excess return (over the risk-free rate), divided by the volatility of the excess return, \[ \frac{\frac{1}{T}\sum_{\tau=1}^T (r^p_t-r^{\text{rf}}_t)}{\Big(\frac{1}{T}\sum_{\tau=1}^T \big( r^p_t - \frac{1}{T}\sum_{\tau=1}^T r^p_t\big)^2\Big)^{1/2}}, \] where $r^p_t$ and $r^{\text{rf}}_t$ are the portfolio and risk-free returns at time $t$. The maximum drawdown is defined as \[ \max_{1 \leq t_1 < t_2\leq T} \frac{V^p_{t_1}}{V^p_{t_2}} - 1, \] where \[ V^p_t = V_0 (1+r^p_1)(1+r^p_2) \cdots (1+r^p_t) \] is the portfolio value at time $t$ (with returns re-invested), starting with value $V_0>0$.
In addition to portfolio performance, we can also examine how well the covariance prediction predicts the portfolio volatility. We compare the realized or ex-post portfolio volatility \[ \left( \frac{1}{T} \sum_{t=1}^T (r_t^Tw_t)^2 \right)^{1/2}, \] to the predicted or ex-ante portfolio volatility \[ \left( \frac{1}{T} \sum_{t=1}^T w_t^T \hat{\Sigma}_t w_t\right)^{1/2}, \] where $w_t\in{\mbox{\bf R}}^n$ are the portfolio weights. This directly measures the ability of the estimated covariance matrix to predict portfolio risk.
\paragraph{Equal weight portfolio.} We take the equal weight or $1/n$ portfolio with $w= (1/n)\mathbf 1$. This portfolio does not depend on the covariance $\hat{\Sigma}_t$, but when we mix it with cash, as described below, it will.
\paragraph{Minimum variance portfolio.} The (constrained) minimum variance portfolio is the solution of the convex optimization problem \[
\] with variable $w\in {\mbox{\bf R}}^n$, where $L_{\max} \geq 1$ is a leverage limit, and $w_{\min}$ and $w_{\max}$ are lower and upper bounds on the weights, respectively.
\paragraph{Risk-parity portfolio.} The portfolio return volatility $\sigma(w) = (w^T \hat{\Sigma}_t w)^{1/2}$ can be broken down into a sum of volatilities (risks) associated with each asset as \[ \frac{\partial \log \sigma(w)}{\partial w_i} = \frac{\partial \sigma(w)}{\sigma(w)} \frac{w_i}{\partial w_i} =\frac{w_i(\hat{\Sigma}_t w)_i} {w^T \hat{\Sigma}_t w}, \quad i=1,\dots,n. \] The risk parity portfolio is the one for which these volatility attributions are equal Qian119. This portfolio can be found by solving the convex optimization problem cvx_book_additional, \[ \mbox{minimize} \quad (1/2)x^T\hat{\Sigma}_t x - \sum_{i=1}^n (1/n)\log x_i, \] with variable $x$, and then taking $w=x^{\star}/(\mathbf 1^Tx^{\star})$.
\paragraph{Maximum diversification portfolio.} The diversification ratio of a long-only portfolio ({\it i.e.}, one with $w \geq 0$) is defined as \[ D(w) =\frac{\hat \sigma_t^T w}{(w^T\hat{\Sigma}_t w)^{1/2}}. \] The diversification ratio tells us how much higher the portfolio volatility would be if all assets were perfectly correlated. The maximum diversification portfolio is the portfolio $w$ that maximizes $D(w)$, possibly subject to constraints choueifaty2008toward. Like the risk-parity portfolio, the maximum diversification portfolio can be found via convex optimization. We let $x^\star$ denote the solution of the convex optimization problem cvx_book_additional \[
\] with variable $x$. The maximum diversification portfolio is $w=x^{\star}/\mathbf 1^T x^{\star}$.
\paragraph{Volatility control with cash.} We mix each of the four portfolios described above with cash to achieve a target value of ex-ante volatility $\sigma^\mathrm{tar}$. To do this we start with the portfolio weight vector $w_t$, and compute its ex-ante volatility $\sigma_t = (w_t^T \hat{\Sigma}_t w_t)^{1/2}$. Then we add a cash component so that the overall ex-ante volatility equals our target, {\it i.e.}, we use the $(n+1)$ weights (with the last component denoting cash) \[ \left[
\right], \qquad \theta = \frac{\sigma^tar}{\sigma_t}. \] This portfolio will have ex-ante volatility $\sigma^\mathrm{tar}$. Note that the cash weight can be either positive (when it dilutes the portfolio volatility) or negative (when it leverages the portfolio volatility to the desired level). The target volatility $\sigma^\text{tar}$ should be chosen so as to avoid portfolios that are either too diluted or too leveraged.
\paragraph{Mean variance portfolio.} The last portfolio we consider is a basic mean-variance portfolio, defined as the solution of the convex optimization problem \[
\] with variable $w$, where $\hat{r}_t$ is the predicted mean return vector at time $t$. The vector $w$ gives the weights of the non-cash assets and $c$ denotes the cash weight. The non-cash and cash weights are limited by $w_{\min}, w_{\max}$ and $c_{\min}, c_{\max}$, respectively. This portfolio does not need cash dilution, since it includes cash in its construction. (If $\sigma^\text{tar}$ is chosen appropriately, it will have ex-ante risk $\sigma^\text{tar}$.) The mean-variance portfolio depends not only on a covariance estimate, but also a return estimate. For this we use one of the simplest possible return estimates, a EWMA of the realized returns.
\chapter{Data sets and experimental setup} We illustrate our method on three different data sets: a set of $49$ industry portfolios, a set of $25$ stocks, and a set of $5$ factor returns, each augmented with cash (with the historical risk-free interest rate). For each data set we show results for six covariance predictors. Everything needed to reproduce the results is available online at
\paragraph{Industry portfolios.} The first data set consists of the daily returns of a universe of $n=49$ daily traded industry portfolios, shown in table (ref), along with cash. The data set spans July 1st 1969 to December 30th, 2022, for a total of 13496 (trading) days. The data was obtained from the Kenneth French Data Library french_data_lib.
\paragraph{Stocks.} The second data set consists of the daily returns of $n=25$ stocks and cash. The stocks were chosen to be the 25 largest stocks in the S&P 500 at the beginning of 2010, listed in table (ref). This data set spans January 4th 2010 to December 30th, 2022, for a total of 3272 (trading) days. The stock data was attained through the Wharton Research Data Services (WRDS) portal WRDS.
\paragraph{Factor returns.} The third data set consists of daily returns of the five Fama-French factors taken from the Kenneth French Data Library french_data_lib, shown in table (ref). The data set spans July 1st 1963 to December 30th, 2022, for a total of 14979 (trading) days.
\paragraph{Cumulative returns.} In figure (ref) we show the cumulative returns of the five factors, and the cumulative returns of five assets chosen from each of the industry and stock data sets.
For each data set we evaluate six covariance predictors, described below.
All the parameters above ({\it e.g.}, half-lives) are chosen as reasonable values that give good overall performance for each predictor. The results are not sensitive to these choices.
For our experiments we use the first two years (500 data points) of each data set to fit the MGARCH predictor and initialize the other predictors. (After this initial MGARCH fit, we re-estimate its parameters annually.) Hence, the evaluation period for our experiments below ranges from June 24th 1971 to December 30th, 2022, for the industry portfolios, from December 28th, 2011, to December 30th, 2022, for the stock portfolios, and from June 28th 1965 to December 30, 2022, for the factor portfolios.
\chapter{Results}
Figure (ref) shows the weights for each of the five components of the CM-IEWMA predictors, averaged yearly, for the three data sets.
We can see how the predictor adapts the weights depending on market conditions. Substantial weight is put on the slower (longer half-life) IEWMAs most years. During and following volatile periods like the 2000 dot.com bubble or 2008 market crash, we see a big increase in weight on the faster IEWMAs. We can illustrate these changes in weights in response to market conditions via the effective half-life of the CM-IEWMA, defined as the weighted average of the five (longer) half-lives, shown in figure (ref), averaged yearly.
Table (ref) shows the average, standard deviation, and maximum of the MSE computed over distinct quarters for the six covariance predictors on the three data sets (with lower being better for all three metrics).
CM-IEWMA and MGARCH do better than the other predictors on all metrics over all data sets, with MGARCH doing slightly better on the industry data and CM-IEWMA slightly better on the stock data. Interestingly, on the factor data set, the CM-IEWMA predictor does better than the prescient predictor.
Figure (ref) shows the average quarterly log-likelihood for the different covariance predictors over the evaluation period. Not surprisingly, the prescient predictor does substantially better than the others. The different predictors follow similar trends, with even the prescient predictor experiencing a drop in log-likelihood during market turbulence. Close inspection shows that the CM-IEWMA and MGARCH predictors almost always have the highest log-likelihood in each quarter.
Figure (ref) shows the average quarterly log-likelihood regret for the different covariance predictors over the evaluation period. Clearly, CM-IEWMA and MGARCH perform best in volatile markets. Figure (ref) illustrates the difference between CM-IEWMA and MGARCH. As seen, CM-IEWMA consistently has lower regret on the industry and stock data sets, while they perform similar on the factor data. More precisely, CM-IEWMA has lower regret than MGARCH in 87% of the quarters for the industry data, 71% for the stock data, and 51% for the factor data.
Table (ref) illustrates the differences in regret further, by showing the average, standard deviation, and the maximum of the average quarterly regret. As we can see, the average quarterly regret is lower for CM-IEWMA than for the other predictors. The regret is also more stable for CM-IEWMA, as the standard deviation is lower. Finally, the maximum average quarterly regret is also lower for CM-IEWMA than for the other predictors. These results are most prominent on the industry and stock data, while MGARCH does similar on the factor data.
Figure (ref) gives a final illustration of these results, by showing the cumulative distribution functions of the average quarterly regret for the different covariance predictors.
Clearly, CM-IEWMA has the lowest regret on the industry and stock data set, and MGARCH does similar on the factor data.
In this section we evaluate the covariance predictors on the portfolios described in \S(ref). In the minimum variance and mean-variance portfolios, we use $L_\text{max}=1.6$ (which corresponds to 130:30 long:short), $w_\text{min} = -0.1$ and $w_\text{max}=0.15$ for the industry and stock return portfolios, and $w_\text{min} = -0.3$ and $w_\text{max}=0.4$ for the factor return portfolio. We use target (annualized) volatilities of 5%, 10%, and 2% for the industry, stock, and factor return portfolios, respectively.
For the mean-variance portfolio, our estimated returns are EWMAs of the trailing realized returns. For the industry and stock data we use 250-day half-life EWMAs, winsorized at the 40th and 60th percentiles (cross-sectionally), and for the factor data a 63-day half-life EWMA (not winsorized).
\paragraph{Equal weight portfolio.} Table (ref) shows the metrics for the equal weight portfolio. All predictors track the volatility targets well. MGARCH attains the highest Sharpe ratios, although the results are very close. The drawdowns are also very similar for all predictors, but MGARCH and CM-IEWMA seem slightly better than the rest.
\paragraph{Minimum variance portfolio.} Table (ref) shows the metrics for the minimum variance portfolio. For the factor data set, MGARCH does best. On the industry and stock data sets, the three EWMA-based predictors track the volatility target fairly well, while RW and MGARCH underestimate volatility. CM-IEWMA and MGARCH both attain a high Sharpe ratio. However, we note that the high Sharpe ratio for MGARCH, as compared to the other predictors, is a consequence of the high volatility. Finally, CM-IEWMA seems to consistently attain a lower drawdown than the other predictors, although the other EWMA-based approaches also do well.
To illustrate how the minimum variance trading strategy has evolved over time, we show the yearly annualized Sharpe ratios for the CM-IEWMA predictor in figure (ref). We can see that the Sharpe ratio achieved by the minimum variance portfolio decreases over time for the industry and stock data sets, with a small upward trend for the factor data set.
\paragraph{Risk parity portfolio.} The results for the risk-parity portfolio are shown in table (ref). Overall the results are similar for the various predictors. There is very little that separates the predictors on the industry data set. On the stock data, CM-IEWMA and MGARCH attain the highest Sharpe ratios and lowest drawdowns. On the factor data set, MGARCH has the best overall performance.
\paragraph{Maximum diversification portfolio.} The maximum diversification portfolio results are illustrated in table (ref). On the industry and stock data sets, CM-IEWMA and MGARCH do best in terms of Sharpe ratio, drawdown, and tracking the volatility target. On the factor data set, MGARCH does best overall.
\paragraph{Mean variance portfolio.} The results for the mean-variance portfolio are given in table (ref). On the industry data set all predictors underestimate volatility. The results are similar across predictors, with CM-IEWMA and MGARCH performing slightly better than the rest in terms of Sharpe ratio and drawdown. On the stock data set, CM-IEWMA seems to do best overall. On the factor data set, the results are almost identical between predictors.
Since we use simple EWMA return predictors, we can expect the mean variance portfolio performance to vary over time. Intuitively it should be better on historical data than more recent data. To illustrate this, figure (ref) shows the yearly annualized Sharpe ratios for the three portfolios.
There is a clear downward trend in the Sharpe ratios for the industry and factor data sets, illustrating the difficulty of predicting returns in recent years. This can be compared to the minimum variance portfolios (figure (ref)) that have a more stable performance over time, and notably do not depend on a mean estimate.
In terms of log-likelihood and regret, CM-IEWMA performs best, followed by MGARCH, which performs better than the simpler covariance predictors. In downstream portfolio optimization experiments, CM-IEWMA and MGARCH again perform better than the other predictors, although in many cases not by much. In these experiments there is more variation in the results, partly explained by the difference between our prediction (of a covariance matrix) and our metrics (such as return, risk, drawdown). Even the simplest covariance predictors do a reasonable job of predicting the portfolio risk.
\chapter{Realized covariance} We have so far focused on predicting the covariance matrix of asset returns using historical return data, {\it i.e.}, we predict $\hat{\Sigma}_t$ from $r_1, \ldots, r_{t-1}$. In this chapter we consider the use of additional data, specifically, intraperiod returns. As an example, suppose the period is (trading) days. The methods described in previous chapters predict the covariance of the daily return from previous daily returns. In so-called realized covariance, we predict the daily return covariance using intraday returns. Instead of single period returns $r_t$, we have multiple returns associated with period $t$. It is not surprising that using multiple realized returns for each period, instead of just one, can improve our covariance estimates.
Recent literature has shown that realized volatility and correlation measurements (based on high-frequency intraperiod data) can improve performance over traditional predictors that rely on a single realization per period. hansen2012realized extend the univariate GARCH model to the joint modeling of returns and realized measures of volatility, and show empirically that this improves performance over the standard GARCH model. In bauwens2012dynamic a multivariate realized GARCH model is proposed. More recently, bollerslev2020multivariate propose a realized semicovariance GARCH model to allow for nuanced responses to positive and negative return shocks.
In this chapter we show that the dynamically weighted prediction combiner of \S(ref) readily handles multiple realized returns per period. For simplicity we will assume each period has the same number of intraperiod returns, equally spaced in time. We redefine the return vector to be a return matrix $r_t\in {\mbox{\bf R}}^{n\times m}$ with columns that are the $m$ intraperiod return vectors, for times $t=1,\ldots,T$. The realized covariance at time $t$ is defined as \[ C_t = r_tr_t^T, \] the same formula for the realized return when $r_t$ is a single (vector) return. The realized covariance matrix $C_t$ has rank $m$ when the $m$ return vectors are linearly independent and $m\leq n$; this can be compared to the realized covariance when we do not have intraperiod returns, which is rank one.
The dynamically weighted prediction combiner of \S(ref) readily handles multiple realized covariance predictors. \paragraph{Realized EWMA.} We define the realized EWMA (REWMA) predictor as \[ \hat{\Sigma}_t = \alpha_{t} \sum_{\tau=1}^{t-1}\beta^{t-1-\tau} C_\tau, \quad t=2,3, \ldots, \] where $C_\tau$ is the realized covariance at time $\tau$, $\beta \in(0,1)$ is the forgetting factor, and $\alpha_t$ is the normalizing constant; see \S(ref) for details. This is the same formula as the usual EWMA covariance, with one return per period, given in (ref), with $r_t$ extended to be a matrix of multiple returns.
\paragraph{Combined multiple realized EWMAs.} The combined multiple realized EWMA (CM-REWMA) predictor starts with a set of $K$ REWMA predictors $\hat{\Sigma}_t^{(k)}$ with half-lives $H^{(k)}$, $k=1,\ldots,K$, and combines them using the dynamically weighted prediction combiner of \S (ref).
\paragraph{Data set.} We consider a universe of $n=39$ assets with five-minute intraday returns corresponding to $m=77$. The assets were taken as a subset of those used by pelger2020understanding, and are available at pelger-website. The data set spans January 2nd 2004 to December 30th 2016, for a total of 252021 data points over 3273 trading days. We list the assets in table (ref).
\paragraph{Four covariance predictors.} We evaluate four covariance predictors, described below.
\paragraph{CM-REWMA component weights.} Figure (ref) shows the component weights for the CM-REWMA predictor, averaged annually.
The weights are fairly stable over time but a weight shift toward faster changing EWMAs is seen in 2008, during the financial crisis.
\paragraph{MSE.} The average, standard deviation, and maximum MSEs, computed over distinct quarters for the four covariance predictors, are given in table (ref).
The REWMA and CM-REWMA do slightly better than the CM-IEWMA predictor, but overall there is not a big difference between the predictors.
\paragraph{Regret.} Figure (ref) shows the average regret over distinct quarters for the CM-IEWMA, REWMA, and CM-REWMA predictors.
The CM-REWMA predictor has the lowest regret in almost all quarters. It has lower regret than the REWMA predictor in 41 out of the 50 quarters, and lower regret than the CM-IEWMA predictor in 39 out of the 50 quarters.
Finally, figure (ref) shows the cumulative distribution functions of the average quarterly regret for the different covariance predictors.
CM-REWMA has lower regret than both the CM-IEWMA and REWMA predictors, while REWMA has lower regret than CM-IEWMA.
\paragraph{Portfolio performance.} Table (ref) shows the portfolio metrics for five different portfolio construction methods.
CM-REWMA does better than, or as well as, REWMA on almost all metrics, and better than CM-IEWMA for all portfolios. However, the difference on portfolio tasks is not large.
\paragraph{Summary.} The results above show that using realized covariance, {\it i.e.}, intraperiod returns instead of just one return per period, gives covariance estimates that are a bit better than those obtained using only one return per period.
\chapter{Large universes} In a practical setting we often encounter a larger number of assets than considered in the previous chapters, which has led to extensive research in high-dimensional covariance estimation. One challenge in large dimensions is ensuring positive definiteness of the covariance matrix, in particular with model-based approaches such as MGARCH garch_survey. Several techniques have been proposed for estimating MGARCH models in large dimensions; see, {\it e.g.}, engle2019large, de2021factor, de2022large. Others have focused on estimating realized covariance matrices in high dimensions; see, {\it e.g.}, oh2016high, vassallo2021dcc, hautsch2015high, debrito2018forecasting, fan2016incorporating, ait2017using. For a detailed review of recent developments in high-dimensional covariance estimation, we recommend bauwens2023modeling.
The methods described in previous chapters can be adapted to handle large universes of assets, say $n$ larger than $100$ or so. In this chapter we describe two closely related methods for improving the performance with large $n$. Both methods end up modeling $\hat \Sigma_t$ as a low rank plus diagonal matrix, in so-called factor form. Before describing these methods, we mention that evaluating log-likelihood regret is complicated with large $n$. For the empirical covariance to be nonsingular (which is needed to evaluate the regret), we need at least $n$ periods; for daily returns with $n=1000$, this amounts to four years. Even if we have $n$ periods of data, we would only be able to evaluate the regret a few times. For example, with $n=1000$ (four years) we need at least 40 years of data to compute the average regret over 10 distinct periods. The log-likelihood, however, can still be evaluated over fewer than $n$ periods.
In practice most return covariance matrices for large universes are constructed from factors, with the model \[ r_t = F_t f_t + z_t, \quad t=1,2, \ldots, \] where $F_t \in {\mbox{\bf R}}^{n \times k}$ is the factor exposure matrix, $f_t \in {\mbox{\bf R}}^k$ is the factor return vector, $z_t \in {\mbox{\bf R}}^n$ is the idiosyncratic return, and $k$ is the number of factors, typically much smaller than $n$. The factor returns are constructed or found by several methods, such as principal component analysis (PCA), or by hand; see, {\it e.g.}, bai2008large, bai2003inferential, lettau2020estimating, lettau2020factors, pelger2022factor, pelger2022interpretable, fama1993common, fama1992cross. Thus we assume that the factor returns are known. Given the factor returns, the rows of the factor exposure matrix are typically found by least squares regression over a rolling or exponentially weighted window cochrane2009asset. The idiosyncratic returns $z_t$ are then found as the residuals in this least squares fit. The factor returns $f_t$ are modeled as $\mathcal N(0,\Sigma_t^\mathrm{f})$, and the idiosyncratic returns $z_t$ are modeled as $\mathcal N(0,E_t)$, where $E_t$ is diagonal. It is also assumed that the factor returns and idiosyncratic returns are independent across time and of each other.
We end up with a covariance matrix in factor form, {\it i.e.}, rank $k$ plus diagonal,
We can easily use the methods described above with a factor model. Simply predict the factor return covariance $\hat \Sigma^\text{f}_t$ (using the factor returns $f_t$) and the idiosyncratic variances $\hat E_t$ (using the entries of $z_t$), using the methods described in this monograph, and then form the covariance estimate \[ \hat{\Sigma}_t = F_t\hat{\Sigma}^{\text{f}}_tF_t^T + \hat{E}_t. \]
The factor model (ref) can be written in a simpler form as
with $\tilde F_t = F_t(\Sigma_t^\text{f})^{1/2}$. This form does not include a factor covariance $\Sigma^\text{f}$, or equivalently, assumes $\Sigma^\text{f}_t =I$, {\it i.e.}, the factors are independent with standard deviation one. (The associated factors are called whitened factors.) We will use the factor model form (ref) in the sequel.
The factor model (ref) has parameters $\tilde F_t$ and $E_t$, which all together include $nk+n$ scalar parameters. (Some of these are redundant; for example we can insist without loss of generality that $F$ is lower triangular.) The factor model contains substantially fewer scalar parameters than a generic $n \times n$ covariance matrix, which contains $n(n+1)/2$ scalar parameters.
The smaller number of parameters is not the only reason for using a factor model. Another is that it often gives better covariance estimates. We can think of the low rank plus diagonal structure as regularization, which can improve out-of-sample performance. In addition, the low rank plus diagonal structure can be exploited in portfolio construction, bringing the computational complexity down from $O(n^3)$ to $O(nk^2)$ operations boyd2004convex. This makes portfolio optimization with $n=1000$ assets and $k=50$ factors extremely fast, and makes possible optimization of portfolios with much larger values of $n$.
In this section we consider the problem of fitting a given covariance matrix $\Sigma$ by one in factor form, $\hat \Sigma = FF^T+ E$, where $F \in {\mbox{\bf R}}^{n \times k}$. This corresponds to the model $r=Ff+z$, with (factor return) $f \sim \mathcal N(0,I)$, and (idiosyncratic return) $r \sim \mathcal (0,E)$, with $E$ diagonal. We let $\theta = (F,E)$ denote the parameters of our factor form model.
We seek $F\in {\mbox{\bf R}}^{n \times k}$ and diagonal $E\in {\mbox{\bf R}}^{n \times n}$ (with positive diagonal entries) that minimize the Kullback-Liebler (KL) divergence between $\mathcal N(0,\Sigma)$ and $\mathcal N(0,\hat \Sigma)$,
The KL divergence can be expressed in terms of the average log-likelihood of $\mathcal N(0,\hat\Sigma)$ under $\mathcal N(0,\Sigma)$ as
where $\ell_{\hat{\Sigma}}(r)$ is the log-likelihood of $r$ under $\mathcal N(0,\hat \Sigma)$. Hence minimizing the KL-divergence (ref) is equivalent to maximizing the expected log-likelihood (ref) of $r$ under the model $\mathcal N(0,\hat \Sigma)$.
\paragraph{Solution via EM.} We can use the expectation-maximization (EM) algorithm to approximately minimize $\mathcal K(\Sigma,FF^T+E)$ over $F$ and $E$ dempster1977maximum,rubin1982algorithms. Usually EM is used to fit a factor model to data, {\it i.e.}, samples; here we use it to fit a given Gaussian distribution $\mathcal N(0,\Sigma)$. The method described below was suggested and derived by Emmanuel Cand\a`es. We are not aware of its appearance in prior literature. A forthcoming paper on this method will include more detail and applications.
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}
The EM algorithm is an iterative method for maximizing (ref). Each iteration consists of two steps: the expectation or E-step, and the maximization or M-step. We use the conventional symbols used to describe EM, and use subscript $j=1,2, \ldots$ to denote iteration number. (A good method for initializing the EM algorithm is provided below.)
\paragraph{E-step.} In the E-step, we find the expected log-likelihood under the current estimate of the parameters $\theta_{j} = (F_j,E_j)$, over the true distribution of $r$:
where $p_{\theta_j}(f \mid r)$ is the density of the conditional distribution of $f$ under the parameter estimates at iteration $j$, and $\ell_{\theta}(r,f)$ is the log likelihood of the joint distribution with variable $\theta=(F,E)$.
With our factor model the complete log-likelihood of $(r,f)$ is
The conditional distribution of $f\mid r$ under $\theta_j$ is bishop2006pattern \[ f \mid r ~\sim~ \mathcal{N}(B_j r, G_j), \] where
Hence, (ref) becomes, up to an additive constant,
where
\paragraph{M-step.} In the M-step (ref) is maximized with respect to $\theta$ to obtain the updated parameters: \[ \theta_{j+1} = \mathop{\rm argmax}_\theta \,\, Q(\theta \, || \, \theta_j). \] The maximizer of (ref) is rubin1982algorithms
where the inner $\mathop{\bf diag}$ extracts the diagonal of its (matrix) argument, and the outer $\mathop{\bf diag}$ creates a diagonal matrix from its (vector) argument.
\paragraph{EM iteration.} The EM iteration has the form
where $C_{rr}$, $C_{rf}$, and $C_{ff}$ come from (ref) and (ref).
\paragraph{Initialization.} To initialize the EM algorithm we use the following method, based on low rank approximation via eigendecomposition. We work with the correlation matrix of $\Sigma$, denoted \[ R = \mathop{\bf diag}(\sigma)^{-1} \Sigma \mathop{\bf diag}(\sigma)^{-1}, \] where $\sigma= \mathop{\bf diag}(\Sigma)^{1/2}$ (entrywise). First we express $R$ in its eigendecomposition $R = \sum_{i=1}^n \lambda_i q_iq_i^T$, with $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n$. We then form the rank $k$ approximation \[ \widehat R = \sum_{i=1}^k \lambda_i q_iq_i^T. \] We only need to compute the $k$ dominant eigenvectors and eigenvalues, which can be done efficiently using for example the Lanczos algorithm golub2013matrix. Let \[ \widehat E = \mathop{\bf diag} \left( \mathop{\bf diag}(R-\widehat R) \right), \] which can be shown to have positive diagonal entries. Our low-rank plus diagonal approximation of $R$ is then $\widehat R+\widehat E$. It is also a correlation matrix, {\it i.e.}, has diagonal entries one. Our final factor approximation of $\Sigma$ is given by \[ \mathop{\bf diag}(\sigma) (\widehat R+\widehat E) \mathop{\bf diag}(\sigma) = FF^T+E, \] with \[ F = \mathop{\bf diag}(\sigma) [\sqrt{\lambda_1} q_1 \cdots \sqrt{\lambda_k} q_k], \qquad E = \mathop{\bf diag}(e \circ \sigma^2), \] where $\circ$ denotes the elementwise (Hadamard) product, and $\sigma^2$ means elementwise.
This initialization alone can serve as a basic method to fit a factor model to a given covariance matrix. We will see below that in terms of portfolio optimization, it serves just as well as a factor model fit using the EM method.
\paragraph{Data set.} We gather the 500 largest NASDAQ stocks (by market capitalization) at the beginning of 2000 from the WRDS portal WRDS, compute the daily returns of these stocks from January 3rd 2000 to December 30th 2022, and remove any stocks with missing return values during this period. This gives us 238 stocks over 5787 (trading) days. We acknowledge that we induce a survivor bias, but the purpose of this empirical study is solely to demonstrate the benefit of regularization in large universes, and not to backtest a trading strategy.
\paragraph{Traditional factor model.} We create a factor model using PCA as follows. Every year, the $k$ principal components of largest explanatory power are computed, using the past two years of returns. These define the columns of the factor exposure matrix $F_t$ for the following year, and the factor returns $f_t$ are the projections of the returns onto these principal components. The idiosyncratic returns $z_t$ are the residuals. We leverage the CM-IEWMA predictor to compute the factor covariance, using three IEWMA components with half-lives (in days) $H^{\text{vol}}/H^{\text{cor}}$ of $\lceil k/2 \rceil/k, k/3k$, and $3k/6k$, where $k$ denotes the number of factors. To estimate the idiosyncratic variances a 21-day EWMA is used. We evaluate the factor models on the average log-likelihood over the evaluation period.
\paragraph{Fitting a factor model to the covariance matrix.} We use a CM-IEWMA covariance predictor with four IEWMA components with half-lives $63/125, 125/250, 250/500$, and $500/1000$ days, respectively, given as $H^{\text{vol}}/H^{\text{cor}}$. Given the CM-IEWMA estimate $\hat{\Sigma}_t$ at time $t$, we approximate it using a factor model as described in \S(ref).
To evaluate the factor models, we look at the average log-likelihood over the evaluation period and several performance metrics for a minimum variance portfolio with $L_\text{max}=1.6$, $w_\text{min} = -0.1$, and $w_\text{max}=0.15$, diluted to a target risk of 10%.
\paragraph{Traditional factor model.} Figure (ref) shows the log-likelihood versus the number of factors $k$ for $k$ between 2 and 75.
A large increase in log-likelihood is attained with around 20 factors, as compared to using the full covariance matrix. Thus using a traditional factor model and applying our covariance estimation method to the factor returns improves our overall covariance prediction.
\paragraph{Fitting a factor model to the covariance matrix.} Figure (ref) shows the log-likelihood versus the number of factors ({\it i.e.}, the rank of the low-rank component) $k$ for various $k$ between 2 and 75, using the eigendecomposition initialization and the EM algorithm.
We see that a rank of about $r=20$ seems optimal for this data set, and achieves a noticeably higher log-likelihood than using the full-rank covariance. Moreover, the EM algorithm does better than just computing the eigendecomposition.
Figure (ref) shows the portfolio metrics for the minimum variance portfolios. We can see that with roughly 10 factors or more, the performance is essentially identical to that obtained using the full covariance matrix. For these experiments we observed no notable difference between the two factor model fitting methods, {\it i.e.}, the simple eigendecomposition based initialization and the more sophisticated EM method. While using the factor model does not improve portfolio performance, it greatly speeds up the computation of the portfolio optimization problems.
\chapter{Smooth covariance predictions} We address here a secondary objective for a covariance prediction $\hat{\Sigma}_t$, which is that it vary smoothly across time. Perhaps the main reason for desiring smoothness of the covariance estimate is that it can lead to reduced trading in portfolio construction methods; it can also lead to improved portfolio performance, even without taking into account transaction costs.
To some extent smoothness happens naturally, since whatever method is used to form $\hat{\Sigma}_t$ from $r_1, \ldots, r_{t-1}$ is likely to yield a similar prediction $\hat{\Sigma}_{t+1}$ from $r_1, \ldots, r_t$. It is also possible to further smooth the predictions over time, perhaps trading off some performance, {\it e.g.}, in log-likelihood regret.
We have already mentioned that the weight optimization problem (ref) can be modified to encourage smoothness of the weights over time. We can also directly smooth the prediction $\hat{\Sigma}_t$, to get a smooth version $\hat \Sigma_t^\text{sm}$. A very simple approach is to let $\hat{\Sigma}_t^\text{sm}$ be a EWMA of $\hat{\Sigma}_t$, with a half-life chosen as a trade-off between smoother predictions and performance. This EWMA post-processing is equivalent to choosing $\hat\Sigma_t^\text{sm}$ to minimize \[ \left\|\hat\Sigma_t^\text{sm}-\hat{\Sigma}_t\right\|_F^2+ \lambda \left\| \hat\Sigma_t^\text{sm}- \hat\Sigma_{t-1}^\text{sm}\right\|_F^2, \] where $\lambda$ is a positive regularization parameter used to control the trade-off between smoothness and performance, or equivalently, the half-life of the EWMA post-processing. Here the first term is a loss, and the second is a regularizer that encourages smoothly varying covariance predictions.
We can create more sophisticated smoothing methods by changing the loss or the regularizer in this optimization formulation of smoothing. For example we can use the Kullback-Liebler (KL) divergence as a loss. With regularizer $\lambda \|\hat \Sigma_t^\text{sm} - \hat \Sigma_{t-1}^\text{sm}\|_F$ (no square in this case), we obtain a piecewise constant prediction, which roughly speaking only updates the prediction when needed. This is a convex optimization problem which can be solved quickly and reliably boyd2004convex.
We consider again the Fama-French factor returns from \S(ref), over the same time horizon. We use the CM-IEWMA covariance predictor with the same parameters as in \S(ref). \paragraph{Smoothly varying covariance.} In the first experiment we smooth the CM-IEWMA covariance estimates by applying a EWMA, which corresponds to the $\lambda \left\| \hat\Sigma_t^\text{sm}- \hat\Sigma_{t-1}^\text{sm}\right\|_F^2$ regularizer. For different EWMA half-lives we attain different levels of smoothness.
\paragraph{Piecewise constant covariance.} In the second experiment we smooth the CM-IEWMA covariance estimates by applying the $\lambda \left\| \hat\Sigma_t^\text{sm}- \hat\Sigma_{t-1}^\text{sm}\right\|_F$ regularizer. For different values of $\lambda$ we attain piecewise constant covariance predictors with different update frequencies.
\paragraph{Smoothly varying covariance.} Figure (ref) shows the regret versus smoothness for various levels of smoothness.
As seen, we can reduce the smoothness by a factor of (roughly) four without losing much performance in terms of regret. This can obviously be useful in practice since a smoother covariance estimate would, for example, reduce trading.
Table (ref) shows the portfolio metrics for various values of $\lambda$ for the minimum variance portfolio with the same parameters as in \S(ref); here the turnover is defined as the average of $252\times\|w_{t+1}-w_{t}\|_1/\|w_{t}\|_1$ over all times $t$ in the evaluation period.
Interestingly, the right amount of smoothing not only reduces turnover, but also improves portfolio performance in terms of Sharpe ratio and drawdown, while keeping the desired volatility level. Too much smoothing, however, leads to reduced portfolio performance. Figure (ref) shows the yearly annualized Sharpe ratios for $\lambda=10^{-4}$, indicating a stable performance over time.
Figure (ref) shows the portfolio weights for three different EWMA half-lives.
As seen, EWMA smoothing leads to smoothly varying portfolio weights, while the weights vary significantly when no smoothing is applied.
\paragraph{Piecewise constant covariance.} Figure (ref) shows the regret versus the update frequency of the covariance estimate.
There is a clear trade-off between the regret and update frequency. Roughly speaking, we could update the covariance matrix weekly without losing much in terms of regret.
As mentioned, a piecewise constant predictor can be desirable in practice, since it encourages us not updating the portfolio weights, which in turn reduces trading costs. Table (ref) shows the portfolio metrics for various values of $\lambda$ for the minimum variance portfolio from \S(ref).
As seen, smoothing can significantly reduce turnover, and interestingly improve the Sharpe ratio and drawdown noticeably while maintaining the correct risk level. Figure (ref) shows the yearly annualized Sharpe ratios for $\lambda=10^{-4}$. The performance is relatively stable over time, with a small downward trend.
To illustrate the impact of smoothing we show the portfolio weights for three different values of $\lambda$ in figure (ref).
Without smoothing the portfolio weights are updated significantly every day. For $\lambda=10^{-4}$ the weights are updated around once or twice a month. Finally, for $\lambda=10^{-5}$ the weights are updated on average every half a year, with only four big weight updates over the whole trading period. Interestingly, the weight updates for $\lambda=10^{-5}$ correspond precisely in time to the volatile regime around 1980, the 2000 dot-com bubble, the 2008 financial crisis, and the 2020 pandemic. In short, we can conclude from table (ref) and figure (ref) that smoothing can lead to less trading and improve the portfolio performance.
Finally, we note that there is some deviation between the regret metric and portfolio performance. As seen from figure (ref) regret increases as we update the covariance matrix less than every other week. However, as seen from table (ref) and figure (ref), portfolio performance can improve notably when updating the covariance matrix only every few months, or even years.
\chapter{Simulating returns}
Our model can be used to simulate future returns, when seeded by past realized ones. To do this, we start with realized returns for periods $1, \ldots, t-1$, and compute $\hat{\Sigma}_t$ using our method. Then we generate or sample $r_t^\text{sim}$ from $\mathcal N(0,\hat{\Sigma}_t)$. We then find $\hat{\Sigma}_{t+1}$ using the returns $r_1, \ldots, r_{t-1}, r^\text{sim}_t$. We generate $r_{t+1}^\text{sim}$ by sampling from $\mathcal N(0,\hat{\Sigma}_{t+1})$. This continues.
This simple method generates realistic return data in the short term. Of course, it does not include shocks or rapid changes in the return statistics that we would see in real data, but the generative return method has several practical applications. To mention just one, we can simulate 100 (say) different realizations over the next quarter (say), and use these to compute 100 performance metrics for our portfolio. This gives us a distribution of the performance metric that we might see over the next quarter.
To illustrate the generative return method, we consider the five Fama-French factor returns from \S(ref). Using the same setup as in \S(ref) we compute CM-IEMWA covariance estimates, using data from January 1st 2011 to December 31 2013, {\it i.e.}, over a three-year period. Returns are then generated for 100 days, using the generative mode described above.
We illustrate the results by looking at the SMB factor, {\it i.e.}, we look at the marginal distribution of this factor. Figure (ref) shows the true SMB factor returns and the simulated returns for two different random number generator seeds.
As seen, we attain realistic returns that could be used to generate scenarios for downstream portfolio optimization tasks, for example.
\chapter{Conclusions} We have introduced a simple method for predicting covariance matrices of financial returns. Our method combines well known ideas such as EWMA, first estimating volatilities and then correlations, and dynamically combining multiple predictions. The method relies on solving a small convex optimization problem (to find the weights used in the combining), which is extremely fast and reliable. The proposed predictor requires little or no tuning or fitting, is interpretable, and produces results better than the popular EWMA estimate, and comparable to MGARCH. Given its interpretability, light weight, and good practical performance, we see it as a practical choice for many applications that require predictions of the covariance of financial returns.
\backmatter
\printbibliography