Multi Anchor Point Shrinkage for the Sample Covariance Matrix (Extended Version)

Next we describe the mathematical setting, motivation, and results in more detail. We restrict attention to a familiar and well-studied baseline model for financial returns: the one-factor, or “market”, model

$$\mathbf{r}=\beta x+\mathbf{z},$$

(1)

where $\mathbf{r}\in\mathbb{R}^{p}$ is a $p$-dimensional random vector of asset (excess) returns in a universe of $p$ assets, $\beta\in\mathbb{R}^{p}$ is an unobserved non-zero vector of parameters to be estimated, $x$ is an unobserved random variable representing the common factor return, and $\mathbf{z}\in\mathbb{R}^{p}$ is an unobserved random vector of residual returns.

With the assumption that the components of $\mathbf{z}$ are uncorrelated with $x$ and each other, the returns of different assets are correlated only through $\beta$, and therefore the covariance matrix of $\mathbf{r}$ is

$$\Sigma=\sigma^{2}\beta\beta^{T}+\Delta,$$

where $\sigma^{2}$ denotes the variance of $x$, and $\Delta$ is the diagonal covariance matrix of $\mathbf{z}$.

Under the further simplifying model assumption¹¹1The assumption of homogeneous residual variance $\delta^{2}$ is a mathematical convenience. If the diagonal covariance matrix $\Delta$ of residual returns can be reasonably estimated, then the problem can be rescaled as $\Delta^{-1/2}\mathbf{r}=\Delta^{-1/2}\beta x+\Delta^{-1/2}\mathbf{z}$, which has covariance matrix $\sigma^{2}\beta_{\Delta}\beta_{\Delta}^{T}+I$, where $\beta_{\Delta}=\Delta^{-1/2}\beta$. that each component of $\mathbf{z}$ has a common variance $\delta^{2}$ (also not observed), we obtain the covariance matrix of returns

$$\Sigma=\sigma^{2}\beta\beta^{T}+\delta^{2}\mathbf{I},$$

(2)

where $\mathbf{I}$ denotes the $p\times p$ identity matrix.

This means that $\beta$, or its normalization $b=\beta/||\beta||$, is the leading eigenvector of $\Sigma$, corresponding to the largest eigenvalue $\sigma^{2}||\beta||^{2}+\delta^{2}$. As estimating $b$ becomes the most significant part of the estimation problem for $\Sigma$, a natural approach is to take as an estimate the first principal component (leading unit eigenvector) $h_{PCA}$ of the sample covariance of returns data generated by the model. This principal component analysis (PCA) estimate is our starting point.

Consider the optimization problem

	$\displaystyle\min_{w\in\mathbb{R}^{p}}w^{T}\Sigma w$
	$\displaystyle e^{T}w=1$

where $e=(1,1,\dots,1)$, the vector of all ones.

The solution, the “minimum variance portfolio”, is the unique fully invested portfolio minimizing the variance of returns. Of course the true covariance matrix $\Sigma$ is not observable and must be estimated from data. Denote an estimate by

$$\hat{\Sigma}=\hat{\sigma}^{2}\hat{\beta}\hat{\beta}^{T}+\hat{\delta}^{2}\mathbf{I}$$

(3)

corresponding to estimated parameters $\hat{\sigma}$, $\hat{\beta}$, and $\hat{\delta}$.

Let $\hat{w}$ denote the solution of the optimization problem

	$\displaystyle\min_{w\in\mathbb{R}^{p}}w^{T}\hat{\Sigma}w$
	$\displaystyle e^{T}w=1.$

It is interesting to compare the estimated minimum variance

$$\hat{V}^{2}=\hat{w}^{T}\hat{\Sigma}\hat{w}$$

with the actual variance of $\hat{w}$:

$$V^{2}=\hat{w}^{T}\Sigma\hat{w},$$

and consider the variance forecast ratio $V^{2}/\hat{V}^{2}$ as one measure of the error made in the estimation of minimum variance, hence of the covariance matrix $\Sigma$.

