Covariance estimation: shrinkage, factor models and denoising.

Too many parameters for so little data

Lesson 4.1.1 ended with a rule: variance is estimated with days; the mean, with years. This lesson deals with the part that can be estimated with days, the covariance, and with why not even that comes for free. With $N$ assets, the covariance matrix has $N(N+1)/2$ distinct parameters. With the 37 risky assets of the universe (all but the money market) that is 703. Six months of daily data provide 128 observations of each asset: the number of data points grows with $N \cdot T$, but the number of parameters, with $N^2$.

The starting estimator is the sample covariance,

unbiased element by element. The problem is not in each element, but in the matrix as a whole: when the ratio $q = N/T$ is not small, its eigenvalues spread out, the large ones come out larger and the small ones smaller than they are. An optimizer looking for minimum variance goes straight to those directions of underestimated variance.

There are four families of answers, and this lesson compares them all against the truth:

flowchart TB
  M[Sample] --> L[Shrinkage]
  M --> F[Factors]
  M --> D[RMT denoising]
  L --> L1[Ledoit-Wolf]
  F --> F1[Observed]
  F --> F2[PCA]
  D --> D1[Marchenko-Pastur]

The window and the truth

The data is that of the level 4 synthetic universe (lesson 4.1.1): 38 indices from 2006 to 2025 generated with a published ten-factor model. The true covariance is in universo4_cov_verdad.csv, one per regime. The working window is the second half of 2024, entirely in the calm regime, so the truth it is compared against is the covariance of the calm. The money market is left out: it is cash, and with a volatility of 0.2% it would swallow any minimum variance portfolio.

import numpy as np
import pandas as pd

precios = pd.read_parquet("universo4_precios.parquet")
precios = precios.set_index("fecha")
fac = pd.read_parquet("universo4_factores_diarios.parquet")
fac = fac.set_index("fecha")
r = precios.pct_change()
r.iloc[0] = precios.iloc[0] / 100 - 1
exceso = r.sub(fac["rf"], axis=0).drop(columns="MONET")
activos = list(exceso.columns)

cov_v = pd.read_csv("universo4_cov_verdad.csv")
calma = cov_v[cov_v["regimen"] == "calma"].set_index("activo")
S_v = calma.loc[activos, activos].to_numpy()

ventana = exceso.loc["2024-07-01":"2024-12-31"]
print("regimes in the window:",
      list(fac.loc[ventana.index, "regimen"].unique()))
