Shrinkage Estimation of Mean and Covariance

The Problem

Every classical allocation recipe — Markowitz, minimum-variance, risk parity — needs two inputs that are never observed: the expected-return vector $\boldsymbol\mu$ and the covariance matrix $\Sigma$. In practice their sample estimates are plugged in as if they were the truth, and the counting alone shows why that hurts: $N$ assets require $N$ means plus $N(N+1)/2$ covariance entries — 1,325 numbers for 50 assets, typically from a few hundred observations. The estimates are unbiased but wildly noisy, and a mean–variance optimiser is an error-maximiser: it loads up precisely on whatever was over-estimated by chance. Shrinkage accepts a little bias to kill a lot of variance, pulling the noisy sample toward a simple, stable target.

Shrinking the mean — Stein's paradox

For the mean, the justification is Stein's paradox: for $N\ge3$ the sample mean is inadmissible — a shrunk estimator beats it in total mean-squared error, always. The notebook demonstrates that by Monte Carlo rather than citation, and the advantage is real but conditional: at $N=10$, $T=30$ the James–Stein estimator cuts total squared error by 23%, and the gain grows as $T$ shrinks and vanishes as data accumulates, with the optimal intensity heading to zero. On the real sector data the intensity hits its cap at $a=1$ — full shrinkage to the grand mean — which is itself the finding: at this sample size the cross-sectional spread in mean returns is indistinguishable from noise.

Shrinking the covariance — conditioning

For the covariance the problem is conditioning. The sample covariance over-disperses its eigenvalues — the largest too large, the smallest crushed toward zero — and it is the inverse that every optimiser needs. Ledoit–Wolf shrinkage repairs the spectrum with a closed-form, tuning-free intensity, and the notebook verifies its own implementation against the derivation, computing the intensity explicitly and matching the engine to within $10^{-8}$ ($a=0.1735$). The effect on conditioning is dramatic: on the sector panel the condition number falls from 96 to 50; on 48 stocks at $N/T=0.8$ it falls from 13,174 to 117.

When Shrinkage Helps — and When It Does Not

What raises this above a standard shrinkage demo is that it reports when shrinkage loses. A minimum-variance backtest sweeping the estimation window from data-scarce to data-rich shows the sample portfolio realising 3.51% weekly volatility when $T$ barely exceeds $N$ against 2.13% for the best shrinkage estimator — a 65% penalty for inverting a near-singular covariance. But at a 52-week window the sample is fine (2.30%) and shrinking toward the identity buys no volatility advantage at all (2.31%). The honest conclusion is drawn: shrinkage is indispensable when data is scarce and stops paying for itself when it is not. The target matters too, and instructively the two experiments disagree about which one to prefer — on the synthetic equicorrelated market a constant-correlation target cuts covariance-estimation error by 20% where the identity gains nothing, yet on the real sector panel the identity target delivers the lower realised volatility at every window. Matching the target to the market's structure helps when you are estimating the matrix; it is not the same objective as minimising a portfolio's out-of-sample variance.

Dimensionality is the amplifier

Dimensionality is the amplifier, and the project makes that the closing argument by moving from ten sectors ($N/T\approx0.1$) to 48 stocks and then to the canonical Fama–French 48 industries. There the verdict is unambiguous on both universes: at a 60-month window the sample minimum-variance portfolio realises 22.7% annualised volatility with turnover of 3.84, against 12–13% and turnover under 0.5 for the shrinkage estimators. Two wins hold unconditionally even where the sample is competitive on volatility — shrinkage roughly halves turnover at every window, a direct trading-cost saving, and it never blows up.

Shrinkage is Bayes

A third notebook makes the connection the formula already hints at: $\hat\theta_{\text{shrunk}}=(1-a)\hat\theta+a\,\theta_{\text{target}}$ is the shape of a Bayesian posterior mean, so shrinkage is Bayes. A hierarchical PyMC model recovers the same partial pooling toward the grand mean that James–Stein produces as a point estimate, but with a full posterior attached. A conjugate Inverse-Wishart prior does the same for the covariance, sweeping the prior strength from sample to target, and remains analytic and instant at scale where NUTS on a full LKJ model would not be practical. One instructive non-result surfaces. Run head-to-head on the same sector panel, the LKJ posterior and Ledoit–Wolf repair conditioning by almost exactly the same amount — sample $51.0$ against Ledoit–Wolf $28.3$ and LKJ $28.0$, so $1.80\times$ and $1.82\times$ respectively. The two routes disagree about how — LKJ shrinks correlations toward zero while an identity target acts directly on the eigenvalue spectrum conditioning actually measures — but on ten well-behaved sectors that distinction does not show up in the condition number. Where Ledoit–Wolf pulls decisively ahead is dimensionality, not method: on the 48-stock panel at $N/T=0.8$ it improves conditioning by a factor of 113, a regime the LKJ sampler was not run in.