The remarkable fact proved in Goldberg et al., [2021] is that, asymptotically as $p$ tends to infinity, the true variance of the estimated portfolio doesn’t depend on $\hat{\sigma}$, $\hat{\delta}$, or $||\hat{\beta}||$, but only on the unit eigenvector $\hat{\beta}/||\hat{\beta}||$. Under some mild assumptions stated later, they show the following.

We use the notation for normalized vectors

$$b=\frac{\beta}{||\beta||},\text{ }q=\frac{e}{\sqrt{p}},\text{ and }h=\frac{\hat{\beta}}{||\hat{\beta}||}.$$

Goldberg, Papanicolaou and Shkolnik call the quantity $\mathcal{E}(h)$ the optimization bias associated to an estimate $h$ of the true vector $b$. They note that the optimization bias $\mathcal{E}(h_{PCA})$ is asymptotically bounded above zero almost surely, and hence the variance forecast ratio explodes as $p\to\infty$.

With this background, the estimation problem becomes focused on finding a better estimate $h$ of $b$ from an observed time series of returns. GPS Goldberg et al., [2021] introduces a shrinkage estimate for $b$ – the GPS estimator $h_{GPS}$ – obtained by “shrinking” the PCA eigenvector $h_{PCA}$ along the unit sphere toward $q$, to reduce excess dispersion. That is, $h_{GPS}$ is obtained by moving a specified distance (computed only from observed data) toward $q$ along the spherical geodesic connecting $h_{PCA}$ and $q$. “Shrinkage” refers to the reduced geodesic distance to the “shrinkage target” $q$.

The GPS estimator $h_{GPS}$ is a significant improvement on $h_{PCA}$. First, $\mathcal{E}(h_{GPS})$ tends to zero with $p$, and in fact $p\mathcal{E}^{2}(h_{GPS})/\log\log(p)$ is bounded (proved in Gurdogan, [2021]). In Goldberg et al., [2021] it is conjectured, with numerical support, that $E[p\mathcal{E}^{2}(h_{GPS})]$ is bounded in $p$, and hence the expected variance forecast ratio remains bounded. Moreover, asymptotically $h_{GPS}$ is closer than $h_{PCA}$ to the true value $b$ in the $\ell_{2}$ norm; and it yields a portfolio with better tracking error against the true minimum variance portfolio.

The purpose of this paper is to generalize the GPS estimator by introducing a way to use additional information about beta to adjust the shrinkage target $q$ in order to improve the estimate.

We can consider the space of all possible shrinkage targets $\tau$ as determined by the family of all nontrivial proper linear subspaces $L$ of $\mathbb{R}^{p}$ as follows. Given $L$ (assumed not orthogonal to $h$), let the unit vector $\tau(L)$ be the normalized projection of $h$ onto $L$. $\tau(L)$ is then a shrinkage target for $h$ determined by $L$ (and $h$). We will describe such a subspace $L$ as the linear span of a set of unit vectors called “anchor points”. In the case of a single anchor point $q$, note that $\tau(\text{span}\{q\})=q$, so this case corresponds to the GPS shrinkage target.

The “MAPS” estimator is a shrinkage estimator with a shrinkage target defined by an arbitrary collection of anchor points, usually including $q$. When $q$ is the only anchor point, the MAPS estimator reduces to the GPS estimator. We can therefore think of the MAPS approach as allowing for the incorporation of additional anchor points when this provides additional information.

In Theorem 2.2, we show that expanding $\text{span}\{q\}$ by adding additional anchor points at random asymptotically does no harm, but makes no improvement.

In Theorem 2.3, we show that if the user has certain mild a priori rank ordering information about groups of components of $\beta$, even with no information about magnitudes, an appropriately constructed MAPS estimator converges exactly to the true vector $b$ in the asymptotic limit.

Theorem 2.4 shows that if the betas have positive serial correlation over recent history, then adding the prior PCA estimator $h$ as an anchor point improves the $\ell_{2}$ error in comparison with the GPS estimator, even if the GPS estimator is computed with the same total data history.