X = ventana.to_numpy() * np.sqrt(252)     # annual scale
T, N = X.shape
print(f"N = {N} assets, T = {T} sessions, N/T = {N / T:.2f}")
print("covariance parameters:", N * (N + 1) // 2)

regimes in the window: ['calma']
N = 37 assets, T = 128 sessions, N/T = 0.29
covariance parameters: 703

Multiplying daily returns by $\sqrt{252}$ puts the covariances directly on the annual scale, the same as the published truth. The eigenvalues of the correlation matrix show where the problem is:

def a_corr(S):
    d = np.sqrt(np.diag(S))
    return S / np.outer(d, d)

S_m = np.cov(X, rowvar=False)
l_v = np.linalg.eigvalsh(a_corr(S_v))[::-1]
l_m = np.linalg.eigvalsh(a_corr(S_m))[::-1]
print("largest, truth:", l_v[:4].round(2),
      " sample:", l_m[:4].round(2))
print("smallest, truth:", l_v[-3:].round(4), " sample:",
      l_m[-3:].round(4))
print(f"condition number: truth {np.linalg.cond(S_v):,.0f}, "
      f"sample {np.linalg.cond(S_m):,.0f}")

largest, truth: [18.21  7.27  2.78  1.13]  sample: [17.08  7.6   2.9   1.16]
smallest, truth: [0.0202 0.013  0.0117]  sample: [0.0139 0.0129 0.0064]
condition number: truth 32,579, sample 31,703

The first eigenvalue, the market's, absorbs half of the trace (37 in a correlation of 37 assets). The small eigenvalues are very small in the truth too: the universe has hedged and unhedged versions of the same stock markets, four Treasury maturity buckets that load on the same rate factor, and sectors that are almost the index. Those combinations really have very little risk. In the sample, the smallest eigenvalue comes out at half the true one (0.0064 against 0.0117): the sample believes there is a portfolio with half the risk it really has.

How to measure an estimator

The distance between matrices (for example, the Frobenius norm of the difference of correlations) is the most cited measure and the least useful: it is dominated by the large elements and says nothing about what happens to a portfolio. This lesson uses two portfolio measures, with the minimum variance portfolio, which only depends on the covariance:

• Real volatility against the optimum: the weights are computed with the estimate and their volatility is evaluated with the truth, divided by that of the true minimum variance portfolio. 1 is the best possible.
• Promised volatility against the real one: the volatility the estimate attributes to its own weights, divided by the true one. Below 1, the portfolio has more risk than it believes.

uno = np.ones(N)
def min_var(S):
    x = np.linalg.solve(S, uno)
    return x / x.sum()

w_opt = min_var(S_v)
vol_opt = np.sqrt(w_opt @ S_v @ w_opt)

def evalua(S):
    w = min_var(S)
    real = np.sqrt(w @ S_v @ w)
    return {"err_corr": np.linalg.norm(a_corr(S) - a_corr(S_v))
            / np.linalg.norm(a_corr(S_v)),
            "vol_real/optima": real / vol_opt,
            "prometida/real": np.sqrt(w @ S @ w) / real}

print(f"vol. of the true minimum variance: {vol_opt:.2%}")
print({k: round(float(v), 3) for k, v in evalua(S_m).items()})

vol. of the true minimum variance: 0.47%
{'err_corr': 0.173, 'vol_real/optima': 1.28, 'prometida/real': 0.657}

The true minimum variance portfolio has a volatility of 0.47% (short bonds and hedged combinations). The one built with the sample covariance has 28% more volatility, and promises 66% of what it has: it believes it has a third less risk. Part of that difference is scale (the half-year was quieter than the average of the calm regime); the Monte Carlo further down, without that effect, gives 72% with 126 sessions. It is the bias of lesson 4.1.1 applied to the covariance: the optimizer picks the directions with the most favorable error.

Shrinkage: Ledoit-Wolf

The idea of shrinkage is to combine the sample covariance, unbiased but noisy, with a structured target, biased but stable, $F$:

Ledoit and Wolf (2004) chose the scaled identity as target, $F = m I$ with $m = \operatorname{tr}(S)/N$, and derived the intensity $\delta$ that minimizes the expected squared error. Their estimator of the intensity compares the dispersion of the sample covariance around the target with the variance of the sample covariance itself:

If the sample covariance is very noisy ($\overline{b}^2$ large relative to $d^2$), it shrinks a lot; with a lot of data, little. It is short to implement and can be checked against scikit-learn:

from sklearn.covariance import LedoitWolf

def ledoit_wolf(X):
    Xc = X - X.mean(axis=0)
    T, N = Xc.shape
    S = Xc.T @ Xc / T
    m = np.trace(S) / N
    d2 = np.sum((S - m * np.eye(N)) ** 2) / N
    b2 = sum(np.sum((np.outer(x, x) - S) ** 2) for x in Xc)
    b2 = min(b2 / N / T**2, d2)
    delta = b2 / d2
    return delta * m * np.eye(N) + (1 - delta) * S, delta

S_lw, delta = ledoit_wolf(X)
skl = LedoitWolf().fit(X)
print(f"intensity: custom {delta:.4f}, sklearn {skl.shrinkage_:.4f}")
print("same matrix:", np.allclose(S_lw, skl.covariance_, atol=1e-12))

intensity: custom 0.0490, sklearn 0.0490
same matrix: True

The custom implementation and scikit-learn's match (both use the covariance with divisor $T$, not $T-1$). The intensity is small, 4.9%, because 128 sessions are not that few for 37 assets. And even so the result is disastrous, because of the target:

def lw_correlacion(X):
    sd = X.std(axis=0)
    C, delta = ledoit_wolf((X - X.mean(axis=0)) / sd)
    return a_corr(C) * np.outer(sd, sd), delta

S_lwc, delta_c = lw_correlacion(X)

from pypfopt import risk_models
S_cc = risk_models.CovarianceShrinkage(
    pd.DataFrame(X), returns_data=True, frequency=1
).ledoit_wolf(shrinkage_target="constant_correlation").to_numpy()

tabla = pd.DataFrame({"sample": evalua(S_m),
                      "LW covariance": evalua(S_lw),
                      "LW correlation": evalua(S_lwc),
                      "LW constant corr.": evalua(S_cc)}).T
print(tabla.round(3))
vols = np.sqrt(np.diag(S_m))
print("min., mean and max. vol.:", vols.min().round(3),
      vols.mean().round(3), vols.max().round(3))

                   err_corr  vol_real/optima  prometida/real
sample                0.173            1.280           0.657
LW covariance         0.226            2.758           1.333
LW correlation        0.190            1.238           1.077
LW constant corr.     0.205            1.222           1.150
min., mean and max. vol.: 0.014 0.112 0.224

The volatilities of the assets range from 1.4% (1-3 year bond) to 22% (energy). The scaled identity pushes all variances toward their mean: the short bond, the key asset of the minimum variance portfolio, looks riskier than it is, and the resulting portfolio has 2.76 times the optimal volatility. Shrinkage is not a “better” estimator in the abstract: it is as good as its target.

There are two ways to fix it. The first, shrink the correlation: standardize the returns, apply Ledoit-Wolf to that matrix (whose natural target is indeed the identity) and put the sample volatilities back. The second, the constant correlation target of Ledoit and Wolf (2004, “Honey, I Shrunk...”): all correlations equal to the mean, each variance its own. PyPortfolioOpt ships it ready-made; the argument frequency=1 stops it from multiplying by 252, because the data is already on the annual scale. Both give about 1.23 times the optimal volatility, somewhat better than the sample covariance (1.28), and promise a little more risk than they have (1.08 and 1.15): they err on the prudent side.

Factor models

A factor model assumes that returns depend on a few common factors and on a residual specific to each asset, independent of the others':

with $D$ diagonal. Instead of 703 parameters there are $N \cdot k$ loadings, $k(k+1)/2$ factor covariances and $N$ residual variances. There are two ways to get the factors:

• Observed: known series (the market, rates, the dollar, the Fama-French factors of lesson 4.1.6). The loadings come from asset-by-asset regressions.
• Statistical: the first $k$ principal components of the sample itself. You do not need to know what they are, but you have to choose $k$.

The universe provides the realized returns of its ten factors in universo4_factores_diarios.parquet:

F = fac.loc[ventana.index,
            [c for c in fac.columns if c.startswith("f_")]]
F = F.to_numpy() * np.sqrt(252)

def factorial(X, F):
    Fc = np.column_stack([np.ones(len(F)), F])
    coef, *_ = np.linalg.lstsq(Fc, X, rcond=None)
    B = coef[1:]
    resid = X - Fc @ coef
    D = np.diag(resid.var(axis=0, ddof=Fc.shape[1]))
    return B.T @ np.cov(F, rowvar=False) @ B + D

def pca(X, k):
    S = np.cov(X, rowvar=False)
    l, V = np.linalg.eigh(S)
    l, V = l[::-1][:k], V[:, ::-1][:, :k]
    Sk = V @ np.diag(l) @ V.T
    return Sk + np.diag(np.diag(S - Sk))

S_f = factorial(X, F)
tabla.loc["observed factors"] = evalua(S_f)
tabla.loc["PCA, 5 components"] = evalua(pca(X, 5))
print(tabla.round(3).iloc[[0, 2, 4, 5]])

                   err_corr  vol_real/optima  prometida/real
sample                0.173            1.280           0.657
LW correlation        0.190            1.238           1.077
observed factors      0.173            1.080           0.982
PCA, 5 components     0.184            1.267           1.113

The observed factor model is the best of the lesson: 1.08 times the optimal volatility and an almost exact promised risk (0.98). It has to be read carefully: the data was generated with exactly these ten factors and independent residuals, so the model is well specified by construction. In real markets the factors are approximations and the residuals have correlations (sector, style) that the model ignores. The model with five principal components illustrates it: it leaves out weak factors, such as the dollar (which separates the hedged versions from the unhedged ones), its residual stops being diagonal and the portfolio ends up at 1.27, like the sample covariance.

Denoising with Marchenko-Pastur

Random matrix theory provides a reference for separating signal from noise. If the $N$ assets were independent with variance $\sigma^2$, the eigenvalues of their sample correlation with $T$ observations would not all come out equal to 1: when $N$ and $T$ are large, they would be distributed according to the Marchenko-Pastur distribution, over the interval

Laloux, Cizeau, Bouchaud and Potters (1999) observed that almost all the eigenvalues of the correlation of S&P 500 stocks fell inside that interval: indistinguishable from noise. The denoising popularized by López de Prado (2020) acts accordingly: it keeps the eigenvalues above $\lambda_+$, replaces the rest by their mean (to preserve the trace) and rebuilds the correlation.

flowchart TB
  S[Covariance] --> C[Correlation]
  C --> A[Eigenvalues]
  A --> R{"Greater than λ+?"}
  R -- yes --> K[Kept]
  R -- no --> P[Averaged]
  K --> N[Rebuild]
  P --> N
  N --> V[Volatilities]

q = N / T
lam_mas, lam_menos = (1 + np.sqrt(q)) ** 2, (1 - np.sqrt(q)) ** 2
print(f"Marchenko-Pastur with σ² = 1: [{lam_menos:.3f}, {lam_mas:.3f}]")
print("eigenvalues above λ+: sample", (l_m > lam_mas).sum(),
      "| truth", (l_v > lam_mas).sum())
print("true eigenvalues below λ-:",
      (l_v < lam_menos).sum(), "of", N)

Marchenko-Pastur with σ² = 1: [0.214, 2.364]
eigenvalues above λ+: sample 3 | truth 3
true eigenvalues below λ-: 21 of 37

With six months, the upper edge is 2.36 and only three sample eigenvalues exceed it; in the truth, also three. So far, the recipe works. But the truth is not “three factors and white noise”: 21 of the 37 true eigenvalues lie below the lower edge. They are not noise: they are the very low risk combinations the universe has (hedged against unhedged, bond maturity buckets). Averaging them with the rest raises their variance, which is the same as erasing those hedges:

def denoise(X, ajusta=True):
    S = np.cov(X, rowvar=False)
    sd = np.sqrt(np.diag(S))
    l, V = np.linalg.eigh(a_corr(S))
    q = X.shape[1] / X.shape[0]
    s2 = 1.0
    for _ in range(50 if ajusta else 0):   # fixed point for σ²
        s2 = l[l <= s2 * (1 + np.sqrt(q)) ** 2].mean()
    ruido = l <= s2 * (1 + np.sqrt(q)) ** 2
    l = l.copy()
    l[ruido] = l[ruido].mean()
    C = a_corr(V @ np.diag(l) @ V.T)
    return C * np.outer(sd, sd), int((~ruido).sum()), s2

S_dn, k1, _ = denoise(X, ajusta=False)
S_dn2, k2, s2 = denoise(X)
print(f"σ² = 1: keeps {k1}; fitted σ² = {s2:.3f}: keeps {k2}")
tabla.loc["denoising σ² = 1"] = evalua(S_dn)
tabla.loc["denoising fitted σ²"] = evalua(S_dn2)
print(tabla.round(3).iloc[[0, 2, 6, 7]])

σ² = 1: keeps 3; fitted σ² = 0.108: keeps 14
                     err_corr  vol_real/optima  prometida/real
sample                  0.173            1.280           0.657
LW correlation          0.190            1.238           1.077
denoising σ² = 1        0.209            1.504           1.239
denoising fitted σ²     0.172            1.291           1.010

With $\sigma^2 = 1$, denoising keeps three eigenvalues, flattens the other 34 and the minimum variance portfolio worsens to 1.50 times the optimum, worse than the sample covariance. Fitting $\sigma^2$ as López de Prado does (here with a simple fixed point: the noise variance is the mean of the eigenvalues that stay inside the edge), the edge drops, 14 eigenvalues are kept and the result matches the sample covariance (1.29), with a much more honest promised risk (1.01 against 0.66). It is a useful tool when the true correlation resembles the assumption (a few strong factors and a residual that looks like noise); when there is structure in the small eigenvalues, you have to check before applying it.

Many windows: Monte Carlo from the truth

A single window is a single sample. To compare with more data, returns are simulated with the true model (loadings from universo4_verdad.csv, calm factor covariance from universo4_factores_verdad.csv and independent residuals), with normals and without GARCH to isolate estimation error, and the comparison is repeated with 30 samples per window length:

verdad = pd.read_csv("universo4_verdad.csv", index_col="activo")
fv = pd.read_csv("universo4_factores_verdad.csv")
fv = fv[fv["regimen"] == "calma"].set_index("factor")
nombres = list(fv.index)
B = verdad.loc[activos, [f"b_{k}" for k in nombres]].to_numpy()
vol_f = fv["vol_anual"].to_numpy()
S_fv = fv[[f"corr_{k}" for k in nombres]].to_numpy()
S_fv = S_fv * np.outer(vol_f, vol_f)
idio = verdad.loc[activos, "vol_idio"].to_numpy()
print("B·Σf·B' + D = calm:",
      np.allclose(B @ S_fv @ B.T + np.diag(idio**2), S_v))

B·Σf·B' + D = calm: True

rng = np.random.default_rng(7)
Lf = np.linalg.cholesky(S_fv)
estimadores = {
    "sample": lambda X, F: np.cov(X, rowvar=False),
    "LW correlation": lambda X, F: lw_correlacion(X)[0],
    "factors": factorial,
    "denoising σ²=1": lambda X, F: denoise(X, ajusta=False)[0],
    "denoising fitted": lambda X, F: denoise(X)[0],
}
filas = []
for T_ in (63, 126, 252, 504):
    for _ in range(30):
        Fs = rng.standard_normal((T_, len(nombres))) @ Lf.T
        Xs = Fs @ B.T + rng.standard_normal((T_, N)) * idio
        for nombre, est in estimadores.items():
            m = evalua(est(Xs, Fs))
            filas.append({"T": T_, "estimador": nombre, **m})
mc = pd.DataFrame(filas).groupby(["T", "estimador"]).mean()
print(mc["vol_real/optima"].unstack().round(2))
print(mc["prometida/real"].unstack().round(2))

estimador  LW correlation  denoising fitted  denoising σ²=1  factors  sample
T
63                   1.26              1.27            1.96     1.09    1.55
126                  1.13              1.19            1.57     1.04    1.17
252                  1.08              1.16            1.46     1.02    1.08
504                  1.04              1.05            1.44     1.01    1.04
estimador  LW correlation  denoising fitted  denoising σ²=1  factors  sample
T
63                   1.10              0.90            1.02     0.91    0.43
126                  1.06              0.84            1.17     0.96    0.72
252                  1.05              0.76            1.27     0.98    0.87
504                  1.02              0.93            1.30     0.99    0.92

The figure at the top summarizes the first table. With three months of data (63 sessions, $N/T = 0.59$), the sample minimum variance portfolio has 55% more volatility than the optimum and promises 43% of what it has. With two years, the difference drops to 4%. Ledoit-Wolf on the correlation wins above all with short windows (1.26 against 1.55 at 63 sessions) and converges with the sample covariance when there is data: its intensity falls on its own. The correctly specified factor model is the best at every length. And denoising with $\sigma^2 = 1$ is the worst with long windows: with 504 sessions the edge narrows, but it still erases the small directions that are real, so its error does not disappear (1.44). Fitted denoising, on the other hand, systematically promises less risk than the real one (from 0.76 to 0.93).

What no estimator sees

All the estimators of this lesson estimate the covariance of the regime in front of them. The 2024 window is calm; the published truth also includes the stress covariance. It is enough to evaluate the same portfolios with it:

estres = cov_v[cov_v["regimen"] == "estres"].set_index("activo")
S_e = estres.loc[activos, activos].to_numpy()
for nombre, S in [("sample", S_m), ("LW correlation", S_lwc),
                  ("factors", S_f)]:
    w = min_var(S)
    print(f"{nombre:15s} promised {np.sqrt(w @ S @ w):.2%}  "
          f"calm {np.sqrt(w @ S_v @ w):.2%}  "
          f"stress {np.sqrt(w @ S_e @ w):.2%}")

sample          promised 0.39%  calm 0.60%  stress 0.81%
LW correlation  promised 0.62%  calm 0.58%  stress 0.78%
factors         promised 0.50%  calm 0.50%  stress 0.70%

Under stress, the volatility of the three portfolios rises by 34% to 40% over that of the calm, and for the sample covariance it is more than double what was promised. It is a different problem from estimation error: it is not solved with more calm data, but with stress scenarios and risk limits (lessons 4.2.4 and 4.2.5).

Which estimator to use

The figure brings together the window of the lesson: the real volatility of the minimum variance portfolio built with each estimator, in multiples of the optimum.

Typical mistakes:

• Comparing estimators with the Frobenius norm. In the lesson window, the correlation error ranges from 0.17 to 0.23 across the eight estimators and does not rank the portfolio results correctly.
• Non-synchronous data: Tokyo, Frankfurt and New York closes on the same day underestimate the correlations between regions. With daily data from markets in different time zones, use weekly returns or correct for the lags.
• Mixing scales: daily and annual, simple and logarithmic, or one asset in local currency and another in USD.
• Windows chosen without thinking about the regime: a window that mixes calm and crisis estimates a covariance that is neither that of one nor that of the other.

Exercises

1. Repeat the estimator table with two years of calm data (2023-01-01 to 2024-12-31). Which estimator wins now and how much does the difference with the sample covariance shrink?

v2 = exceso.loc["2023-01-01":"2024-12-31"]
X2 = v2.to_numpy() * np.sqrt(252)
F2 = fac.loc[v2.index,
             [c for c in fac.columns if c.startswith("f_")]]
F2 = F2.to_numpy() * np.sqrt(252)
print("T =", len(X2), "| regimes:",
      list(fac.loc[v2.index, "regimen"].unique()))
t2 = pd.DataFrame({
    "sample": evalua(np.cov(X2, rowvar=False)),
    "LW correlation": evalua(lw_correlacion(X2)[0]),
    "factors": evalua(factorial(X2, F2)),
    "denoising fitted": evalua(denoise(X2)[0])}).T
print(t2.round(3))

T = 502 | regimes: ['calma']
                  err_corr  vol_real/optima  prometida/real
sample               0.113            1.053           0.847
LW correlation       0.120            1.068           0.984
factors              0.113            1.025           0.913
denoising fitted     0.113            1.058           0.861

2. A manager rebalances to the minimum variance portfolio every six months. Compute the turnover (half the sum of the absolute weight changes) between the first and second halves of 2024 with the sample covariance and with Ledoit-Wolf on the correlation.

def pesos_mv(desde, hasta, est):
    Xv = exceso.loc[desde:hasta].to_numpy() * np.sqrt(252)
    return min_var(est(Xv))

muestral = lambda Xv: np.cov(Xv, rowvar=False)
lwc = lambda Xv: lw_correlacion(Xv)[0]
for nombre, est in [("sample", muestral), ("LW correlation", lwc)]:
    w1 = pesos_mv("2024-01-01", "2024-06-30", est)
    w2 = pesos_mv("2024-07-01", "2024-12-31", est)
    print(f"{nombre:15s} turnover {0.5 * np.abs(w2 - w1).sum():.0%}  "
          f"gross {np.abs(w2).sum():.0%}")

sample          turnover 60%  gross 196%
LW correlation  turnover 28%  gross 160%

3. Explain in your own words why a very small eigenvalue of the true correlation, which corresponds to “hedged version minus unhedged version of Europe”, is not noise, and what a denoising procedure would have to do to respect it.