Notebooks

Downloads

Shrinkage Module — Source Code

"""
shrinkage.py  --  Shrinkage estimators of mean and covariance for asset returns.

This engine backs the notebooks in
    Shrinkage Estimation of Mean and Covariance

It is deliberately self-contained (NumPy only) and reads like a small textbook:
every function carries the formula it implements and a one-line reference to the
place in Meucci, "Risk and Asset Allocation" (Springer, 2005) or the original
Ledoit-Wolf / James-Stein papers.

The three big ideas
-------------------
1. SAMPLE moments (mean, covariance) are unbiased but *high variance*: with N
   assets and T observations you must estimate N means and N(N+1)/2 covariance
   entries, and when N is not tiny relative to T the estimates are extremely
   noisy. Plugging them into a portfolio optimiser produces wild, unreliable
   weights (the "error-maximisation" of Michaud).

2. SHRINKAGE trades a little bias for a large variance reduction by pulling the
   noisy sample estimate toward a low-variance, structured TARGET:

        theta_shrunk = (1 - a) * theta_sample + a * theta_target,   a in [0, 1].

   The optimal intensity `a` can be derived analytically (James-Stein for the
   mean; Ledoit-Wolf for the covariance).

3. BAYESIAN reading: the shrinkage estimator IS the posterior mean when the
   target plays the role of the prior and `a` encodes the prior's confidence.
   That is the bridge the *_pymc notebook makes explicit.

Conventions
-----------
X : array (T, N)   T observations (rows) of N assets (columns), in return units.
mu    : (N,)       mean vector.
Sigma : (N, N)     covariance matrix.
We use the MLE covariance (divide by T) unless stated, to match Meucci.
"""

import numpy as np

# --------------------------------------------------------------------------- #
#  Sample (plug-in) moments                                                     #
# --------------------------------------------------------------------------- #

def sample_moments(X, mle=True):
    """Sample mean and covariance.

    mle=True divides the covariance by T (maximum-likelihood, Meucci's default);
    mle=False divides by T-1 (unbiased).
    """
    X = np.asarray(X, float)
    T = X.shape[0]
    mu = X.mean(axis=0)
    Xc = X - mu
    denom = T if mle else T - 1
    Sigma = Xc.T @ Xc / denom
    return mu, Sigma


# --------------------------------------------------------------------------- #
#  1. Shrinkage of the MEAN  (James-Stein / Meucci Ch.4)                         #
# --------------------------------------------------------------------------- #

def james_stein_mean(X, target=None):
    """Classic James-Stein shrinkage of the sample mean toward a target `b`.

    Stein's paradox (1956): for N >= 3 the sample mean is *inadmissible* -- the
    estimator that shrinks it toward any fixed point beats it in total
    mean-squared error, uniformly. The optimal intensity is

        a = min(1, (N - 2) / T * s2bar / ((mu_hat - b)' Sigma^-1 (mu_hat - b)))

    Here we use the version with the estimated covariance (Jorion 1986 style).
    `target` defaults to the grand mean (average across assets), a common choice.

    Returns (mu_shrunk, a, b).
    """
    X = np.asarray(X, float)
    T, N = X.shape
    mu_hat, Sigma_hat = sample_moments(X, mle=True)
    b = np.full(N, mu_hat.mean()) if target is None else np.asarray(target, float)

    diff = mu_hat - b
    Sinv = np.linalg.inv(Sigma_hat)
    quad = diff @ Sinv @ diff                       # (mu-b)' S^-1 (mu-b)
    a = (N - 2) / T / quad if quad > 0 else 0.0
    a = float(np.clip(a, 0.0, 1.0))
    mu_shr = (1 - a) * mu_hat + a * b
    return mu_shr, a, b


def meucci_shrink_location(X, target=None):
    """Meucci's location-shrinkage formula (S_ShrinkageEstimators.m).

        a = (1/T) * (sum(lam) - 2*max(lam)) / ((mu_hat - b)'(mu_hat - b))

    where lam are the eigenvalues of the sample covariance. This is the same
    James-Stein idea written with the eigenvalue spectrum; note it uses the
    plain (not Mahalanobis) distance to the target. Target defaults to 0.
    """
    X = np.asarray(X, float)
    T, N = X.shape
    mu_hat, Sigma_hat = sample_moments(X, mle=True)
    b = np.zeros(N) if target is None else np.asarray(target, float)

    lam = np.linalg.eigvalsh(Sigma_hat)
    diff = mu_hat - b
    denom = diff @ diff
    a = (lam.sum() - 2 * lam.max()) / T / denom if denom > 0 else 0.0
    a = float(np.clip(a, 0.0, 1.0))
    mu_shr = (1 - a) * mu_hat + a * b
    return mu_shr, a, b