The benefit of improving the $\ell_{2}$ error in addition to the optimization bias is that it also allows us to reduce the tracking error of the estimated minimum variance fully invested portfolio, discussed in Section 3 and Theorem 3.1.

In the next sections we present the main results. The framework, assumptions, and statements of the main theorems are presented in Sections 2 and 3. Some simulation experiments are presented in Section 4 to illustrate the impact of the main results for some specific situations.

Proofs of the theorems of Section 2 are organized in Section 5, with the proofs of some of the needed technical propositions and lemmas appearing in Section 6. Additional details and computations may be found in Gurdogan, [2021].

We consider a simple random sample history generated from the basic model (1). The sample data can be summarized as

$$R=\beta X^{T}+Z$$

(5)

where $R\in\mathbb{R}^{p\times n}$ holds the observed individual (excess) returns of $p$ assets for a time window that is set by $n$ consecutive observations. We may consider the observables $R$ to be generated by non-observable random variables $\beta\in\mathbb{R}^{p}$, $X\in\mathbb{R}^{n}$ and $Z\in\mathbb{R}^{p\times n}$.

The entries of $X$ are the market factor returns for each observation time; the entries of $Z$ are the specific returns for each asset at each time; the entries of $\beta$ are the exposure of each asset to the market factor, and we interpret $\beta$ as random but fixed at the start of the observation window of times $1,2,3,...,n$ and remaining constant throughout the window. Only $R$ is observable.

In this paper we are interested in asymptotic results as $p$ tends to infinity with $n$ fixed. Therefore we consider equation (5) as defining an infinite sequence of models, one for each $p$.

To specify the relationship between models with different values of $p$, we need a more precise notation. We’ll let $\beta$ refer to an infinite sequence $(\beta(1),\beta(2),\dots)\in\mathbb{R}^{\infty}$, and $\beta^{p}=(\beta(1),\dots,\beta(p))\in\mathbb{R}^{p}$ the vector obtained by truncation after $p$ entries. When the value $p$ is understood or implied, we will frequently drop the superscript and write $\beta$ for $\beta^{p}$.

Similarly, $Z\in\mathbb{R}^{\infty\times n}$ is a vector of $n$ sequences (the columns), and $Z^{p}\in\mathbb{R}^{p\times n}$ is obtained by truncating the sequences at $p$.

With this setup, passing from $p$ to $p+1$ amounts to simply adding an additional asset to the model without changing the existing $p$ assets. The $p$th model is denoted

$$R^{p}=\beta^{p}X^{T}+Z^{p},$$

but for convenience we will often drop the superscript $p$ in our notation when there is no ambiguity, in favor of equation (5).

Let $\mu_{p}(\beta)$ and $d_{p}(\beta)$ denote the mean and dispersion of $\beta^{p}$, given by

$$\mu_{p}(\beta)=\frac{1}{p}\sum\limits_{i=1}^{p}\beta(i)\text{ }\text{ and }\text{ }d_{p}(\beta)^{2}=\frac{1}{p}\sum\limits_{i=1}^{p}(\frac{\beta(i)-\mu_{p}(\beta)}{\mu_{p}(\beta)})^{2}.$$

(6)

We make the following assumptions regarding $\beta$, $X$ and $Z$:

1. A1.

   (Regularity of beta) The entries $\beta(i)$ of $\beta$ are uniformly bounded, independent random variables, fixed at time 1. The mean $\mu_{p}(\beta)$ and dispersion $d_{p}(\beta)$ converge to limits $\mu_{\infty}(\beta)\in(0,\infty)$ and $d_{\infty}(\beta)\in(0,\infty)$.

2. A2.

   (Independence of beta, X, Z) $\beta$, $X$ and $Z$ are jointly independent of each other.

3. A3.

   (Regularity of X) The entries $X_{i}$ of $X$ are iid random variables with mean zero, variance $\sigma^{2}$ .

