Lab | Shrinkage

24.6 Shrinkage

Key points

Shrinkage estimators aim to reduce the variance of an estimator, and hence to improve its efficiency.
The most popular shrinkage estimator for the expectation is the James-Stein estimator, which averages the sample mean with a constant estimator (24.71).
Several covariance shrinkage estimators are available (Ledoit-Wolf (24.84), graphical lasso (24.102), factor analysis (24.105), etc.).

Shrinkage estimators average a high-variance estimate, say the sample mean, with a high-bias, subjective estimate of the value of the unknown parameter. The resulting shrinkage estimators trade off bias for efficiency, resulting in better estimators, according to the estimation assessment criteria that we discuss in Section 27.1.

In theory, we can apply shrinkage to the estimation of any property of a distribution. In practice, shrinkage is used mostly with expectations, covariances/correlations, and with factor loadings in linear factor models (14.14).

We cover the classical James-Stein shrinkage of the expectation parameters in Section 24.6.1.

In the next sections we present different techniques to construct shrinkage estimators of the covariance or correlation matrix: Ledoit-Wolf shrinkage ( Section 24.6.2); random matrix theory (Section 24.6.3); sparse eigenvectors rotations (Section 24.6.4); glasso over Markov networks (Section 24.6.5); factor analysis ( Section 24.6.6).

Another shrinkage methodology for the correlation matrix is the homogenous clusters correlation shrinkage, which we cover in the context of machine learning techniques, as an application of k-means clustering. Refer to Section 70.

We discuss the shrinkage of factor loadings in Section 25.5, fully dedicated to linear factor models.

Shrinkage is also a feature of Bayesian estimation, whereby data-driven estimates are naturally shrunk towards the Bayesian prior. We discuss Bayesian estimation for expectation and covariances in Section 24.5; and for factor loadings in Section 25.3

In Section 33.6 we discuss generalized shrinkage in the more general context of views processing.

24.6.1 Mean shrinkage: James-Stein

In this section we address the shrinkage estimation of the expectation (8.13) of a set of i.i.d. variables.

Historically, the first shrinkage estimator was the James-Stein estimator of the expectation, refer to the original article [Stein, 1955] and see also [Lehmann and Casella, 1998].

To understand the James-Stein shrinkage estimator, let us start from the sample mean $ˆ μ$ (24.116) of a time series of i.i.d. variables ${ϵ_{1}, \dots, ϵ_{¯ t}}$ (36.2).

The sample mean is unbiased but highly inefficient, see Example 27.10. In words, when we use as input different time series ${ϵ_{1}, \dots, ϵ_{¯ t}}$ , the output estimate varies widely, see Figure 23.15. As a result, “optimal” mean-variance allocations based on the sample mean vary wildly when different time series are used as input into the allocation process, see Section 60b.1.3.

Next, let us consider an estimator that is extremely efficient, namely a target vector $μ_{t a r g e t}$ which completely disregards any historical input data. A choice can be $μ_{t a r g e t} = 0$ as in [Alexander, 1998]. Alternatively, forward-looking targets implied by the notion of carry, see Section 55.1.2.

Then, as it turns out, we can improve the sample mean by “shrinking” it towards the shrinkage target: this way we obtain the James-Stein estimator

¯ ¯ ¯ μ \equiv (1 - γ) ˆ μ + γ μ_{t a r g e t},

(24.71)

where $γ \in (0, 1)$ is a suitable constant that represents the confidence in our shrinkage target. The limit case $γ = 0$ implies no shrinkage (full confidence in the sample mean), whereas the limit case $γ = 1$ represents full confidence in the shrinkage target, as we disregard the sample mean.

For normally distributed variables $ε_{t} \sim N (μ, σ^{2})$ a sub-optimal confidence reads 70.29