# --------------------------------------------------------------------------- #
#  2. Shrinkage of the COVARIANCE                                               #
# --------------------------------------------------------------------------- #

def _scaled_identity_target(Sigma):
    """Target C = mean(eigenvalues) * I  -- a sphere with the average variance."""
    N = Sigma.shape[0]
    mu = np.trace(Sigma) / N
    return mu * np.eye(N)


def meucci_shrink_scatter(X):
    """Meucci's scatter-shrinkage toward the scaled identity (Ch.4).

        C = mean(lam) * I
        a = (1/T) * [ (1/T) sum_t trace( (x_t x_t' - S)^2 ) ] / trace( (S - C)^2 )

    with x_t the demeaned observations and S the MLE covariance. Returns
    (Sigma_shrunk, a, C).
    """
    X = np.asarray(X, float)
    T, N = X.shape
    mu_hat, S = sample_moments(X, mle=True)
    Xc = X - mu_hat
    C = _scaled_identity_target(S)

    num = 0.0
    for t in range(T):
        M = np.outer(Xc[t], Xc[t]) - S
        num += np.trace(M @ M) / T
    den = np.trace((S - C) @ (S - C))
    a = num / T / den if den > 0 else 0.0
    a = float(np.clip(a, 0.0, 1.0))
    Sigma_shr = (1 - a) * S + a * C
    return Sigma_shr, a, C


def ledoit_wolf_identity(X):
    """Ledoit & Wolf (2004) 'A well-conditioned estimator for large-dimensional
    covariance matrices' -- shrinkage of the sample covariance toward

        F = m * I,     m = trace(S)/N   (average variance),

    with the ANALYTIC optimal intensity (their Theorem, delta-hat):

        S       = MLE covariance
        m       = <S, I>/N            (mean variance)
        d2      = ||S - m I||_F^2 / N (dispersion of S around the sphere)
        b2bar   = (1/N) (1/T^2) sum_t ||x_t x_t' - S||_F^2   (capped at d2)
        a*      = b2bar / d2          (shrinkage intensity, in [0,1])
        Sigma*  = a* * m I + (1 - a*) * S

    Returns (Sigma_star, a_star, F).
    """
    X = np.asarray(X, float)
    T, N = X.shape
    mu_hat, S = sample_moments(X, mle=True)
    Xc = X - mu_hat

    m = np.trace(S) / N
    F = m * np.eye(N)
    d2 = np.sum((S - F) ** 2) / N                      # ||S - F||^2 / N

    b2 = 0.0
    for t in range(T):
        M = np.outer(Xc[t], Xc[t]) - S
        b2 += np.sum(M ** 2)
    b2 = b2 / (N * T * T)
    b2 = min(b2, d2)                                   # cap (Ledoit-Wolf)

    a = b2 / d2 if d2 > 0 else 0.0
    a = float(np.clip(a, 0.0, 1.0))
    Sigma_star = a * F + (1 - a) * S
    return Sigma_star, a, F


def constant_correlation_target(X):
    """Ledoit & Wolf (2003) 'Honey, I shrunk the sample covariance matrix'
    target: keep the sample variances, replace all pairwise correlations by
    their common average rbar.

        F_ii = S_ii ;  F_ij = rbar * sqrt(S_ii S_jj),  rbar = mean off-diag corr.

    Returns the target F only (intensity handled by `shrink_to_target`).
    """
    X = np.asarray(X, float)
    _, S = sample_moments(X, mle=True)
    s = np.sqrt(np.diag(S))
    R = S / np.outer(s, s)
    N = S.shape[0]
    off = R[np.triu_indices(N, 1)]
    rbar = off.mean()
    F = rbar * np.outer(s, s)
    np.fill_diagonal(F, np.diag(S))
    return F, rbar


def shrink_to_target(X, F):
    """Generic Ledoit-Wolf shrinkage of S toward an arbitrary target F, using the
    same b2/d2 intensity estimator as `ledoit_wolf_identity`. Lets you shrink
    toward the constant-correlation target (or any structured F)."""
    X = np.asarray(X, float)
    T, N = X.shape
    mu_hat, S = sample_moments(X, mle=True)
    Xc = X - mu_hat
    d2 = np.sum((S - F) ** 2) / N
    b2 = 0.0
    for t in range(T):
        M = np.outer(Xc[t], Xc[t]) - S
        b2 += np.sum(M ** 2)
    b2 = min(b2 / (N * T * T), d2)
    a = float(np.clip(b2 / d2 if d2 > 0 else 0.0, 0.0, 1.0))
    return a * F + (1 - a) * S, a