4. A4.

   (Regularity of Z) The entries $Z_{ij}$ of $Z$ have mean zero, finite variance $\delta^{2}$, and uniformly bounded fourth moment. In addition, the $n$-dimensional rows of $Z$ are mutually independent, and within each row the entries are pairwise uncorrelated.²²2Note we do not assume $\beta,X$, or $Z$ are Normal or belong to any specific family of distributions.

We carry out our analysis with the projection of the vectors on the unit sphere $\mathbb{S}^{p-1}\subset\mathbb{R}^{p}$. To that end we define

$$b=\frac{\beta}{||\beta||}\text{ , }q=\frac{e}{\sqrt{p}},$$

(7)

where $e=e^{p}=(1,1,\dots,1)\in\mathbb{R}^{p}$, and $||.||$ denotes the usual Euclidean norm. With the given assumptions the covariance matrix $\Sigma_{\beta}$ of $R$, conditional on $\beta$, is

$$\Sigma_{\beta}=\sigma^{2}\beta\beta^{T}+\delta^{2}I.$$

(8)

Since $\beta$ stays constant over the $n$ observations, the sample covariance matrix $\frac{1}{n}RR^{T}$ converges to $\Sigma_{\beta}$ almost surely if $n$ is taken to $\infty$, and is the maximum likelihood estimator of $\Sigma_{\beta}$.

Since $b$ is a leading eigenvector of $\Sigma_{\beta}$ (corresponding to the largest eigenvalue), then the PCA estimator $h$ (the unit leading eigenvector $h$ of the sample covariance matrix $\frac{1}{n}RR^{T}$) is a natural estimator of $b$. (We always select the choice of unit eigenvector $h$ such that $(h,q)\geq 0$.)

Since $\beta$ and $X$ only appear in the model $R=\beta X+Z$ as a product, there is a scale ambiguity that we can resolve by combining their scales into a single parameter $eta$:

$$\eta^{p}=\frac{1}{p}|\beta^{p}|^{2}\sigma^{2}.$$

It is easy to verify that

$$\eta^{p}=\mu_{p}(\beta)^{2}(d_{p}(\beta)^{2}+1)\sigma^{2},$$

and therefore by our assumptions $\eta^{p}$ tends to a positive, finite limit $\eta^{\infty}$ as $p\to\infty$.

Our covariance matrix becomes

$$\Sigma_{\beta}\equiv\Sigma_{b}=p\eta bb^{T}+\delta^{2}I,$$

(9)

where we drop the superscript $p$ when convenient. The scalars $\eta,\delta$ and the unit vector $b$ are to be estimated by $\hat{\eta}$, $\hat{\delta}$, and $h$. As described above, asymptotically only the estimate $h$ of $b$ will be significant. Improving this estimate is the main technical goal of this paper.

In Goldberg et al., [2021] the PCA estimate $h$ is replaced by an estimate $h_{GPS}$ that is “data driven”, meaning that it is computable solely from the observed data $R$. We henceforth use the notation $h_{GPS}=\hat{h}_{q}$, for a reason that will be clear shortly. As an intermediate step we also consider a non-observable “oracle” version $h_{q}$, defined as the orthogonal projection in $\mathbb{S}^{p-1}$ of $b$ onto the geodesic joining $h$ to $q$. The oracle version is not data driven because it requires knowledge of the unobserved vector $b$ that we are trying to estimate, but it is a useful concept in the definition and analysis of the data driven version. Both the data driven estimate $\hat{h}_{q}$ and the oracle estimate $h_{q}$ can be thought of as obtained from the eigenvector $h$ via “shrinkage” along the geodesic connecting $h$ to the anchor point, $q$.

The GPS data-driven estimator $\hat{h}_{q}$ is successful in improving the variance forecast ratio, and in arriving at a better estimate of the true variance of the minimum variance portfolio. In this paper we have the additional goal of reducing the $l_{2}$ error of the estimator, which, for example, is helpful in reducing tracking error. To that end, we introduce the following new data driven estimator, denoted $\hat{h}_{L}$.