γ = min {1, \frac{2}{¯ t} \frac{tr (σ^{2}) - 2 λ_{1}}{{(ˆ μ - μ_{t a r g e t})}_{t a r g e t}^{'} (ˆ μ - μ_{t a r g e t})}},

(24.72)

where $λ_{1}$ denotes the largest eigenvalue (43.412) of $σ^{2}$ . Notice how the confidence in the target decreases as the number of available observations increases: $lim ¯ t \to \infty γ = 0$ . Also note that the sub-optimal confidence (24.72) only applies when the dimension satisfies $¯ ı ≫ 2$ [W], or else the shrinkage factor $γ$ becomes negative. In practice, we can replace the unknown population parameters $(σ^{2}, λ_{1})$ by their sample counterparts $({ˆ σ}^{2}, {ˆ λ}_{1})$ (24.3).

The shrinkage of the mean toward a prior target (24.71) also follows from a purely Bayesian approach to estimation, see Example 24.21.

Example 24.21. In Example 24.13 on Bayesian theory, we consider a (univariate) normal i.i.d. variable $ε_{t}$ (24.12), we compute the conditional likelihood $f (i | μ, σ^{2})$ (24.40) and we assign a normal-inverse-Wishart model (24.45)-(24.47) to the prior distribution of the parameters $(M, Σ^{2})$ .
In Example 24.15 we then use the Bayesian approach to obtain the posterior distribution (23.172) of $(M, Σ^{2})$ , which turns out to be normal-inverse-Wishart itself (24.53)-(24.54). In this context, notice that the posterior location parameter $μ^{p o s}$ (24.56) mirrors the structure of the James-Stein shrinkage estimator (24.71), where we take as target the prior parameter $μ^{p r i} \equiv μ_{t a r g e t}$ and we set the confidence as $γ \equiv \frac{t^{p r i}}{t^{p r i} + ¯ t}$ : also note that, like in the optimal confidence expression (24.72), the confidence $γ$ in the target $μ^{p r i}$ decreases as the number of observations $¯ t$ increases (24.66).

Example 24.22. Sample mean shrinkage: James-Stein estimator

We consider the time series of log-returns $ϵ_{t}$ of the stocks in the S&P 500 and we compute their sample mean $ˆ μ$ (24.116). We partition the S&P into the $¯ k = 10$ MSCI US sector indices $s_{1} \equiv$ “Energy”, $s_{2} \equiv$ “Materials”, $s_{3} \equiv$ “Industrials”, $s_{4} \equiv$ “Consumer Discretionary”, $s_{5} \equiv$ “Consumer Staple”, $s_{6} \equiv$ “Healthcare”, $s_{7} \equiv$ “Financials”, $s_{8} \equiv$ “IT”, $s_{9} \equiv$ “Telecommunications”, $s_{10} \equiv$ “Utilities”.
As a simple example of target mean $μ_{t a r g e t}$ in (24.71) we can use the grand mean, which is the cross-sectional average of the entries of the sample mean

{[μ_{t a r g e t}]}_{i} = \frac{\sum_{j \in s_{k}} {ˆ μ}_{j}}{\sum_{j \in s_{k}} 1}, if i \in s_{k} .

(24.73)

Then we take $γ = 1$ in formula (24.71) so that the shrinkage mean coincides with the target mean, $¯ ¯ ¯ μ \equiv μ_{t a r g e t}$ .
Figure 24.14 shows the sample mean $ˆ μ$ and the shrinkage mean $¯ ¯ ¯ μ$ of the returns of two stocks in the same sector. Both means are computed over 100 different random small sub-samples of the whole sample. From the whole sample we compute the global mean, which represents the true “hidden” return, in the spirit of cross-validation, see Section 30.5. We observe that the sample mean varies widely across different sub-samples, i.e. it is highly inefficient. More formally, the average loss (27.47), or risk (27.46), for the sample mean is larger than for the shrinkage estimator. Figure 24.14 also shows the mean-variance optimal weights using shrinkage mean and sample mean. Notice that, using the sample mean, the allocation varies widely.

24.6.2 Covariance shrinkage: Ledoit-Wolf

The Ledoit-Wolf shrinkage in [Ledoit and Wolf, 2004] follows from the observation that the eigenvalues of the estimated covariance or correlation matrix tend to be more dispersed than the eigenvalues of the true data generating process.

To illustrate this effect, let us consider the case where the data generating process $ε_{t}$ is standardized noise (36.2), i.e. all its entries are independent and identically distributed, with zero expectation and unitary variance

ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'} \sim f_{ε} i.i.d., with ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} E {ε_{i, t}} = 0 V {ε_{i, t}} = 1 C v {ε_{i, t}, ε_{j, t}} = 0, for i \neq j \end{matrix} .

(24.74)

We do not need to assume that the entries $ε_{i, t}$ are normally distributed. However, we can consider normally distributed (49.1) noise as an example

ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'} \sim N (0_{¯ ı \times 1}, I_{¯ ı}) i.i.d..

(24.75)

In this case, if we consider the spectral decomposition (43.372) of the true covariance

{λ, e} = e i g (C v {ε_{t}}),

(24.76)

we obtain that all eigenvalues are equal to $1$

C v {ε_{t}} = I_{¯ ı} \Rightarrow λ = 1_{¯ ı} .

(24.77)

Suppose that we have a time series of realizations $i = {ϵ_{1}, \dots, ϵ_{¯ t}}$ (23.1) of the standardized noise (24.74). Let us now consider the simplest estimator of the covariance, namely the sample covariance ${ˆ σ}_{ε}^{2}$ (24.117), and its spectral decomposition (43.372)

{ˆ λ, ˆ e} = e i g ({ˆ σ}_{ε}^{2}) .

(24.78)

Empirically, the sample spectrum $ˆ λ$ (24.78) is more dispersed, or steeper, than the true spectrum $λ$ (24.76). To quantify this statement, let us define the condition number [W] of a dispersion matrix as the ratio of the largest to the smallest eigenvalue

κ (σ^{2}) \equiv \frac{λ_{1}}{λ_{¯ ı}} .

(24.79)

Then the sample condition number $κ ({ˆ σ}^{2})$ is larger than the true condition number

κ ({ˆ σ}^{2}) \equiv \frac{{ˆ λ}_{1}}{{ˆ λ}_{¯ ı}} > κ (σ^{2}) \equiv \frac{λ_{1}}{λ_{¯ ı}} = 1 .

(24.80)

In particular, let us consider the extreme, although in financial applications not rare, case where the number of observations $¯ t$ is less than the dimension $¯ ı$ of the variables $ε_{t}$ . In this case, the sample covariance ${ˆ σ}^{2}$ becomes singular 70.24

¯ t < ¯ ı \Rightarrow r a n k ({ˆ σ}^{2}) = ¯ t < ¯ ı,

(24.81)

which implies that at least the last eigenvalue is zero ${ˆ λ}_{¯ ı} = 0$ . Hence the sample condition number goes to infinity, i.e. ${ˆ σ}^{2}$ is ill-conditioned

¯ t < ¯ ı \Rightarrow κ ({ˆ σ}^{2}) = + \infty .

(24.82)

The degeneration of the condition number bears a geometrical interpretation, see Figure 24.15. The sample eigenvalues are more dispersed than the true eigenvalues, or equivalently the sample spectrum is steeper than the true spectrum. Hence, because the reference axes of a generic ellipsoid $\partial E (μ, σ^{2})$ have length equal to the square root of the eigenvalues of $σ^{2}$ (8.54), the ellipsoid $\partial E (0, σ^{2})$ (12.42) of the true covariance $σ^{2} \equiv C v {ε_{t}}$ is more spherical than the ellipsoid $\partial E (0, {ˆ σ}^{2})$ (12.42) of the sample covariance ${ˆ σ}^{2}$ (24.78), or intuitively

{λ_{i}}_{i = 1}^{¯ ı} \Leftrightarrow sphere estimation \to {{ˆ λ}_{i}}_{i = 1}^{¯ ı} \Leftrightarrow cigar.

(24.83)

Example 24.23. Estimated spectrum vs. true spectrum

In Figure 24.15 we illustrate the average estimation error on the spectrum.

We consider a multivariate standard normal white noise $ε_{t}$ (24.75) of large dimensions.

In the bottom right subplot we display the input, namely an $¯ ı \times ¯ t$ panel $i = (ϵ_{1} | \dots | ϵ_{¯ t})$ of joint i.i.d. realizations of the multivariate standard normal white noise (24.75). We draw one such panel $i^{(j)}$ multiple times for $j = 1, 2, \dots, ¯ ȷ$ . For any given panel, we compute the sample covariance matrix ${ˆ σ}^{2 (j)}$ and extract the respective sample eigenvalues ${{ˆ λ}_{i}^{(j)}}_{i}$ (24.78).

In the bottom left subplot we show the $¯ ı \times ¯ ı$ sample covariance matrix ${ˆ σ}^{2 (j)}$ (24.117) estimated from any given panel $j$ .

In the top left plot we show the true eigenvalues ${λ_{i}^{(j)}}_{i}$ (24.77), sorted in decreasing order, as a function of the entry number $i$ . Besides true spectrum, we also show the sample spectrum ${{ˆ λ}_{i}^{(j)}}_{i}$ 24.78) sorted in decreasing order,and its average across different simulations $j$ . The true spectrum is flat because all the entries are ones. However, the sample spectrum tilts, consistently with the Ledoit-Wolf effect.

In the top right plot we show the mean-covariance ellipsoid (8.47) of the first and last principal factors (8.82). The true principal factors are white noise, and therefore their mean-covariance ellipsoid is a circle (24.83). The mean-covariance ellipsoid of the sample principal factors is squeezed (24.83), because the largest principal axis is the square root of the largest sample eigenvalue, and similar for the smallest axis and sample eigenvalue, which are distorted by the Ledoit-Wolf effect.

As we change the input panel of realized white noise, the true covariance remains unaffected, but the sample covariance changes. However, as long as we do not change the number of joint realizations $¯ t$ in the panel, the average sample spectrum and the average covariance ellipsoid remain unaffected. When we decrease the number of joint realizations $¯ t$ , the average sample spectrum steepens and the average ellipsoid squeezes. In particular, when the number of joint realizations $¯ t$ becomes less than the dimension of the white noise, we obtain the first zero sample eigenvalue: the sample covariance matrix is no longer invertible (24.81). Then, as we lower further the number of joint realizations $¯ t$ , we obtain more zero eigenvalues, which is reflected in the right end of the average sample spectrum matching the horizontal axis, and the sample covariance ellipsoid flattening to a degenerate line.

The analysis is robust to distributional assumptions: if we perform a similar analysis with exponentially distributed white noise, we observe the precise same phenomenon, both qualitatively and quantitatively, in both this sample spectrum and the sample covariance ellipsoid.

A shrinkage estimator that addresses the steepening of the spectrum (24.83) for an arbitrary covariance matrix is developed in [Ledoit and Wolf, 2004]. Similarly to the James-Stein shrinkage estimator of the mean (24.71), the shrinkage estimator of the covariance is a convex combination [W] of the sample covariance ${ˆ σ}^{2}$ and a target covariance, set as a multiple of the identity matrix $σ_{t a r g e t}^{2} \equiv \frac{1}{¯ ı} tr ({ˆ σ}^{2}) I_{¯ ı}$

{¯ ¯ ¯ σ}^{2} = (1 - γ) {ˆ σ}^{2} + γ σ_{t a r g e t}^{2} .

(24.84)

In the case of normally distributed i.i.d. variables $ε_{t} \sim N (μ, σ^{2})$ , the optimal coefficient $γ \in [0, 1]$ can be derived by minimizing the error $∥ {¯ ¯ ¯ σ}^{2} - σ^{2} ∥_{F}^{2}$ , where $∥ \cdot ∥_{F}$ is the Frobenius norm (43.256), as explained in [Ledoit and Wolf, 2004] 70.30 . However it requires the knowledge of the true population covariance $σ^{2}$ , and therefore cannot be computed. [Ledoit and Wolf, 2004] propose an estimator for the theoretical optimal $γ$

ˆ γ = min {1, \frac{\frac{1}{{¯ t}^{2}} \sum_{t = 1}^{¯ t} tr {{(ϵ_{t} ϵ_{t}^{'} - {ˆ σ}^{2})}_{t}^{2}}}{tr {{({ˆ σ}^{2} - \frac{1}{¯ ı} tr ({ˆ σ}^{2}) I_{¯ ı})}_{¯ ı}^{2}}}},

(24.85)

which asymptotically converges to the true coefficient. Notice how $ˆ γ$ decreases as the number of observations available in the sample mean increases: $lim ¯ t \to \infty ˆ γ = 0$ .

The shrinkage estimator for the covariance (24.84) causes the spectrum to flatten. Hence, the shrinkage estimator is better conditioned [W] than the sample covariance 70.31 . The intuition behind the covariance shrinkage is as follows: the covariance shrinkage estimator (24.84) attempts to reverse the squeezing of the sphere caused by the estimation process (24.83).

24.6.3 Correlation shrinkage: random matrix theory

In Section 24.6.2 we used shrinkage to provide a more accurate spectral decomposition (43.413) of the true, unknown covariance matrix $σ^{2} \equiv C v {ε_{t}}$ of a set of i.i.d. variables $ε_{t}$ . More precisely, we flattened the spuriously steep spectrum of the sample estimate of the covariance matrix by averaging the sample estimate with a multiple of the identity matrix (24.84). A different approach to flatten spuriously steep spectrum of an estimate of the correlation matrix rests on random matrix theory, as discussed in Section 26.

Following [Rao, 2004], consider the Marchenko-Pastur distribution (26.300)

f_{Λ}^{M P} (x) \equiv (1 - \frac{1}{q}) δ^{(0)} (x) 1_{q > 1} + \frac{\sqrt{(x - λ^{-}) (λ^{+} - x)}}{2 π x q} 1_{x \in (λ^{-}, λ^{+})}^{-},

(24.86)

where we set

q \equiv \frac{¯ ı}{¯ t}, λ^{\pm} \equiv {(1 \pm \sqrt{q})}^{2}, σ^{2} \equiv 1 .

(24.87)

Consider now a standardized noise $ε_{t}$ that satisfies the assumptions (24.74). Given the vector $ˆ λ = {({ˆ λ}_{1}, \dots, {ˆ λ}_{¯ ı})}_{1}^{'}$ of sample eigenvalues, we can build its empirical pdf (26.36)

{ˆ f}_{Λ} (x | ˆ λ) \equiv \frac{1}{¯ ı} \sum_{i = 1}^{¯ ı} δ_{i}^{({ˆ λ}_{i})} (x) .

(24.88)

If the dimensions $¯ ı$ and $¯ t$ are large enough, then the empirical pdf (24.88) is approximated by the Marchenko-Pastur distribution (26.315)-(24.87)

{ˆ f}_{Λ} (x | ˆ λ) ¯ ı \equiv q ¯ t \to \infty \approx f_{Λ}^{M P} (x),

(24.89)

see Example 24.24.

The approximation (24.89) holds regardless of the distribution of the white noise $ε_{t}$ . Note that, when $q \leq 1$ , the (closed) support [W] of the density $f_{Λ}^{M P}$ (24.86) is the interval $[λ^{-}, λ^{+}]$ ; when $q > 1$ , the density $f_{Λ}^{M P}$ (24.86) has the same support, plus a spike at zero. Therefore, according to the density $f_{Λ}^{M P}$ , for any value of $q$ the estimated sample eigenvalues $ˆ λ$ are scattered away from the true eigenvalues $λ$ , which by construction are concentrated at $1$ .

It is possible to generalize the Marchenko-Pastur approximation (24.89) to the case of pure noise with general variance through its scaling property 70.22

{ˆ f}_{Λ} (x | ˆ λ) ¯ ı \equiv q ¯ t \to \infty \approx \frac{1}{λ} f_{Λ}^{M P} (\frac{x}{λ}),

(24.90)

and also to the non-zero correlation case [Bun et al., 2017]. In the context of random matrix theory another notable limit distribution is Wigner’s semicircle distribution (26.210),as discussed in Section 26.5.

Example 24.24. Correlation shrinkage: spectrum estimation via random matrix theory

In Figure 24.16 we illustrate the power of the Marchenko-Pastur distribution $f_{Λ}^{M P}$ (24.86).

Continuing from Example 24.23, in the top left plot we show the distribution of the eigenvalues, as approximated by a normalized histogram (49.281). The true distribution is peaked at $1$ with no other eigenvalues. The sample distribution is dispersed around the peak of the true spectrum according to the Marchenko-Pastur distribution (24.89).

When we change the number of joint realizations $¯ t$ , the Marchenko-Pastur distribution of the eigenvalues changes accordingly. In particular, when the number of joint realizations $¯ t$ becomes less than the dimension of the white noise $¯ ı$ , and we obtain more and more zero eigenvalues, the distribution of the eigenvalues peaks at zero.

This result is robust to distributional assumptions: when we perform a similar analysis with exponentially distributed white noise the Marchenko-Pastur distribution does not change.

Following [Bouchaud and Potters, 2009a] and [Bun et al., 2017], we can leverage random matrix theory for shrinkage purposes, as we proceed to explain. Consider a general data generating process $ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'}$ which generate data that are jointly i.i.d. across time and are a linear combinations of

- $¯ k < ¯ ı$ signals ${ε_{k}^{s i g n l}}_{k = 1}^{¯ k}$ independent of each other, with variances sorted in decreasing order (without loss of generality)

λ_{1} \geq \dots \geq λ_{¯ k};

(24.91)

- $¯ ı - ¯ k$ sources of noise ${ε_{k}^{s i g n l}}_{k = ¯ k + 1}^{¯ ı}$ independent of each and of the signals with smaller, equal variances $λ_{¯ k + 1} = \dots = λ_{¯ ı} = λ < λ_{¯ k}$ .

We can write the i.i.d. process $ε_{t}$ as follows

⎛ ⎜ ⎝ \begin{matrix} ε_{1} \cdot ε_{¯ ı} \end{matrix} ⎞ ⎟ ⎠ = (e_{1}, \dots, e_{¯ k}) \times ⎛ ⎜ ⎜ ⎝ \begin{matrix} ε_{1}^{s i g n l} \cdot ε_{¯ k}^{s i g n l} \end{matrix} ⎞ ⎟ ⎟ ⎠ + (e_{¯ k + 1}, \dots, e_{¯ ı}) \times o \times ⎛ ⎜ ⎜ ⎝ \begin{matrix} ε_{¯ k + 1}^{n o i s e} \cdot ε_{¯ ı}^{n o i s e} \end{matrix} ⎞ ⎟ ⎟ ⎠,

(24.92)

where

- $e \equiv (e_{1}, \dots, e_{¯ ı})$ is an $¯ ı \times ¯ ı$ rotation/reflection (43.231) in $R^{¯ ı}$ and thus it satisfies $e e^{'} = I_{¯ ı}$ ;

- $o$ an $(¯ ı - ¯ k) \times (¯ ı - ¯ k)$ rotation/reflection (43.231) in $R^{¯ ı - ¯ k}$ and thus it satisfies $o o^{'} = I_{(¯ ı - ¯ k)}$ .

The true spectrum (24.76) of the covariance matrix $C v {ε_{t}}$ of the i.i.d. process (24.92) coincides with the juxtaposition of the variances of both signals and noise, i.e.

e i g v a l (C v {ε_{t}}) = {(λ_{1}, \dots, λ_{¯ k}, λ, \dots, λ)}_{1}^{'},

(24.93)

while the eigenvectors (43.369) turn out to be the columns of matrix $e$ 70.25 .

Now consider a time series for the i.i.d. process $i = {ϵ_{1}, \dots, ϵ_{¯ t}}$ . For the signals, the sample eigenvectors and the sample eigenvalues are close to their true counterparts ${ˆ λ}_{k} \approx λ_{k}$ and ${ˆ e}_{k} \approx e_{k}$ for $k = 1, \dots, ¯ k$ . Instead, the sample eigenvalues corresponding to the many sources of noise ( $k = ¯ ı - ¯ k + 1, \dots, ¯ ı$ ) will be affected by the elliptical distortion (24.83) illustrated in Figure 24.15.

Then the empirical pdf of the sample eigenvalues (24.88) splits into two components: the distribution of the large number of small eigenvalues $({ˆ λ}_{¯ k + 1}, \dots, {ˆ λ}_{¯ ı})$ stemming from the noise, which follows the rescaled Marchenko-Pastur law (24.90); and the spikes $δ_{k}^{({ˆ λ}_{k})} (x)$ corresponding to the large eigenvalues $({ˆ λ}_{1}, \dots, {ˆ λ}_{¯ k})$ stemming from the signals

{ˆ f}_{Λ} (x | ˆ λ) \approx w \times \frac{1}{¯ k} \sum_{k \leq ¯ k} δ_{k}^{({ˆ λ}_{k})} (x) + (1 - w) \times \frac{1}{λ_{n o i s e}} f_{Λ}^{M P} (\frac{x}{λ_{n o i s e}}),

(24.94)

where

- $δ^{(λ)} (x)$ denotes the Dirac delta (46.95) centered in $λ$ , which spikes in $x = λ$ ;

- $w = ¯ k ∕ ¯ ı$ is the relative weight of the two distributions: $¯ k$ signal/spike eigenvalues, versus $¯ ı - ¯ k$ Marchenko-Pastur/noise eigenvalues over a total of $¯ ı$ eigenvalues;

- $λ_{n o i s e}$ is the average of the noise eigenvalues $λ_{n o i s e} \equiv \frac{1}{¯ ı - ¯ k} \sum_{k > ¯ k} {ˆ λ}_{k}$ ; this number is the counterpart of $λ$ in the true spectrum (24.93);

- $f_{Λ}^{M P} (x)$ is the Marchenko-Pastur distribution (24.86), where $q = (¯ ı - ¯ k) ∕ ¯ t$ is the ratio of number of eigenvalues over observations, adapted from the ratio $q \equiv \frac{¯ ı}{¯ t}$ (24.87) to the current situation of $¯ ı - ¯ k$ noise eigenvalues.

At this point we have all the theoretical instruments required to understand and perform spectrum shrinkage in practice. In particular, the decomposition (24.94) can be used to model the histogram of the sample spectrum $ˆ λ$ stemming from any input dataset of observations $i = {ϵ_{1}, \dots, ϵ_{¯ t}}$ .

First of all, let us assume the observations $ϵ_{t}$ to be already standardized, as explained in the joint time series (53b.56), in order to obtain correlation-like covariances that are immune to computational problems.

Then, we need to search for threshold $¯ k$ such that the empirical distribution ${ˆ f}_{Λ}$ of the smaller eigenvalues ${ˆ λ}_{k} < {ˆ λ}_{¯ k}$ best fits the rescaled Marchenko-Pastur law in the empirical pdf (24.94), at which point we cannot reject that the smaller sample eigenvalues stem from the noise component of process $ε_{t}$ , and therefore they must be all equal

V : λ_{¯ k + 1} = \dots = λ_{¯ ı} .

(24.95)

The statistical interpretation of the view (24.95) is the following. The larger eigenvalues cannot have been generated by noise. However, for the smaller eigenvalues we cannot rule out that they stem from noise, i.e. a data generating process whose true covariance is a multiple $λ_{n o i s e}$ of the identity matrix and thus has all equal eigenvalues.

The geometric interpretation of the view (24.95) is the following: the principal directions (8.71) $e_{k}$ identified by the smallest eigenvalues $λ_{k}$ ( $k > ¯ k$ ) of the covariance matrix $σ^{2}$ are spurious and thus it makes more sense to assume isotropy, or spherical symmetry. Indeed, the decomposition is not affected by arbitrary rotations $o$ applied to the noise component in the i.i.d. process (24.92).

To summarize, the spectrum shrinkage stemming from random matrix theory can be performed in four steps.

First, we perform the PCA decomposition (43.413) on the estimated sample covariance ${ˆ σ}^{2}$ (24.117), in order to compute the eigenvectors matrix $ˆ e$ (43.369) and the sample spectrum $ˆ λ$ (43.356).

Second, we determine the index $¯ k$ such that the rescaled Marchenko-Pastur law $f_{Λ}^{M P}$ (24.90) best fits the empirical distribution of the smaller eigenvalues ${ˆ λ}_{k} < {ˆ λ}_{¯ k}$ : the fit can be measured by any divergence $D$ between the two distributions.

Third, we shrink the spectrum by imposing isotropy (24.95).

Fourth, we combine the eigenvectors $ˆ e$ with the spectrum $¯ ¯ ¯ λ$ to obtain the new covariance matrix ${¯ ¯ ¯ σ}^{2}$ .

The corresponding routine reads:

The spectrum shrinkage routine (Table 24.9) can be applied at first impression to any covariance matrix ${ˆ σ}^{2}$ . However, if the input is not the sample covariance ${ˆ σ}^{2}$ (24.117) of a time series of i.i.d variables ${ϵ_{1}, \dots, ϵ_{¯ t}}$ , the optimization in Step 3 must be modified accordingly, see e.g. [Potters et al., 2005].

Example 24.25. Spectrum shrinkage via random matrix theory

In Figure 24.17 we highlight the main features of the covariance shrinkage based on the Marchenko-Pastur distribution from random matrix theory.

In the top-left plot we show the sample spectrum (43.389) computed from the sample covariance (24.78). Then we assume that the smallest eigenvalues are noisy and thus their relative size bears no meaning and we shrunk their value to their mean (24.95). Here we show the empirical distribution of the small sample eigenvalues.

In the top-center plot we also show the rescaled Marchenko-Pastur distribution (24.90) that is supposed to fit the noisy sample eigenvalue distribution (24.94). Furthermore, we highlight the spike toward which the noisy eigenvalues are shrunk which is supposed to be their true population spectrum.

In the top-right plot we show the mean-covariance ellipsoids (8.47) defined by the sample covariance and by the shrinkage covariance respectively, sliced along the plane spanned by a high-volatility direction, namely the 10-th sample eigenvector; and a low-volatility direction, namely the last eigenvector. The direction of the last eigenvector is statistically the same as the direction of all the eigenvectors whose value has been shrunk, due to the isotropy implied by the shrinkage to the mean (24.95).

In the bottom plot we display all the original and filtered correlations sorted in decreasing order. The overall effect of the eigenvalue shrinkage is a mild shrinkage of the correlations toward zero.

As we vary the cutoff that divides noise from the signal, the fit of the Marchenko-Pastur distribution to the distribution of the noisy empirical eigenvalues varies. We need to select a value that can be accepted as consistent Marchenko-Pastur distribution.

When all the eigenvalues are shrunk to isotropy, we obtain the Ledoit-Wolf shrinkage target (24.84).

24.6.4 Covariance shrinkage: sparse eigenvector rotations

In Sections 24.6.2-24.6.3 we used shrinkage to provide an accurate spectral decomposition (43.413) of the true, unknown covariance matrix of a set of i.i.d. variables $ε_{t}$ (36.2)

C v {ε_{t}} \equiv σ^{2} = e D i a g (λ) e^{'} .

(24.96)

More precisely, we flattened the spuriously steep spectrum of the sample estimate of the covariance matrix by averaging the sample estimate with a suitable target matrix (24.84), or by averaging the smallest eigenvalues (24.95).

The approach in [Cao et al., 2011] to improve the estimate of the spectral decomposition (24.96) acts on the eigenvectors matrix $e$ , rather than on the spectrum of eigenvalues $λ$ (43.356).

We start with the sample covariance ${ˆ σ}^{2} \equiv ˆ e D i a g (ˆ λ) {ˆ e}^{'}$ (24.117) of a set of i.i.d. variables $ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'}$ , computed from a time series $i = {ϵ_{1}, \dots, ϵ_{¯ t}}$ of independent realizations of $ε_{t}$ . We stress the fact that it is always advisable to standardize (54b.83) the observations $ϵ_{t}$ before focusing on the estimation of the covariance, in order to obtain correlation-like estimates that are immune to computational problems.

To define the shrinkage views $V$ , let us recall that the eigenvectors matrix $e$ is orthogonal (43.398) and hence represents a rotation (or a reflection) [W]. Then we set shrinkage views imposing that the eigenvectors matrix $e$ is a sequence of $¯ k$ simple rotations along the coordinate planes

V : e = r (i_{1}, i_{1}^{'}, α_{1}) \times \dots \times r (i_{¯ k}, i_{¯ k}^{'}, α_{¯ k}),

(24.97)

where each rotation is a Givens rotation [W], that only involves two indices and reads entrywise

r (i, j, α) \equiv ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} cos α & entries i, i and j, j sin α & entry i, j - sin α & entry j, i 1 & all remaining diagonal elements 0 & all remaining elements. \end{matrix} .

(24.98)

We then define the shrinkage estimate of the covariance as the matrix that maximizes the conditional likelihood $f (i | σ^{2})$ (24.38) of the time series $i$ under the normal distribution assumption $ε_{t} \sim N (μ, σ^{2})$ , all while satisfying the constraint (24.97) on the eigenvectors

\begin{matrix} {¯ ¯ ¯ σ}^{2} = & argmax σ^{2} f (i | σ^{2}) & (24.99) subject to {\begin{matrix} σ^{2} \equiv e D i a g (λ) e^{'} e \in V . \end{matrix} \end{matrix}

Maximizing the target in (24.99) is equivalent to finding the eigenvectors $e$ that minimize the following “distance” with respect to the sample eigenvectors $ˆ e$ 70.26

D (e, ˆ e) \equiv det (D i a g (d i a g (e^{'} {ˆ σ}^{2} e))) .

(24.100)

Then the eigenvectors that determine shrinkage estimate $σ^{2}$ , which satisfies the shrinkage views $V$ (24.97) and at the same time is as close as possible to the original estimate ${ˆ σ}^{2}$ , become

¯ e \equiv argmin e \in V D (e, ˆ e) .

(24.101)

The shrinkage (24.101) under the rotation views (24.97) is not computationally tractable. However, [Cao et al., 2011] provide a heuristic approximation to the true optimum.

24.6.5 Covariance shrinkage: glasso

In this section we present an approach that shrinks the covariance $C v {ε}$ (8.16) and correlation $C r {ε}$ (10.15) matrices of a set of i.i.d. variables (36.2) to a parsimonious undirected graphical model or Markov random field (21.82).

As explained in Section 21.7, a set of simultaneous i.i.d. variables $ε_{t} = (ε_{1}, \dots, ε_{¯ ı})$ is a parsimonious Markov random field when the number of $0$ -entries in the inverse correlation $ϕ^{2} \equiv {(c^{2})}^{- 1}$ is large.

Example 24.26. Refer to Example 21.5 to visualize a simple graphical model.

We start from the sample inverse correlation with flexible probabilities ${ˆ ϕ}^{2 H F P}$ , i.e. the inverse of the correlation matrix ${ˆ c}^{2 H F P}$ (24.8), estimated on a time series $i = {ϵ_{1}, \dots, ϵ_{¯ t}}$ of realizations of the i.i.d. variables $ε_{t}$ .

In order to shrink the market to a parsimonious network, we impose the view that the inverse correlations generate a parsimonious network with few connections. Hence, we want the null entries of the inverse correlations to be at least a given lower boundary $k$

V \equiv {ϕ^{2} : (num of 0 -entries in ϕ^{2}) \geq k} .

(24.102)

We then find the parameter ${¯ ¯ ¯ ϕ}^{2}$ that minimizes the distance from ${ˆ ϕ}^{2 H F P}$ , all while satisfying the view (24.102)

{¯ ¯ ¯ ϕ}^{2} \equiv argmin ϕ^{2} \in V D (ϕ^{2}, {ˆ ϕ}^{2 H F P}) .

(24.103)

where $D$ is the pseudo-distance (E.22.208). The target in (24.103) is equivalent to the relative entropy under normal assumptions (5.75).

The shrinkage views (24.102) are a cardinality constraint (45.138). Minimizing the distance (E.22.208) under such a constraint is impossible in practice. However, we can resort to quasi-optimal heuristics. One such heuristic is via greedy algorithms [W], in particular forward selection, see Section 47.2.3.

A moderately less optimal, although more efficient, solution is the graphical lasso or glasso, see [Hastie et al., 2009]. Similarly to the regular lasso discussed in Section 25.5, the graphical lasso minimizes the distance (E.22.208) from a generic symmetric and positive definite matrix ${c - -}^{2}$ , while penalizing the sum of the absolute values of the entries $∥ ϕ^{2} ∥_{1} \equiv \sum_{i, j} | {[ϕ^{2}]}_{i, j}^{2} |$ , as follows

ϕ_{λ}^{2} \equiv argmin ϕ^{2} tr (ϕ^{2} {c - -}^{2}) - ln det (ϕ^{2}) + λ ∥ ϕ^{2} ∥_{1} .

(24.104)

As in the standard lasso, the penalty $λ$ forces some entries of $ϕ^{2}$ to shrink to zero, and the number of zero entries increases as the penalty $λ$ increases, see Example 24.27.

Example 24.27. Stocks Markov network structure via graphical lasso of correlation

We consider the same setting as in Example 24.19. We apply the glasso shrinkage, for different penalty values $λ$ , on the inverse ${ˆ ϕ}^{2}$ of the sample correlation ${ˆ c}^{2}$ stemming from ${ˆ σ}^{2}$ (10.22), as described in Table 24.10.
In the top plot of Figure 24.18 we show the evolution of the graph of the shrunk inverse correlation matrix $ϕ_{λ}^{2}$ (24.104).
In the bottom-left plots we display the heat map of the HFP correlation matrix ${ˆ s}^{2} \leftarrow c o r r ({ˆ s}^{2})$ (24.3)-(10.22) and its inverse ${ˆ ϕ}^{2}$ .
In the bottom-right plots we display the heat map of the correlation matrix stemming from the regularized inverse covariance $ϕ_{λ}^{2}$ (24.104).
The elements of the inverse-correlation terms become null one at a time as the penalty increases. When all the elements of a column/row (except for the diagonal element) of the inverse matrix have become null, then the same column/row disappears also in the covariance matrix. This is a consequence of the formula on partitioned matrix inversion (43.615).
To implement the glasso we used the code adapted by H. Karshenas from the original code provided by [Friedman et al., 2008].
Refer to Example 24.19 to see how the same parsimonious structure of the Gaussian Markov random field (21.85) can be obtained through Bayesian shrinkage.

To summarize, the routine to perform shrinkage to a parsimonious network proceeds as follows

24.6.6 Covariance shrinkage: factor analysis

In this section we present factor analysis, used to shrink the covariance matrix $C v {ε}$ (8.16) or correlation matrix $C r {ε}$ (10.15) of a set of i.i.d. variables (36.2) to a low-rank-diagonal structure (43.608).

More precisely, let us define the low-rank diagonal function

l r d (b, d) \equiv b b^{'}    low-rank + D i a g (d)      diagonal,

(24.105)

where $b$ is a full-rank $¯ n \times ¯ k$ matrix (typically $¯ k ≪ ¯ n$ ) and $d$ is an $¯ n \times 1$ positive vector.

Then, given a reference $¯ n \times ¯ n$ symmetric, positive-definite matrix $σ^{2}$ , e.g. the HFP covariance (24.3), the output ${¯ ¯ ¯ σ}^{2}$ is a low-rank diagonal matrix $l r d (β, δ)$ (24.105) determined by the matrix $β$ and vector $δ$ that solve the following optimization problem

(β, δ) \equiv argmin b, d D (l r d (b, d), σ^{2}) .

(24.106)

for a suitable divergence $D$ .

If the divergence $D$ is specified as the distance induced (43.263) by the Frobenius norm (43.256), the optimization is solved via the principal axis factorization, as discussed in Section 14.4.2.

If $D$ is specified as the relative entropy between normal variables (5.75), the optimization is solved by the maximum likelihood factorization, as discussed in Section 21.3.2.

Refer also to [MacCallum et al., 2007] and [Cudeck and MacCallum, 2012] for a discussion on the performance of the two algorithms.

Example 24.28. Let us consider the following $3 \times 3$ positive-definite covariance matrix

σ^{2} \equiv ⎛ ⎜ ⎝ \begin{matrix} 3.10 & 0.36 & 0.16 0.36 & 4.31 & 0.28 0.16 & 0.28 & 0.88 \end{matrix} ⎞ ⎟ ⎠ .

(24.107)

We shrink $σ^{2}$ towards a low-rank diagonal matrix with $¯ k = 1$ (24.105) using the principal axis factorization algorithm (PAF) (Table 14.3)

σ^{2; P A F} \equiv l r d (β^{P A F}, δ^{P A F}) \approx σ^{2};

(24.108)

and the maximum likelihood factorization algorithm (MLF) (Table 21.1)

σ^{2; M L F} \equiv l r d (β^{M L F}, δ^{M L F}) \approx σ^{2} .

(24.109)

Then we compare the distance of each factor analysis matrix $σ^{2; P A F}$ and $σ^{2; M L F}$ with respect to the original covariance matrix $σ^{2}$ using both i) the Frobenius distance (43.263)-(43.256) and ii) the (normal) relative entropy (5.75)

\begin{matrix} Frobenius distance & Relative Entropy - ---- - & - ------------------- - & - ----------------- - PAF & 1 0^{- 6} \times 0.18 & 1 0^{- 14} \times 0.13 MLF & 1 0^{- 6} \times 0.19 & 1 0^{- 15} \times 0.9 \end{matrix}

(24.110)

¯ ¯ ¯ μ \equiv (1 - γ) ˆ μ + γ μ_{t a r g e t},

{¯ ¯ ¯ σ}^{2} = (1 - γ) {ˆ σ}^{2} + γ σ_{t a r g e t}^{2} .

V \equiv {ϕ^{2} : (num of 0 -entries in ϕ^{2}) \geq k} .

l r d (b, d) \equiv b b^{'}    low-rank + D i a g (d)      diagonal,

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} X_{1} \cdot X_{n} \cdot X_{¯ n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠      X = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} α_{1} \cdot α_{n} \cdot α_{¯ n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠      α + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} β_{1, 1} & \dots & β_{1, k} & \dots & β_{1, ¯ k} \cdot & \cdot & \cdot β_{n, 1} & \dots & β_{n, k} & \dots & β_{n, ¯ k} \cdot & \cdot & \cdot β_{¯ n, 1} & \dots & β_{¯ n, k} & \dots & β_{¯ n, ¯ k} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠      β ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} Z_{1} \cdot Z_{k} \cdot Z_{¯ k} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠      Z + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} {˚ ε}_{1} \cdot {˚ ε}_{n} \cdot {˚ ε}_{¯ n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠      ˚ ε,

E {X} \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} E {X_{1}} \cdot E {X_{n}} \cdot E {X_{¯ n}} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠,

ˆ μ \equiv e s_{ˆ μ} (i) \equiv \frac{1}{¯ t} \sum_{t = 1}^{¯ t} ϵ_{t};

ε_{t} \sim F (ε) i.i.d..

λ_{1} \geq \dots \geq λ_{¯ n} \geq 0 .

γ = min {1, \frac{2}{¯ t} \frac{tr (σ^{2}) - 2 λ_{1}}{{(ˆ μ - μ_{t a r g e t})}_{t a r g e t}^{'} (ˆ μ - μ_{t a r g e t})}},

\begin{matrix} {ˆ s}_{ε}^{2 H F P} & \equiv {ˆ C v}^{H F P} {ε} = \sum_{t = 1}^{¯ t} p_{t} (ϵ_{t} - {ˆ m}_{ε}^{H F P}) {(ϵ_{t} - {ˆ m}_{ε}^{H F P})}_{t}^{'} & (24.3) = \sum_{t = 1}^{¯ t} p_{t} ϵ_{t} ϵ_{t}^{'} - {ˆ m}_{ε}^{H F P} {({ˆ m}_{ε}^{H F P})}^{'} . \end{matrix}

ε_{t} \sim N (μ, σ^{2}) .

\begin{matrix} f (i | θ) & = exp (¯ t (θ_{μ}^{'} {ˆ η}_{μ} + tr (θ_{σ} {ˆ η}_{σ}) - ψ (θ) + ln h)) & (24.40) = {(2 π)}^{- \frac{¯ ı ¯ t}{2}} det {(σ^{2})}^{- \frac{¯ t}{2}} exp (- \frac{¯ t}{2} [tr ({ˆ s}^{2} {(σ^{2})}^{- 1}) + {(ˆ m - μ)}^{'} {(σ^{2})}^{- 1} (ˆ m - μ)]), \end{matrix}

{(Σ^{2})}^{- 1} \sim W i s h a r t (ν^{p r i}, \frac{1}{ν^{p r i}} {(σ^{2 p r i})}^{- 1});

M | σ^{2} \sim N (μ^{p r i}, \frac{1}{t^{p r i}} σ^{2}),

f^{p o s} (θ) \equiv f (θ | i) = \frac{f (i | θ) f^{p r i} (θ)}{\int f (i | ϑ) f^{p r i} (ϑ) d ϑ},

{(Σ^{2})}^{- 1} | i \sim W i s h a r t (ν^{p o s}, \frac{1}{ν^{p o s}} {(σ^{2 p o s})}^{- 1}),

M | σ^{2}, i \sim N (μ^{p o s}, \frac{1}{t^{p o s}} σ^{2}) .

μ^{p o s} = γ_{μ} μ^{p r i} + (1 - γ_{μ}) ˆ m,

ˆ m large dataset ¯ t \leftarrow μ^{c l_e q} = μ^{p o s} high confidence t^{p r i} \to μ^{p r i} .

ℓ o s s (θ, ˆ q) \equiv ∥ q (θ) - ˆ q ∥_{2}^{2} .

r (θ, e s) \equiv E {D (q (Θ) ∥ e s (I)) | θ} .

X \sim N (μ, σ^{2}),

e i g (q) \equiv {e i g v a l (q), e i g v e c (q)} = {λ, e},

i \equiv (ϵ_{1}, \dots, ϵ_{t}, \dots, ϵ_{¯ t}) \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} ϵ_{1, 1} & ϵ_{1, t} & ϵ_{1, ¯ t} \cdot & \cdot & \cdot ϵ_{i, 1} & \dots & ϵ_{i, t} & \dots & ϵ_{i, ¯ t} \cdot & \cdot & \cdot ϵ_{¯ ı, 1} & ϵ_{¯ ı, t} & ϵ_{¯ ı, ¯ t} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠,

ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'} \sim f_{ε} i.i.d., with ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} E {ε_{i, t}} = 0 V {ε_{i, t}} = 1 C v {ε_{i, t}, ε_{j, t}} = 0, for i \neq j \end{matrix} .

{ˆ σ}^{2} \equiv e s_{{ˆ σ}^{2}}^{2} (i) \equiv \frac{1}{¯ t} \sum_{t = 1}^{¯ t} (ϵ_{t} - ˆ μ) {(ϵ_{t} - ˆ μ)}_{t}^{'} .

{ˆ λ, ˆ e} = e i g ({ˆ σ}_{ε}^{2}) .

{λ, e} = e i g (C v {ε_{t}}),

\sqrt{e i g v a l (C v {X})} = principal semi-axes length \Leftrightarrow shape of ellipsoid.

\partial E (M o d {X}, M o D i s^{2} {X}) \equiv {x : {(x - M o d {X})}^{'} {(M o D i s^{2} {X})}^{- 1} (x - M o d {X}) = 1} .

ε_{t} \equiv {(ε_{1, t}, \dots, ε_{¯ ı, t})}_{1, t}^{'} \sim N (0_{¯ ı \times 1}, I_{¯ ı}) i.i.d..

C v {ε_{t}} = I_{¯ ı} \Rightarrow λ = 1_{¯ ı} .

\partial E (E {X}, C v {X}) = {x \in R^{¯ n} : {(x - E {X})}^{'} {(C v {X})}^{- 1} (x - E {X}) = 1} .

Z_{n}^{P C} = argmax Z \in C_{{0, . . ., n - 1}} V {Z},

{λ_{i}}_{i = 1}^{¯ ı} \Leftrightarrow sphere estimation \to {{ˆ λ}_{i}}_{i = 1}^{¯ ı} \Leftrightarrow cigar.

¯ t < ¯ ı \Rightarrow r a n k ({ˆ σ}^{2}) = ¯ t < ¯ ı,

∥ a ∥_{F} \equiv ∥ v e c (a) ∥_{2} = \sqrt{tr (a^{'} a)} = \sqrt{\sum_{n = 1}^{min (¯ n, ¯ ¯ ¯ m)} γ_{n}^{2}},

s^{2} = e \times D i a g (λ) \times e^{'},

{¯ f}_{X}^{W O E} (λ) = f_{Λ}^{M P} (λ) \equiv (1 - \frac{1}{q}) δ (λ) 1_{q > 1} + \frac{\sqrt{(λ - λ^{-}) (λ^{+} - λ)}}{2 π λ q σ^{2}} 1_{λ \in (λ^{-}, λ^{+})}^{-},

{ˆ f}_{Λ} (y | λ) \equiv \frac{1}{¯ n} \sum_{n = 1}^{¯ n} δ (y - λ_{n}),

{ˆ f}_{Λ} (x | ˆ λ) \equiv \frac{1}{¯ ı} \sum_{i = 1}^{¯ ı} δ_{i}^{({ˆ λ}_{i})} (x) .

lim ¯ n \to \infty ˆ f (λ | λ (x^{M P})) = f_{Λ}^{M P} (λ) .

q \equiv \frac{¯ ı}{¯ t}, λ^{\pm} \equiv {(1 \pm \sqrt{q})}^{2}, σ^{2} \equiv 1 .

{ˆ f}_{Λ} (x | ˆ λ) ¯ ı \equiv q ¯ t \to \infty \approx f_{Λ}^{M P} (x),

f_{Λ}^{M P} (x) \equiv (1 - \frac{1}{q}) δ^{(0)} (x) 1_{q > 1} + \frac{\sqrt{(x - λ^{-}) (λ^{+} - x)}}{2 π x q} 1_{x \in (λ^{-}, λ^{+})}^{-},

{¯ f}_{X}^{G O E} (λ) = f_{Λ}^{S C} (λ) \equiv \frac{1}{π σ^{2}} \sqrt{2 σ^{2} - λ^{2}} 1_{λ \in (- \sqrt{2} σ, \sqrt{2} σ)} .

ˆ h i s t_{X} (x) \equiv \sum_{j} {ˆ f}_{Δ}^{(j)} 1_{x \in X^{(j)} .}^{(j)}

o \in O (¯ n) \Leftrightarrow o^{'} o = o o^{'} = I_{¯ n} .

⎛ ⎜ ⎝ \begin{matrix} ε_{1} \cdot ε_{¯ ı} \end{matrix} ⎞ ⎟ ⎠ = (e_{1}, \dots, e_{¯ k}) \times ⎛ ⎜ ⎜ ⎝ \begin{matrix} ε_{1}^{s i g n l} \cdot ε_{¯ k}^{s i g n l} \end{matrix} ⎞ ⎟ ⎟ ⎠ + (e_{¯ k + 1}, \dots, e_{¯ ı}) \times o \times ⎛ ⎜ ⎜ ⎝ \begin{matrix} ε_{¯ k + 1}^{n o i s e} \cdot ε_{¯ ı}^{n o i s e} \end{matrix} ⎞ ⎟ ⎟ ⎠,

e \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} e_{1, 1} ⋮ e_{n, 1} ⋮ e_{¯ n, 1} \end{matrix} ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ \dots ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ \begin{matrix} e_{1, n} ⋮ e_{n, n} ⋮ e_{¯ n, n} \end{matrix} ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ \dots ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ \begin{matrix} e_{1, ¯ n} ⋮ e_{n, ¯ n} ⋮ e_{¯ n, ¯ n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ \equiv (e_{1} | \dots | e_{n} | \dots | e_{¯ n}),

{ˆ f}_{Λ} (x | ˆ λ) ¯ ı \equiv q ¯ t \to \infty \approx \frac{1}{λ} f_{Λ}^{M P} (\frac{x}{λ}),

\int_{T} δ^{(u)} (t) g (t) d μ (t) = g (u) .

e i g v a l (C v {ε_{t}}) = {(λ_{1}, \dots, λ_{¯ k}, λ, \dots, λ)}_{1}^{'},

{ˆ f}_{Λ} (x | ˆ λ) \approx w \times \frac{1}{¯ k} \sum_{k \leq ¯ k} δ_{k}^{({ˆ λ}_{k})} (x) + (1 - w) \times \frac{1}{λ_{n o i s e}} f_{Λ}^{M P} (\frac{x}{λ_{n o i s e}}),

{˜ ϵ}_{i, t} \equiv Φ_{ν}^{- 1} (u_{i, t}),

V : λ_{¯ k + 1} = \dots = λ_{¯ ı} .

e_{n} \equiv argmax v \in V_{{0, . . ., n - 1}} V {v^{'} X},

λ \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} λ_{1} ⋮ λ_{¯ k} λ_{¯ k + 1} ⋮ λ_{¯ n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} λ_{1} ⋮ λ_{¯ k} λ_{1}^{Re} + i λ_{1}^{Im} λ_{1}^{Re} - i λ_{1}^{Im} ⋮ λ_{¯ ȷ}^{Re} + i λ_{¯ ȷ}^{Im} λ_{¯ ȷ}^{Re} - i λ_{¯ ȷ}^{Im} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

C v {X} = e \times D i a g (λ) \times e^{'},

λ_{n} \in R,

C v {ε_{t}} \equiv σ^{2} = e D i a g (λ) e^{'} .

{˜ ε}_{i, t} \equiv \frac{ε_{i, t} - {ˆ μ}_{i}}{{ˆ σ}_{i}},

e e^{'} = e^{'} e = I_{¯ n} .

f (i | θ) i.i.d. = \prod_{t = 1}^{¯ t} f (ϵ_{t} | θ),

V : e = r (i_{1}, i_{1}^{'}, α_{1}) \times \dots \times r (i_{¯ k}, i_{¯ k}^{'}, α_{¯ k}),

\begin{matrix} {¯ ¯ ¯ σ}^{2} = & argmax σ^{2} f (i | σ^{2}) & (24.99) subject to {\begin{matrix} σ^{2} \equiv e D i a g (λ) e^{'} e \in V . \end{matrix} \end{matrix}

¯ e \equiv argmin e \in V D (e, ˆ e) .

C v {X} \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} V {X_{1}} & C v {X_{1}, X_{2}} & \cdot & C v {X_{1}, X_{¯ n}} C v {X_{2}, X_{1}} & V {X_{2}} & \cdot & \cdot \cdot & \cdot & \cdot & \cdot C v {X_{¯ n}, X_{1}} & \cdot & \cdot & V {X_{¯ n}} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ .

c_{X}^{2} \equiv C r {X} \equiv ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 1 & C r {X_{1}, X_{2}} & \cdot & C r {X_{1}, X_{¯ n}} C r {X_{2}, X_{1}} & 1 & \cdot & C r {X_{2}, X_{¯ n}} \cdot & \cdot & \cdot & \cdot C r {X_{¯ n}, X_{1}} & \cdot & \cdot & 1 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ .

(Y_{m} ⊥ ⊥ Y_{n}) | Y_{- {m, n}} \Leftrightarrow (m \sim n) \notin E,

{ˆ c}_{ε}^{2 H F P} = c o r r ({ˆ s}_{ε}^{2 H F P}) (10.22) = D i a g (1_{¯ ı \times 1} . ∕ {ˆ s}_{ε; s t d e v}^{H F P}) {ˆ s}_{ε}^{2 H F P} D i a g (1_{¯ ı \times 1} . ∕ {ˆ s}_{ε; s t d e v}^{H F P}) .

D (ϕ^{2}, {ϕ - - -}_{X}^{2}) \equiv tr (ϕ^{2} ({ϕ - - -}_{X}^{2})^{- 1}) - ln (det (ϕ^{2})) .

{¯ ¯ ¯ ϕ}^{2} \equiv argmin ϕ^{2} \in V D (ϕ^{2}, {ˆ ϕ}^{2 H F P}) .

D_{K L} (μ_{1}, σ_{1}^{2} ∥ μ_{0}, σ_{0}^{2}) = \frac{1}{2} (tr (σ_{1}^{2} {(σ_{0}^{2})}^{- 1}) - ln (det (σ_{1}^{2} {(σ_{0}^{2})}^{- 1})) + {(μ_{1} - μ_{0})}_{1}^{'} {(σ_{0}^{2})}^{- 1} (μ_{1} - μ_{0}) - ¯ n) .

k - - \leq c a r d ({n such that x_{n} \neq 0}) \leq ¯ k,

c o r r (σ^{2}) \equiv D i a g (\frac{1_{¯ n \times 1}}{v o l (σ^{2})}) \times σ^{2} \times D i a g (\frac{1_{¯ n \times 1}}{v o l (σ^{2})}),

ϕ_{λ}^{2} \equiv argmin ϕ^{2} tr (ϕ^{2} {c - -}^{2}) - ln det (ϕ^{2}) + λ ∥ ϕ^{2} ∥_{1} .

{[a^{- 1}]}_{\neq n, n}^{- 1} = - {[a^{- 1}]}_{n, n}^{- 1} {[a]}_{\neq n, \neq n}^{- 1} {[a]}_{\neq n, n} .

Y \sim N (μ, σ^{2}),

s_{¯ n}^{2} = β s_{¯ k}^{2} β^{'} + D i a g (δ ⊙ δ),

D (x, y) \equiv ∥ x - y ∥,


Inputs	${ˆ σ}^{2}$

1. Perform spectral decomposition	$(ˆ e, ˆ λ) spectral dec. \Leftarrow {ˆ σ}^{2} = ˆ e D i a g (ˆ λ) {ˆ e}^{'}$	(8.42)
for $k = 1, . . ., ¯ ı$
2a. Compute empirical eigenvalue distribution	${ˆ f}_{Λ}^{(k)} \Leftarrow {ˆ λ}_{k + 1}, \dots, {ˆ λ}_{¯ ı}$
2b. Compute rescaled Marchenko-Pastur distribution	${˜ f}_{Λ}^{(k)} \leftarrow \frac{1}{λ_{n o i s e}} f_{Λ}^{M P} (\frac{x}{λ_{n o i s e}})$
end
3. Assess divergence for optimal cutoff	$¯ k \leftarrow argmin k \in {1, . . ., ¯ ı - 1} D ({˜ f}_{Λ}^{(k)} ∥ {ˆ f}_{Λ}^{(k)})$
4. Impose isotropy	$⎧ ⎨ ⎩ \begin{matrix} {¯ λ}_{i} \leftarrow {ˆ λ}_{i} & for all i \leq ¯ k {¯ λ}_{i} \leftarrow \frac{1}{¯ ı - ¯ k} \sum_{k > ¯ k} {ˆ λ}_{k} & for all i > ¯ k \end{matrix}$
5. Compute new covariance	${¯ ¯ ¯ σ}^{2} \leftarrow ˆ e D i a g (¯ ¯ ¯ λ) {ˆ e}^{'}$

Outputs	${¯ ¯ ¯ σ}^{2}$


Inputs	${σ - - -}^{2}, k, {λ_{ι}}$

1. Compute correlation	${c - -}^{2} \leftarrow D i a g (1_{¯ ı \times 1} . ∕ {σ - - -}_{v o l}) {σ - - -}^{2} D i a g (1_{¯ ı \times 1} . ∕ {σ - - -}_{v o l})$
for $ι = 1, . . ., ¯ ı$
2. Perform graphical lasso	${˜ ϕ}_{λ_{ι}}^{2} graphical lasso \Leftarrow ({c - -}^{2}, λ_{ι})$	(24.104)
3. Extract correlation	$⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} {˜ c}_{λ_{ι}}^{2} \leftarrow {({˜ ϕ}_{λ_{ι}}^{2})}^{- 1} c_{λ_{ι}}^{2} \leftarrow D i a g (1_{¯ ı \times 1} . ∕ {[{˜ c}_{λ_{ι}}]}_{v o l}) {˜ c}_{λ_{ι}}^{2} D i a g (1_{¯ ı \times 1} . ∕ {[{˜ c}_{λ_{ι}}]}_{v o l}) ϕ_{λ_{ι}}^{2} \leftarrow {(c_{λ_{ι}}^{2})}^{- 1} = D i a g ({[{˜ c}_{λ_{ι}}]}_{v o l}) {˜ ϕ}_{λ_{ι}}^{2} D i a g ({[{˜ c}_{λ_{ι}}]}_{v o l}) \end{matrix}$
end
4. Compute optimal penalty	$¯ λ \leftarrow min {λ_{ι}$ such that $(num of 0 -entries in ϕ_{λ_{ι}}^{2}) \geq k}$
5. Shrunk inv. correlation	${¯ ¯ ¯ ϕ}^{2} \leftarrow ϕ_{¯ λ}^{2}$
6. Shrunk correlation	${¯ c}^{2} \leftarrow {({¯ ¯ ¯ ϕ}^{2})}^{- 1}$
7. Shrunk covariance	${¯ ¯ ¯ σ}^{2} \leftarrow D i a g ({σ - - -}_{v o l}) {¯ c}^{2} D i a g ({σ - - -}_{v o l})$

Outputs	${¯ ¯ ¯ σ}^{2}, {¯ c}^{2}, {¯ ¯ ¯ ϕ}^{2}$