# --------------------------------------------------------------------------- #
#  3. Diagnostics: eigenvalue dispersion & conditioning                         #
# --------------------------------------------------------------------------- #

def eigenvalue_dispersion(Sigma):
    """Return sorted eigenvalues (descending) and a dispersion ratio
    max(lam)/min(lam) (= condition number). Sample covariances systematically
    OVER-disperse eigenvalues -- the largest are biased up, the smallest down --
    which is exactly what shrinkage toward a sphere corrects."""
    lam = np.linalg.eigvalsh(Sigma)[::-1]
    cond = lam[0] / lam[-1] if lam[-1] > 0 else np.inf
    return lam, cond


def condition_number(Sigma):
    lam = np.linalg.eigvalsh(Sigma)
    return lam[-1] / lam[0] if lam[0] > 0 else np.inf


# --------------------------------------------------------------------------- #
#  4. Bayesian bridge (used by the *_pymc notebook narrative)                    #
# --------------------------------------------------------------------------- #

def bayes_mean_posterior(mu_hat, Sigma, T, mu_0, tau2):
    """Posterior mean of a Normal-mean model with a Normal prior, showing that
    the Bayesian posterior mean is EXACTLY a shrinkage estimator.

    Likelihood : mu_hat | mu ~ N(mu, Sigma/T)
    Prior      : mu        ~ N(mu_0, tau2 * I)
    Posterior mean = (1-A) mu_hat + A mu_0  with matrix weight
        A = (Sigma/T) (Sigma/T + tau2 I)^-1.

    Returns the posterior mean vector.
    """
    N = len(mu_hat)
    SoT = Sigma / T
    A = SoT @ np.linalg.inv(SoT + tau2 * np.eye(N))
    return (np.eye(N) - A) @ mu_hat + A @ mu_0


# --------------------------------------------------------------------------- #
#  5. Simulation helper (synthetic ground-truth experiments)                    #
# --------------------------------------------------------------------------- #

def make_true_market(N, rho=0.6, vol_lo=0.10, vol_hi=0.40, mu_scale=0.5, seed=0):
    """Construct a known (mu_true, Sigma_true) 'market' in the style of Meucci's
    scripts: equicorrelation rho, volatilities ramped from vol_lo to vol_hi, and
    a mean proportional to Sigma*1 (so a mean-variance optimum exists). Use this
    to generate synthetic samples whose truth you know, so estimation error is
    measurable."""
    rng = np.random.default_rng(seed)
    vols = np.linspace(vol_lo, vol_hi, N)
    C = (1 - rho) * np.eye(N) + rho * np.ones((N, N))
    Sigma = np.outer(vols, vols) * C
    mu = mu_scale * Sigma @ np.ones(N) / N
    return mu, Sigma, rng


def pru_loss(Sigma_est, Sigma_true):
    """Percentage Relative improvement Utility-agnostic loss: squared Frobenius
    distance between estimate and truth, normalised by the truth's norm. Lower
    is better. A simple scalar to compare estimators in synthetic experiments."""
    return np.sum((Sigma_est - Sigma_true) ** 2) / np.sum(Sigma_true ** 2)

References

• Ledoit, O. & Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 88(2), 365–411. — the closed-form optimal shrinkage intensity derived and verified here
• Stein, C. (1956). Inadmissibility of the usual estimator for the mean of a multivariate normal distribution. Proceedings of the Third Berkeley Symposium. — the paradox demonstrated by Monte Carlo in the first notebook
• James, W. & Stein, C. (1961). Estimation with quadratic loss. Proceedings of the Fourth Berkeley Symposium. — the shrinkage estimator for the mean vector
• Meucci, A. (2005). Risk and Asset Allocation. Springer. — chapter 4, estimation and the shrinkage of location and scatter; the framework of this arc
• Michaud, R. O. (1989). The Markowitz optimization enigma: is optimized optimal? Financial Analysts Journal 45(1), 31–42. — the optimiser as error-maximiser, the motivation for the whole exercise
• DeMiguel, V., Garlappi, L. & Uppal, R. (2009). Optimal versus naive diversification. Review of Financial Studies 22(5), 1915–1953. — how badly estimation error degrades optimised portfolios in practice
• Fama, E. F. & French, K. R. — 48 industry portfolios, Kenneth French Data Library. — the canonical high-dimensional benchmark used for the closing backtest