Let $L_{p}\subset\mathbb{R}^{p}$ denote a nontrivial proper linear subspace of $\mathbb{R}^{p}$. We will sometimes drop the dimension $p$ from the notation. Denote by $k_{p}$ the dimension of $L_{p}$, with $1\leq k_{p}\leq p-1$.

Let $h^{p}$ denote the normalized leading eigenvector of $\frac{1}{n}R^{p}(R^{p})^{T}$, $s_{p}^{2}$ its largest eigenvalue, and $l_{p}^{2}$ the average of the remaining eigenvalues. Then we define the data driven “MAPS” (Multi Anchor Point Shrinkage) estimator by

$$\hat{h}_{L}=\frac{\tau_{p}h+\underset{L}{\operatorname{proj}}(h)}{||\tau_{p}h+\underset{L}{\operatorname{proj}}(h)||}\text{ }\text{ where }\text{ }\tau_{p}=\frac{\psi_{p}^{2}-||\underset{L}{\operatorname{proj}}(h)||^{2}}{1-\psi_{p}^{2}}$$

(10)

and

$$\psi_{p}=\sqrt{\frac{s_{p}^{2}-l_{p}^{2}}{s_{p}^{2}}}$$

(11)

is the relative gap between $s_{p}^{2}$ and $l_{p}^{2}$.

When $L$ is the one-dimensional subspace spanned by the vector $q$, then $\hat{h}_{L}$ is precisely the GPS estimator $\hat{h}_{q}$, located along the spherical geodesic connecting $h$ to $q$. The phrase “multi anchor point” comes from thinking of $q$ as an “anchor point” shrinkage target in the GPS paper, and $L$ as a subspace spanned one or more anchor points. The new shrinkage target determined by $L$ is the normalized orthogonal projection of $h$ onto $L$. When $L$ is the one-dimensional subspace spanned by $q$, the normalized projection of $h$ onto $L$ is just $q$ itself. In the event that $L$ is orthogonal to $h$, the MAPS estimator $\hat{h}_{L}$ reverts to $h$ itself.

Does adding anchor points to create a MAPS estimator from a higher-dimensional subspace improve the estimation? The answer depends on whether there is any relevant information in the added anchor points.

We need the concept of a random linear subspace of $\mathbb{R}^{p}$. Let $k_{p}$ be a positive integer such that $1\leq k_{p}\leq p-1$. Let $\xi^{p}$ be an $O(p)$-valued random variable, where $O(p)$ denotes the orthogonal group in $\mathbb{R}^{p}$. Let $\{e^{p}_{1},e^{p}_{2},\dots,e^{p}_{p}\}$ denote the standard Cartesian basis of $\mathbb{R}^{p}$.

We say that $L_{p}$ is a random linear subspace of $\mathbb{R}^{p}$ with dimension $k_{p}$ if, for some $\xi^{p}$ as above,

$$L_{p}=\text{span}_{p}\{\xi^{p}e^{p}_{i}|i=1,2,\dots,k_{p}\},$$

where $\text{span}_{p}$ denotes the linear span of a set of vectors in $\mathbb{R}^{p}$.

We say $L_{p}$ is independent of a random variable $X$ if the generator $\xi^{p}$ is independent of $X$. Moreover, we say $H_{p}$ is a Haar random subspace of $\mathbb{R}^{p}$ if it is a random linear subspace as above, and the random variable $\xi^{p}$ induces the (uniform) Haar measure on $O(p)$.

For example, any non-decreasing sequence satisfying $k_{p}\leq Cp^{\alpha}$ for $\alpha<1/2$ is square root dominated.

Theorem 2.2 says adding random anchor points to form a MAPS estimator does no harm, but also makes no improvement asymptotically. Equation (13) says that the GPS estimator is neither improved nor harmed by adding extra anchor points uniformly at random. Therefore the goal will be to find useful anchor points that take advantage of additional information that might be available.

As with stocks grouped by sector, it may be that the betas can be separated into ordered groups, where the rank ordering of the groups is known, but not the ordering within groups. This turns out to be enough information for the MAPS estimator to converge asymptotically to the true value almost surely.

