""" 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 = /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)