--- title: "4.1.2 Covariances: shrinkage, factors and denoising" description: "Sample covariance, Ledoit-Wolf, factor models and Marchenko-Pastur denoising, measured against the published true covariance." url: "https://pythonparatrading.com/en/python-trading-course/level-4/covariance-estimation-shrinkage-and-denoising/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-4/covariance-estimation-shrinkage-and-denoising.png" --- # 4.1.2 · Covariance estimation: shrinkage, factor models and denoising > With 37 assets, a covariance has 703 parameters, and six months of data give 128 observations per asset. This lesson compares four families of estimators (sample, Ledoit-Wolf, factor models and random matrix theory denoising) against the true covariance of the synthetic universe, and measures what matters to a portfolio: how much risk it really has and how much it thinks it has. - Course: [Level 4: Portfolios, risk, execution and live trading](https://pythonparatrading.com/en/python-trading-course/level-4/) - Module: 4.1 · Portfolio theory and optimisation - Duration: 70 min - Prerequisites: [4.1.1 · Markowitz: efficient frontier and real constraints](https://pythonparatrading.com/en/python-trading-course/level-4/markowitz-efficient-frontier-and-constraints/) ## By the end of this lesson - Explain why the sample covariance degrades when the number of assets approaches the number of observations and how this shows up in the eigenvalues. - Implement the Ledoit-Wolf estimator, check it against scikit-learn and choose the right shrinkage target (identity, correlation, constant correlation). - Build a covariance with an observed factor model and with principal components, and know when each one is well specified. - Apply Marchenko-Pastur denoising and recognize when it erases real structure instead of noise. - Measure the estimators with portfolio criteria: real volatility of the minimum variance portfolio and ratio between promised and real risk. ## 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, $$ S = \frac{1}{T-1} \sum_{t=1}^{T} (x_t - \overline{x})(x_t - \overline{x})^\top , $$ 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: ```mermaid 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. ```python 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) ``` ```output 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: ```python 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}") ``` ```output 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. ```python 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()}) ``` ```output 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$: $$ \hat\Sigma = \delta F + (1 - \delta) S, \qquad 0 \le \delta \le 1 . $$ 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: $$ d^2 = \tfrac{1}{N}\lVert S - mI \rVert_F^2, \qquad \overline{b}^2 = \tfrac{1}{N T^2} \sum_{t} \lVert x_t x_t^\top - S \rVert_F^2, \qquad \delta = \frac{\min(\overline{b}^2, d^2)}{d^2} . $$ 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: ```python 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)) ``` ```output 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: ```python 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)) ``` ```output 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': $$ x_t = B f_t + \varepsilon_t, \qquad \Sigma = B\, \Sigma_f\, B^\top + D , $$ 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`: ```python 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]]) ``` ```output 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 $$ \lambda_{\pm} = \sigma^2 \left(1 \pm \sqrt{N/T}\right)^2 . $$ 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. ```mermaid 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] ``` ```python 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) ``` ```output Marchenko-Pastur with σ² = 1: [0.214, 2.364] eigenvalues above λ+: sample 3 | truth 3 true eigenvalues below λ-: 21 of 37 ``` *Fig. — Eigenvalues of the 37-asset correlation, true and from six months of data, with the Marchenko-Pastur interval for pure noise.* Eigenvalues sorted from largest to smallest on a log scale. The first three, true and sample, exceed the upper Marchenko-Pastur edge (2.36, highlighted); the rest fall to 0.01, and 21 of the 37 true ones lie below the lower edge (0.214): real low-risk structure, not noise. The sample eigenvalues follow the true ones, with the smallest below them. 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: ```python 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]]) ``` ```output σ² = 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: ```python 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)) ``` ```output B·Σf·B' + D = calm: True ``` ```python 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)) ``` ```output 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: ```python 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%}") ``` ```output 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. *Fig. — Actual volatility of the minimum-variance portfolio built with each estimator in the second half of 2024, as a multiple of the optimum.* Horizontal bars. Sample 1.28; Ledoit-Wolf on the covariance 2.76 (highlighted: the scaled-identity target destroys the structure); Ledoit-Wolf on the correlation 1.24; constant correlation 1.22; observed factors 1.08; PCA with five components 1.27; denoising with σ² = 1 1.50; with fitted σ² 1.29. | Estimator | When it works | Main risk | |---|---|---| | Sample | $T$ much larger than $N$ (more than 10 times) | Promises less risk than the real one when $N/T$ grows | | Ledoit-Wolf on the correlation | By default, above all with short windows | Little, if the target is reasonable | | Ledoit-Wolf on the covariance | Assets with similar volatilities | Destroys the structure if they are not | | Observed factors | When the factors explain the structure | Misspecification (missing factors) | | PCA | A few strong factors | Choosing $k$; leaves out the weak ones | | Denoising (Marchenko-Pastur) | A few strong eigenvalues and a rest that looks like noise | Erases real hedges if there are very small eigenvalues | 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? **Solution.** ```python 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)) ``` ```output 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 ``` With 502 sessions, the sample covariance ends up 5% from the optimum and the estimators converge: Ledoit-Wolf barely shrinks and the factor model is still the best because it is well specified. The promised risk comes out below the real one for the sample covariance, the factor model and denoising: in those two years the joint volatility was 14% lower than the calm average, a scale error that no estimator corrects. Ledoit-Wolf, which in the Monte Carlo promises somewhat more risk than the real one (it shrinks the correlations and, with them, the hedges), offsets it here by chance (0.98). 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. **Solution.** ```python 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%}") ``` ```output sample turnover 60% gross 196% LW correlation turnover 28% gross 160% ``` The unconstrained minimum variance portfolio is a leveraged long-short portfolio, and its weights move a lot from one half-year to the next with either estimator. Ledoit-Wolf cuts the turnover by more than half (from 60% to 28%) and the gross leverage (from 196% to 160%), but the portfolio is still long-short and leveraged: shrinkage stabilizes the covariance, not the portfolio's decisions. For that you need constraints (lesson 4.1.1) and a turnover penalty (lesson 4.1.5). 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. **Solution.** The two versions of Europe have the same exposure to the European stock market and only differ in the dollar: their difference is, almost exactly, a currency position, with much less volatility than each of them separately. That low risk combination exists because of how the product is built, not because of sampling chance, and it produces a small eigenvalue of the correlation. Marchenko-Pastur describes pure noise matrices, with no directions like that; when the truth has them, the lower edge does not separate noise from signal. A procedure that respects them has to fit the noise level (σ²) instead of assuming 1, or act only on the intermediate eigenvalues, or start from a factor model that includes the factor that generates them (the dollar). ## Self-assessment 1. With six months of data, the minimum variance portfolio built with the sample covariance promises a volatility that is 66% of the real one. Why is it always wrong in that direction? - A. Because of a rounding error - B. Because the data has fat tails - C. Because the optimizer looks for the directions of lowest estimated variance, and those are precisely the ones the sample underestimates: the small eigenvalues of the sample covariance come out smaller than the true ones - D. Because the sample covariance is biased upward *Answer C: Because the optimizer looks for the directions of lowest estimated variance, and those are precisely the ones the sample underestimates: the small eigenvalues of the sample covariance come out smaller than the true ones.* The sample covariance is unbiased element by element, but its eigenvalues are not: the extremes spread out (the large ones come out larger and the small ones, smaller). Minimizing the estimated variance selects the directions with the most favorable error, and the promised risk ends up below the real one. 2. Ledoit-Wolf with the scaled identity target, applied to the covariance of 37 assets with volatilities from 1.4% to 22%, gives a minimum variance portfolio with 2.76 times the optimal volatility. What goes wrong? - A. The target: the scaled identity pushes all variances toward the average variance, so the short bond looks much riskier than it is and the portfolio stops using it - B. The intensity formula - C. scikit-learn has a bug - D. More data is needed *Answer A: The target: the scaled identity pushes all variances toward the average variance, so the short bond looks much riskier than it is and the portfolio stops using it.* The target has to be reasonable for the structure of the data. With such different volatilities, you shrink the correlation (with the sample variances) or use the constant correlation target; both give about 1.23. 3. In the course universe, 21 of the 37 eigenvalues of the true correlation lie below the lower Marchenko-Pastur edge with six months of data. What does that imply for denoising? - A. That denoising works better than ever - B. That more principal components should be used - C. That the true correlation is miscalculated - D. That the real structure is not “a few factors plus white noise”: there are very low risk directions that are real (hedged and unhedged versions, bonds of different maturities), and replacing all the small eigenvalues by their mean erases them *Answer D: That the real structure is not “a few factors plus white noise”: there are very low risk directions that are real (hedged and unhedged versions, bonds of different maturities), and replacing all the small eigenvalues by their mean erases them.* Marchenko-Pastur describes the eigenvalues of a pure noise matrix. If the truth has very small eigenvalues, as here, clipping everything that falls inside the edge flattens real directions. With σ² = 1 the minimum variance portfolio worsens to 1.50 times the optimum; fitting σ², it ends up like the sample covariance. 4. The ten observed factor model is the best estimator of the lesson (1.08 times the optimal volatility with six months). Why should you be careful about generalizing it? - A. Because factor models always overfit - B. Because here the data was generated with exactly those ten factors: the model is well specified by construction, which does not happen in real markets - C. Because it has more parameters than the sample covariance - D. Because it is not investable *Answer B: Because here the data was generated with exactly those ten factors: the model is well specified by construction, which does not happen in real markets.* A factor model wins when its factors capture the true structure. With wrong or incomplete factors (such as PCA with five components, which leaves out weak factors), the residual is not diagonal and the estimator loses its advantage: 1.27. 5. Which estimator anticipates the covariance of a crisis from six calm months? - A. Ledoit-Wolf - B. The factor model - C. None: they all estimate the covariance of the calm; under stress, the real volatility of the same portfolios rises by 34% to 40% - D. Denoising *Answer C: None: they all estimate the covariance of the calm; under stress, the real volatility of the same portfolios rises by 34% to 40%.* Estimators reduce estimation error within a regime. Regime change is a different problem: it has to be handled with stress scenarios and limits (lessons 4.2.4 and 4.2.5). ## Key points - The sample covariance has N(N+1)/2 parameters. When N/T is not small, its extreme eigenvalues spread out and the minimum variance portfolio promises less risk than it has: 66% of the real one with 37 assets and 128 sessions. - Ledoit-Wolf combines the sample covariance with a structured target, using the intensity that minimizes expected error. The custom implementation matches scikit-learn; the target matters: with heterogeneous volatilities, you shrink the correlation, not the covariance. - A well-specified factor model is the best estimator (1.08 times the optimal volatility with six months); with incomplete factors, such as five principal components, it loses the advantage. - Marchenko-Pastur denoising assumes the structure lies in a few large eigenvalues. If the truth has very low risk directions, it erases them: with σ² = 1 it worsens the portfolio; with fitted σ², it matches the sample covariance. - Estimators are judged by the portfolio, not by the distance between matrices: the Frobenius error barely tells them apart. And none of them anticipates the covariance of a crisis. ## Resources ### Research - [Price Discovery & Arbitrage Price Theory](https://pythonparatrading.com/en/arbitrage-pricing-theory-and-price-discovery-2/) · Research-0006 · 6 min ### Documentation - [scikit-learn: LedoitWolf](https://scikit-learn.org/stable/modules/generated/sklearn.covariance.LedoitWolf.html) · scikit-learn.org - [PyPortfolioOpt: risk models and shrinkage](https://pyportfolioopt.readthedocs.io/en/latest/RiskModels.html) · pyportfolioopt.readthedocs.io ### Books - [Ledoit, O. and Wolf, M. (2004): A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 88(2), 365-411](https://doi.org/10.1016/S0047-259X(03)00096-4) · doi.org - [Ledoit, O. and Wolf, M. (2004): Honey, I Shrunk the Sample Covariance Matrix. The Journal of Portfolio Management 30(4), 110-119](https://doi.org/10.3905/jpm.2004.110) · doi.org - [Laloux, L., Cizeau, P., Bouchaud, J.-P. and Potters, M. (1999): Noise Dressing of Financial Correlation Matrices. Physical Review Letters 83(7), 1467-1470](https://doi.org/10.1103/PhysRevLett.83.1467) · doi.org - [López de Prado, M. (2020): Machine Learning for Asset Managers. Cambridge University Press (chapter 2, denoising)](https://doi.org/10.1017/9781108883658) · doi.org ### Data - [Level 4 data: true covariance of the universe by regime (CSV)](https://pythonparatrading.com/curso/datos/nivel-4/universo4_cov_verdad.csv) · Download · universo4_cov_verdad.csv Previous lesson: [4.1.1 · Markowitz: efficient frontier and real constraints](https://pythonparatrading.com/en/python-trading-course/level-4/markowitz-efficient-frontier-and-constraints/) Next lesson: [4.1.3 · Black-Litterman: from market equilibrium to views](https://pythonparatrading.com/en/python-trading-course/level-4/black-litterman-model/)