Theorem 2.3 says that if we have certain prior information about the ordering of the $\beta$ elements in the sense of finding an ordered partition (but with no prior information about the magnitudes of the elements or their ordering within partition atoms), then asymptotically we can estimate $b$ exactly.

Having in hand a genuine ordered partition a priori is likely only approximately possible in the real world. Theorem 2.3 is suggestive that even partial grouped order information about the betas can be helpful in strictly improving the GPS estimate. This is confirmed empirically in section 4.

The next theorem shows that even with no a priori information beyond the observed time series of returns, we can still use MAPS to improve the GPS estimator.

In the analysis above we have treated $\beta$ as a constant throughout the sampling period, but in reality we expect $\beta$ to vary slowly over time. To capture this in a simple way, let’s now assume that we have access to returns observations for $p$ assets over a fixed number of $2n$ periods. The first $n$ periods we call the first (or previous) time block, and the second $n$ periods the second (or current) time block. We then have returns matrices $R_{1},R_{2}\in\mathbb{R}^{p\times n}$ corresponding to the two time blocks, and $R=[R_{1}R_{2}]\in\mathbb{R}^{p\times 2n}$ the full returns matrix over the full set of $2n$ observation times.

Define the sample covariance matrices $S,S_{1},S_{2}$ as $\frac{1}{2n}RR^{T}$, $\frac{1}{n}R_{1}R_{1}^{T}$, and $\frac{1}{n}R_{1}R_{1}^{T}$, respectively. Let $h,h_{1},h_{2}$ denote the respective (normalized) leading eigenvectors (PCA estimators) of $S,S_{1},S_{2}$. (Of the two choices of eigenvector, we always select the one having non-negative inner product with $q$.)

Instead of a single $\beta$ for the entire observation period, we suppose there are random vectors $\beta_{1}$ and $\beta_{2}$ that enter the model during the first and second time blocks, respectively, and are fixed during their respective blocks. We assume both $\beta_{1}$ and $\beta_{2}$ satisfy assumptions (1) and (2) above, and denote by $b_{1}$ and $b_{2}$ the corresponding normalized vectors. The vectors $\beta_{1}$ and $\beta_{2}$ should not be too dissimilar in the mild sense that $(\beta_{1},\beta_{2})\geq 0$.

The Cauchy-Schwartz inequality shows $-1\leq\rho_{p}(\beta_{1},\beta_{2})\leq 1$. Furthermore, it is straightforward to verify that

$$(b_{1},b_{2})-(b_{1},q)(b_{2},q)=\frac{d_{p}(\beta_{1},\beta_{2})}{\sqrt{1+d_{p}(\beta_{1})^{2}}\sqrt{1+d_{p}(\beta_{2})^{2}}}.$$

(19)

and hence $d_{p}(\beta_{1},\beta_{2})$, and $\rho_{p}(\beta_{1},\beta_{2})$ have limits $d_{\infty}(\beta_{1},\beta_{2})$, and $\rho_{\infty}(\beta_{1},\beta_{2})$ as $p\to\infty$.

Our motivation for this model is the intuition that the betas for different time periods are noisy representations of a fundamental beta, and that the beta from a recent time block provides some useful information about beta in the current time block. To make this precise in support of the following theorem, we make the following additional assumptions.

Theorem 2.4 says that the MAPS estimator obtained by adding the PCA estimator $h$ from the previous time block as a second anchor point outperforms the GPS estimator asymptotically, as measured by $l_{2}$ error, even if the latter estimated with the full $2n$ (double) data set. This works when the previous time block carries some information about the current beta (non-zero correlation). In the case of perfect correlation $\rho_{\infty}(\beta_{1},\beta_{2})=1$ the two betas are equal, and we then return to the GPS setting where beta is assumed constant, so no improved performance is expected.

The cost of implementing this “dynamic MAPS” estimator is comparable to that of the GPS estimator, so should generally be preferred when no rank order information is available for beta.