The full notebook behind the video: reproducing Marcos Lopez de Prado's Hierarchical Risk Parity with a Monte Carlo experiment, and comparing it against Markowitz minimum variance and the inverse variance portfolio. Part 1 of 2 (real data comes in Part 2).
The paper on his own simulated data, his exact inputs (distance $\sqrt{\tfrac12(1-\rho)}$, single linkage, quasi-diagonalization, inverse-variance bisection). CLA is replaced by a long-only min-variance QP (same objective). The essential exhibits are reproduced below, and each step names the exact function that produces the result.
These are risk-based allocators, not return-seeking strategies. IVP, HRP and min-variance/CLA all take only the covariance matrix as input, with no expected-return forecast. The goal is to minimise risk and diversify, not to pick winners, which is exactly why they are judged on out-of-sample variance, not return.
Everything below compares three ways to turn a covariance matrix into portfolio weights. Keep this table handy, most exhibits are really about the "Inverts Σ?" and "Out-of-sample variance" rows.
| IVP, Inverse-Variance | HRP, Hierarchical Risk Parity | CLA / Min-Variance | |
|---|---|---|---|
| Weight rule | $w_i \propto 1/\sigma_i^2$ directly | inverse-variance between clusters, down the tree | minimize $w^{\top}\Sigma\,w$ (global) |
| Uses correlations? | No, only the diagonal (variances) | Yes, via the clustering / hierarchy | Yes, the full covariance matrix |
| Inverts Σ? | No | No | Yes, the fragile step |
| How it's built | flat, one shot | cluster → quasi-diagonalize → recursive bisection | constrained quadratic optimization |
| Diversification | by variance only; double-counts redundant assets | spread across the correlation structure | concentrates, zeroes assets (~93% top-5) |
| Numerical stability | very stable | stable (never inverts) | breaks when Σ is ill-conditioned (Markowitz's curse) |
| Estimation-error sensitivity | low | low | high, amplified by the condition number κ |
| Out-of-sample variance (this run, 2,000 sims) | 0.098 (middle) | 0.072, lowest | 0.122 (highest) |
| Which function | getIVP |
getHRP (getRecBipart) |
getMV |
| In the paper | benchmark | the proposed method | CLA (Critical Line Algorithm) |
| Nickname | "naive risk parity" | "IVP organized by the hierarchy" | "the optimizer" |
The three stages plus two helpers, this is the whole method:
correlDist(corr), turns the correlation matrix into a proper distance $d=\sqrt{\tfrac12(1-\rho)}$: correlated assets are close, uncorrelated are far.getQuasiDiag(link), Stage 2: reads the clustering tree (link) and returns the asset order that puts similar assets next to each other.getClusterVar / getIVP, the inverse-variance math used inside a cluster (getIVP = weights ∝ 1/variance).getRecBipart(cov, sortIx), Stage 3: walks the tree top-down and splits weight between the two halves inversely to their variance → the final HRP weights.getHRP(cov, corr), one call that chains all three: correlDist → sch.linkage(...,'single') → getQuasiDiag → getRecBipart.getMV(cov), long-only minimum-variance (our stand-in for the paper's CLA; identical objective).import random, numpy as np, pandas as pd, matplotlib.pyplot as plt
import scipy.cluster.hierarchy as sch
from scipy.cluster.hierarchy import ClusterWarning
from scipy.optimize import minimize
from scipy.linalg import block_diag
import networkx as nx, warnings
warnings.filterwarnings("ignore", category=ClusterWarning)
plt.rcParams["figure.figsize"] = (7, 4.5)
def getIVP(cov, **kw):
ivp = 1. / np.diag(cov.values if hasattr(cov, "values") else cov); return ivp / ivp.sum()
def getClusterVar(cov, cItems):
cov_ = cov.loc[cItems, cItems]; w_ = getIVP(cov_).reshape(-1, 1); return float((w_.T @ cov_.values @ w_)[0, 0])
def getQuasiDiag(link):
link = link.astype(int); sortIx = pd.Series([link[-1, 0], link[-1, 1]]); numItems = link[-1, 3]
while sortIx.max() >= numItems:
sortIx.index = range(0, sortIx.shape[0]*2, 2); df0 = sortIx[sortIx >= numItems]
i = df0.index; j = df0.values - numItems; sortIx[i] = link[j, 0]
df0 = pd.Series(link[j, 1], index=i+1); sortIx = pd.concat([sortIx, df0]).sort_index(); sortIx.index = range(sortIx.shape[0])
return sortIx.tolist()
def getRecBipart(cov, sortIx):
w = pd.Series(1.0, index=sortIx); cItems = [sortIx]
while len(cItems) > 0:
cItems = [i[j:k] for i in cItems for j, k in ((0, len(i)//2), (len(i)//2, len(i))) if len(i) > 1]
for i in range(0, len(cItems), 2):
c0, c1 = cItems[i], cItems[i+1]; a = getClusterVar(cov, c0); b = getClusterVar(cov, c1)
alpha = 1 - a/(a+b); w[c0] *= alpha; w[c1] *= 1 - alpha
return w
def correlDist(corr): return ((1 - corr) / 2.) ** .5
def getHRP(cov, corr):
corr, cov = pd.DataFrame(corr), pd.DataFrame(cov)
sortIx = corr.index[getQuasiDiag(sch.linkage(correlDist(corr).values, "single"))].tolist()
return getRecBipart(cov, sortIx).sort_index()
def getMV(cov, **kw):
V = np.asarray(cov, float); V = V / np.mean(np.diag(V)); n = V.shape[0]
res = minimize(lambda w: w @ V @ w, np.repeat(1/n, n), method="SLSQP", bounds=[(0, 1)]*n,
constraints=({"type": "eq", "fun": lambda w: w.sum()-1},), options={"ftol": 1e-12, "maxiter": 500})
return res.x
print("HRP engine ready")HRP engine ready
The paper has no real data, everything is simulated. Four builders, and they all share one trick: start with a few independent "source" series, then bolt on noisy copies of them. The copies are what create the correlation (and the redundancy you saw collapse in the eigenvalues).
| Function | What it builds | Feeds |
|---|---|---|
generateData_ex(nObs, size0, size1, sigma1) |
size0 independent series + size1 noisy copies of them (seed 12345, so it's exactly his data) |
paper's A.3 example (reference only, not plotted) |
generateData_mc(...) |
same idea, tiny daily-scale returns, plus 2 random shocks (one common to two assets, one asset-specific) | the Monte Carlo + Exhibit 8 |
formBlockMatrix(nBlocks, bSize, bCorr) |
a clean block-diagonal correlation matrix: nBlocks groups of bSize assets, each group correlated at bCorr, zero across groups |
Exhibit 9 |
corr_with_m_correlated(N, m, sigma) |
N assets where the last m are noisy copies of the first N−m "sources" (so there are only N−m real factors) |
Exhibit 1 |
def generateData_ex(nObs, size0, size1, sigma1):
# size0 independent series, then size1 noisy copies of random ones (the paper's A.3)
np.random.seed(12345); random.seed(12345)
x = np.random.normal(0, 1, size=(nObs, size0))
cols = [random.randint(0, size0-1) for _ in range(size1)] # which source each copy tracks
y = x[:, cols] + np.random.normal(0, sigma1, size=(nObs, len(cols))) # copy = source + noise
return pd.DataFrame(np.append(x, y, axis=1), columns=range(1, size0+size1+1)), cols
def generateData_mc(nObs, sLength, size0, size1, mu0, sigma0, sigma1F):
# same, at daily scale, plus a common shock and a specific shock (the paper's A.4)
x = np.random.normal(mu0, sigma0, size=(nObs, size0))
cols = [random.randint(0, size0-1) for _ in range(size1)]
y = x[:, cols] + np.random.normal(0, sigma0*sigma1F, size=(nObs, len(cols)))
x = np.append(x, y, axis=1)
p = np.random.randint(sLength, nObs-1, size=2); x[np.ix_(p, [cols[0], size0])] = np.array([[-.5, -.5], [2, 2]]) # common shock
p = np.random.randint(sLength, nObs-1, size=2); x[p, cols[-1]] = np.array([-.5, 2]) # specific shock
return x, cols
def formBlockMatrix(nBlocks, bSize, bCorr):
# perfect block-diagonal correlation: nBlocks groups of bSize assets, each correlated at bCorr
block = np.ones((bSize, bSize))*bCorr; block[range(bSize), range(bSize)] = 1.
return block_diag(*([block]*nBlocks))
def corr_with_m_correlated(N, m, sigma=0.3, T=2500, seed=1):
# N assets; the last m are noisy copies of the first (N-m) sources -> only N-m real factors
rng = np.random.default_rng(seed); x = rng.normal(size=(T, N))
if m > 0:
src = rng.integers(0, N-m, size=m); x[:, N-m:] = x[:, src] + rng.normal(0, sigma, size=(T, m))
return np.corrcoef(x, rowvar=False)
print("data generators ready")data generators ready
Before the mechanics, one picture frames everything. Classic Markowitz plots every portfolio in risk vs return space; the upper-left boundary is the efficient frontier. The methods in this notebook forecast no returns, so they cannot choose a point along that frontier for its return, they all aim at its leftmost tip: the minimum-variance portfolio (lowest risk). Below: a cloud of simulated long-only portfolios, the frontier, and that min-variance vertex (red star), the corner HRP, IVP and CLA all target.
# a tiny investable universe with expected returns + covariance (2-factor structure)
rng = np.random.default_rng(7); Nef, T = 8, 1500
loads = rng.normal(0, 1, size=(Nef, 2)); fac = rng.normal(0, 0.012, size=(T, 2))
Rd = 0.0004 + fac @ loads.T + rng.normal(0, 0.010, size=(T, Nef)) # daily returns
mu = Rd.mean(0) * 252; Sann = np.cov(Rd, rowvar=0) * 252 # annualised mean & covariance
Wr = rng.dirichlet(np.ones(Nef), size=6000) # random long-only portfolios
p_ret = Wr @ mu; p_vol = np.sqrt(np.einsum("ij,jk,ik->i", Wr, Sann, Wr))
wmv = getMV(Sann); mv_vol = np.sqrt(wmv @ Sann @ wmv); mv_ret = wmv @ mu # minimum-variance portfolio
def ef_point(t): # min variance at a target return (long-only)
r = minimize(lambda w: w @ Sann @ w, np.repeat(1/Nef, Nef), method="SLSQP", bounds=[(0, 1)]*Nef,
constraints=({"type": "eq", "fun": lambda w: w.sum() - 1},
{"type": "eq", "fun": lambda w, t=t: w @ mu - t}), options={"ftol": 1e-12, "maxiter": 300})
return np.sqrt(r.x @ Sann @ r.x) if r.success else np.nan
tg = np.linspace(mv_ret, mu.max(), 30); ef_vol = np.array([ef_point(t) for t in tg])
fig, ax = plt.subplots(figsize=(8, 5))
sc = ax.scatter(p_vol*100, p_ret*100, c=p_ret/np.maximum(p_vol, 1e-9), cmap="viridis", s=8, alpha=.5)
ax.plot(ef_vol*100, tg*100, "k-", lw=2, label="efficient frontier")
ax.scatter([mv_vol*100], [mv_ret*100], color="#B2182B", s=180, marker="*", zorder=5, edgecolor="white",
label=f"minimum-variance ({mv_vol*100:.1f}% vol) <- HRP / IVP / CLA aim here")
ax.set_xlabel("volatility (%, annualised)"); ax.set_ylabel("expected return (%, annualised)")
ax.set_title("Efficient frontier of simulated portfolios; risk-based methods target the min-variance vertex", fontsize=10, weight="bold")
fig.colorbar(sc, label="return / vol (Sharpe-like)"); ax.legend(loc="lower right", fontsize=9); plt.tight_layout(); plt.show()
print(f"minimum-variance portfolio: vol {mv_vol*100:.1f}%, return {mv_ret*100:.1f}% (the leftmost point of the frontier)")
print("None of these methods used the return vector mu; they only minimise risk.")
minimum-variance portfolio: vol 5.8%, return 9.5% (the leftmost point of the frontier) None of these methods used the return vector mu; they only minimise risk.
How we get it: build correlation matrices with corr_with_m_correlated(60, m) for growing m, take eigenvalues with np.linalg.eigvalsh, sort them descending, plot on a log scale. The legend's condition number is eig[0] / eig[-1].
As you add correlated assets, the smallest eigenvalue collapses toward 0 → condition number explodes → inverting the matrix (Markowitz) amplifies noise into wild weights. The number of eigenvalues above the collapse = the real factors (here, 60 − m sources).
N = 60
fig, ax = plt.subplots(figsize=(8, 5))
for m in [0, 20, 40, 55]:
ev = np.sort(np.linalg.eigvalsh(corr_with_m_correlated(N, m)))[::-1]
ax.plot(range(1, N+1), ev, lw=1.8, label=f"{m} correlated ({N-m} real factors, cond = {ev[0]/max(ev[-1],1e-9):.0f})")
ax.set_yscale("log"); ax.set_xlabel("eigenvalue rank (sorted big → small)"); ax.set_ylabel("eigenvalue (log-scale)")
ax.set_title("Exhibit 1, Markowitz's curse", weight="bold"); ax.legend(); plt.tight_layout(); plt.show()
C = corr_with_m_correlated(60,55) ev = np.sort(np.linalg.eigvalsh(C))[::-1] print(ev)
[1.48351376e+01 1.21344198e+01 1.12217269e+01 1.01912266e+01 7.46328104e+00 1.06861684e-01 1.05315080e-01 1.02520862e-01 1.00536888e-01 9.90895298e-02 9.83289199e-02 9.78254464e-02 9.62379670e-02 9.50237862e-02 9.43878219e-02 9.37506214e-02 9.31766637e-02 9.21792876e-02 9.11313693e-02 8.95498728e-02 8.83274131e-02 8.71519519e-02 8.61167592e-02 8.59379589e-02 8.53952176e-02 8.50560581e-02 8.48273243e-02 8.31447757e-02 8.19465089e-02 8.18307034e-02 8.10942360e-02 8.07242110e-02 8.01610240e-02 7.93944343e-02 7.84942774e-02 7.79322204e-02 7.74166213e-02 7.69178670e-02 7.67913272e-02 7.57048414e-02 7.46436011e-02 7.32733036e-02 7.24878749e-02 7.09137399e-02 7.06284666e-02 6.97456464e-02 6.90144000e-02 6.81596365e-02 6.78117051e-02 6.68225047e-02 6.57500112e-02 6.52174134e-02 6.47852127e-02 6.45507336e-02 6.23265280e-02 1.10528986e-02 7.53220342e-03 7.17300763e-03 6.58990192e-03 5.44769382e-03]
Exhibit 1 shows the smallest eigenvalue collapsing as redundant assets pile up. Here is why that is fatal for Markowitz: minimum-variance weights need the inverse covariance,
$$w = \frac{\Sigma^{-1}\mathbf{1}}{\mathbf{1}^{\top}\Sigma^{-1}\mathbf{1}},$$
and in the eigenbasis $\Sigma^{-1}=\sum_i \frac{1}{\lambda_i}\,v_i v_i^{\top}$, so the smallest eigenvalue of $\Sigma$ becomes the largest term of $\Sigma^{-1}$. Writing $\lambda_{\min}=\min_i\lambda_i$ and $\lambda_{\max}=\max_i\lambda_i$ for the smallest and largest eigenvalues of $\Sigma$ (both from np.linalg.eigvalsh), as $\lambda_{\min}\to 0$:
Below we rebuild the same matrices, invert them, and watch the inverse-norm and the leverage explode (left), while a tiny 0.01% perturbation to the correlations moves the weights ever more (right), the weight shift tracks $\kappa$ almost perfectly (log-log corr ≈ 0.97). HRP sidesteps all of this: it never inverts $\Sigma$.
def mv_weights(C): # minimum-variance weights need the inverse covariance
Ci = np.linalg.inv(C); one = np.ones(len(C))
return Ci @ one / (one @ Ci @ one)
N = 60; ms = np.arange(0, 57, 3); rng = np.random.default_rng(0)
kappa, inv_norm, leverage, wsens = [], [], [], []
for m in ms:
C = corr_with_m_correlated(N, m, seed=1); lam = np.linalg.eigvalsh(C)
kappa.append(lam[-1] / max(lam[0], 1e-12)) # condition number k = lam_max / lam_min
inv_norm.append(1.0 / max(lam[0], 1e-12)) # spectral norm of the inverse = 1/lam_min
w = mv_weights(C); leverage.append(np.abs(w).sum()) # gross leverage sum|w_i| (1.0 = fully long-only)
E = rng.normal(0, 1e-4, size=(N, N)); E = (E + E.T) / 2; np.fill_diagonal(E, 0) # 0.01% symmetric noise
wsens.append(np.linalg.norm(mv_weights(C + E) - w) / np.linalg.norm(w)) # relative weight move
kappa = np.array(kappa); wsens = np.array(wsens)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 4.6))
a1.plot(ms, inv_norm, "o-", color="#B2182B", label=r"$\|\Sigma^{-1}\|_2 = 1/\lambda_{\min}$")
a1.plot(ms, leverage, "s-", color="#2A4E7E", label=r"gross leverage $\sum_i|w_i|$")
a1.set_yscale("log"); a1.set_xlabel("# redundant assets m")
a1.set_title("Inverting the covariance explodes as the smallest eigenvalue -> 0", fontsize=10, weight="bold"); a1.legend()
a2.loglog(kappa, wsens, "o", color="#1B7837", label="min-var weights")
kf = np.array([kappa.min(), kappa.max()]); c = np.median(wsens / kappa)
a2.loglog(kf, c * kf, "k:", label=r"slope 1 ($\propto\kappa$)")
a2.set_xlabel(r"condition number $\kappa=\lambda_{\max}/\lambda_{\min}$")
a2.set_ylabel(r"relative weight shift $\|\Delta w\|/\|w\|$")
a2.set_title("A 0.01% input error -> weight error grows with the condition number", fontsize=10, weight="bold"); a2.legend()
plt.tight_layout(); plt.show()
print(f"m=0 : cond={kappa[0]:.0f} ||inv||={inv_norm[0]:.1f} gross leverage={leverage[0]:.2f}")
print(f"m={ms[-1]} : cond={kappa[-1]:.0f} ||inv||={inv_norm[-1]:.0f} gross leverage={leverage[-1]:.1f} (long-only would be 1.00)")
print(f"log-log corr(condition number, weight shift) = {np.corrcoef(np.log(kappa), np.log(wsens))[0,1]:.3f}")
m=0 : cond=2 ||inv||=1.4 gross leverage=1.00 m=54 : cond=1617 ||inv||=144 gross leverage=7.0 (long-only would be 1.00) log-log corr(condition number, weight shift) = 0.972
The same numbers, level by level, every column is built from the two extreme eigenvalues of $\Sigma$:
| redundant assets $m$ | $\lambda_{\min}=\min\operatorname{eig}(\Sigma)$ | $\lambda_{\max}=\max\operatorname{eig}(\Sigma)$ | $\lVert\Sigma^{-1}\rVert=1/\lambda_{\min}$ (red) | $\kappa=\lambda_{\max}/\lambda_{\min}$ (cond) | leverage $\sum_i\lvert w_i\rvert$ |
|---|---|---|---|---|---|
| 0 | 0.726 | 1.33 | 1.4 | 2 | 1.00 |
| 3 | 0.039 | 2.09 | 25.4 | 53 | 1.00 |
| 6 | 0.037 | 2.16 | 27.0 | 58 | 1.01 |
| 15 | 0.020 | 3.86 | 49.3 | 190 | 1.08 |
| 30 | 0.020 | 3.92 | 50.1 | 196 | 1.33 |
| 45 | 0.013 | 6.61 | 78.8 | 521 | 2.20 |
| 54 | 0.007 | 11.24 | 143.9 | 1,617 | 6.98 |
Notice $1/\lambda_{\min}$ (red) leaps 1.4 → 25 with just 3 redundant assets, one near-duplicate already pushes $\lambda_{\min}$ near its noise floor, then plateaus, while $\kappa$ keeps climbing to 1,617 because $\lambda_{\max}$ keeps growing as copies stack onto the same few factors.
The eigenvalue blow-up has a very visible symptom: where the weight actually goes. Min-variance inverts Σ, so it pours capital into the few assets that look lowest-variance and, because it cannot tell near-duplicates apart, zeroes the redundant ones. HRP never inverts Σ; it only splits weight between clusters, so it stays spread out.
Two views on a 30-asset correlated universe (half the assets are noisy copies of the others): - left , the actual weights sorted: CLA/Min-Variance spikes into a handful and zeroes the rest; HRP and IVP stay broad; - right , as we add redundant assets, CLA's top-5 concentration marches toward 100% and it zeroes ever more names, while HRP barely moves.
def make_universe(N, m, T=2500, seed=1):
# N assets with heterogeneous vols; the last m are noisy copies of earlier ones
rng = np.random.default_rng(seed); vol = rng.uniform(0.10, 0.40, N); x = rng.normal(size=(T, N))
if m > 0:
src = rng.integers(0, N - m, size=m); x[:, N - m:] = x[:, src] + rng.normal(0, 0.3, size=(T, m))
r = x * vol
return np.cov(r, rowvar=0), np.corrcoef(r, rowvar=0)
def weights(cov, corr):
return {"HRP": np.asarray(getHRP(cov, corr)), "IVP": getIVP(cov), "CLA/MinVar": getMV(cov)}
top5 = lambda w: np.sort(w)[::-1][:5].sum() * 100 # % held by the 5 biggest positions
nzero = lambda w: int((w < 1e-4).sum()) # how many assets get ~0 weight
col = {"HRP": "#B2182B", "IVP": "#2A4E7E", "CLA/MinVar": "#1B7837"}
N = 30
covD, corrD = make_universe(N, 15, seed=1); W = weights(covD, corrD) # 30 assets, 15 redundant
fig, (a1, a2) = plt.subplots(1, 2, figsize=(13, 4.8), gridspec_kw={"width_ratios": [1.1, 1]})
# (left) sorted weights -- CLA spikes and flatlines at zero; HRP/IVP stay broad
for k in ["IVP", "HRP", "CLA/MinVar"]:
a1.plot(np.sort(W[k])[::-1] * 100, "o-", ms=4, color=col[k],
label=f"{k} (top-5 {top5(W[k]):.0f}%, {nzero(W[k])} zeroed)")
a1.set_xlabel("asset rank (largest weight first)"); a1.set_ylabel("weight (%)")
a1.set_title("Sorted weights on a 30-asset correlated universe", fontsize=10, weight="bold"); a1.legend(fontsize=8)
# (right) concentration vs how many assets are redundant
ms = list(range(0, N - 1, 2)); c5 = {k: [] for k in W}; zc = []
for mm in ms:
cov_, corr_ = make_universe(N, mm, seed=1); Wm = weights(cov_, corr_)
for k in Wm: c5[k].append(top5(Wm[k]))
zc.append(nzero(Wm["CLA/MinVar"]))
for k in ["IVP", "HRP", "CLA/MinVar"]:
a2.plot(ms, c5[k], "o-", ms=3, color=col[k], label=k)
a2.set_xlabel("# redundant assets m"); a2.set_ylabel("top-5 concentration (%)")
a2.set_title("More redundancy -> Markowitz concentrates more", fontsize=10, weight="bold"); a2.legend(fontsize=8)
a2.annotate(f"CLA zeroes {zc[0]} -> {zc[-1]} of {N} assets", xy=(0.5, 0.07), xycoords="axes fraction",
fontsize=8, color="#1B7837", ha="center", weight="bold")
plt.tight_layout(); plt.show()
print("30-asset universe, 15 redundant:")
print(f" top-5 concentration: CLA {top5(W['CLA/MinVar']):.0f}% HRP {top5(W['HRP']):.0f}% IVP {top5(W['IVP']):.0f}%")
print(f" assets zeroed: CLA {nzero(W['CLA/MinVar'])}/{N} HRP {nzero(W['HRP'])} IVP {nzero(W['IVP'])}")
30-asset universe, 15 redundant: top-5 concentration: CLA 59% HRP 44% IVP 40% assets zeroed: CLA 14/30 HRP 0 IVP 0
Now the mechanism, on a small correlated universe so every plot is legible: generateData_ex(10000, 5, 5, .25) builds 10 assets (5 independent sources + 5 noisy copies), and we give them realistic, heterogeneous volatilities (roughly 7% to 30% annualised, from calm "bond-like" up to volatile "tech-like"). Rescaling by volatility leaves the correlations unchanged, so the clustering is identical, but now the risk differences actually drive the allocation. .cov() / .corr() give the matrices the algorithm uses.
nObs, size0, size1, sigma1 = 10000, 5, 5, .25
x, cols = generateData_ex(nObs, size0, size1, sigma1)
# give the assets realistic, heterogeneous volatilities (rescaling leaves correlations -> clustering unchanged)
rng = np.random.default_rng(7)
src_vol = np.array([0.18, 0.08, 0.28, 0.14, 0.11]) # the 5 sources, annualised
copy_vol = np.array([src_vol[cols[i]] * rng.uniform(0.75, 1.15) for i in range(size1)]) # copies ~ their source's vol
asset_vol = np.concatenate([src_vol, copy_vol])
x = x * asset_vol # rescale each asset by its vol
cov, corr = x.cov(), x.corr()
print("asset vols (% annual):", {a: round(v*100) for a, v in zip(x.columns, asset_vol)})
print("correlated pairs (source -> copy):", [(j+1, size0+i) for i, j in enumerate(cols, 1)])asset vols (% annual): {1: 18, 2: 8, 3: 28, 4: 14, 5: 11, 6: 14, 7: 20, 8: 30, 9: 24, 10: 7}
correlated pairs (source -> copy): [(4, 6), (1, 7), (3, 8), (3, 9), (2, 10)]
The calculation, step by step:
- correlDist(corr) maps the 10x10 correlation to a distance $d=\sqrt{\tfrac12(1-\rho)}$: correlated assets are close, uncorrelated are far.
- sch.linkage(..., 'single') builds the tree, 9 merges joining the 10 assets bottom-up (left, the dendrogram). Low merges = tight pairs (a source and its copy); the tallest merges join unrelated clusters.
- getQuasiDiag reads that tree and returns sortIx, the reordering that places similar assets next to each other; corr.loc[sortIx, sortIx] is the same 10x10 matrix with the ~5 blocks now sitting on the diagonal (right), the structure that was invisible in the raw order.
print(f"x: {x.shape[0]:,} daily returns x {x.shape[1]} assets -> correlation matrix {corr.shape[0]}x{corr.shape[1]}")
link = sch.linkage(correlDist(corr).values, "single")
print(f"single-linkage tree: {link.shape[0]} merges joining {corr.shape[0]} assets")
sortIx = corr.index[getQuasiDiag(link)].tolist()
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
sch.dendrogram(link, labels=corr.columns.tolist(), ax=axes[0]); axes[0].set_title("Stage 1, dendrogram (single linkage)", weight="bold")
im = axes[1].imshow(corr.loc[sortIx, sortIx].values, cmap="RdYlBu_r", vmin=-1, vmax=1)
axes[1].set_xticks(range(10)); axes[1].set_xticklabels(sortIx); axes[1].set_yticks(range(10)); axes[1].set_yticklabels(sortIx)
axes[1].set_title("Stage 2, quasi-diagonalized correlation", weight="bold"); fig.colorbar(im, ax=axes[1], fraction=.046)
plt.tight_layout(); plt.show()x: 10,000 daily returns x 10 assets -> correlation matrix 10x10 single-linkage tree: 9 merges joining 10 assets
The diagonal blocks come from one fact: the distance $d=\sqrt{\tfrac12(1-\rho)}$ makes same-cluster assets (a source and its copy, correlation near 1) close, and different-cluster assets (correlation near 0) far. Single linkage fuses each cluster internally before it reaches across to another one. Below, all 45 asset-pairs as correlation vs distance, coloured by same vs different cluster.
member = {s + 1: s for s in range(size0)} # the 5 sources
member.update({size0 + 1 + i: cols[i] for i in range(size1)}) # the 5 noisy copies
D = correlDist(corr)
pairs = [(a, b, corr.loc[a, b], D.loc[a, b], member[a] == member[b])
for a in corr.index for b in corr.index if a < b]
same = [p for p in pairs if p[4]]; diff = [p for p in pairs if not p[4]]
es = max(same, key=lambda p: p[2]); ed = min(diff, key=lambda p: abs(p[2]))
print(f"same cluster (asset {es[0]} & {es[1]}): corr {es[2]:+.2f} -> distance {es[3]:.3f} (close)")
print(f"different cluster (asset {ed[0]} & {ed[1]}): corr {ed[2]:+.2f} -> distance {ed[3]:.3f} (far)")
fig, ax = plt.subplots(figsize=(7.5, 4.6))
ax.scatter([p[2] for p in diff], [p[3] for p in diff], s=70, color="#9AA0A6", edgecolor="w", label=f"different clusters ({len(diff)} pairs)")
ax.scatter([p[2] for p in same], [p[3] for p in same], s=70, color="#B2182B", edgecolor="w", label=f"same cluster ({len(same)} pairs)")
rr = np.linspace(-1, 1, 100); ax.plot(rr, np.sqrt(0.5 * (1 - rr)), color="k", lw=1, ls=":", label=r"$d=\sqrt{(1-\rho)/2}$")
ax.set_xlabel(r"correlation $\rho$"); ax.set_ylabel("distance d")
ax.set_title("Same cluster = high correlation, short distance; different = near 0, long distance", fontsize=10, weight="bold")
ax.legend(); plt.tight_layout(); plt.show()same cluster (asset 2 & 10): corr +0.97 -> distance 0.120 (close) different cluster (asset 3 & 6): corr +0.00 -> distance 0.707 (far)
Now the allocation, and the answer to by cluster or by asset: by cluster first. Cut the dendrogram into its 5 natural clusters (each a source and its copies, sizes 1 to 3), then share capital in two levels:
- across the 5 clusters by risk parity, cluster weight $\propto 1/\mathrm{Var}_{\text{cluster}}$ (the inverse-variance-weighted getClusterVar);
- within each cluster, inverse-variance across its assets (getIVP).
Each asset's final weight = (its cluster's share) × (its share inside the cluster). Now that the clusters carry different risk, the shares diverge sharply: the low-vol cluster commands roughly half the capital, while the high-vol 3-asset cluster gets only a few percent so its members land near 1% each. That is HRP's core move: tilt toward the calmer clusters, and never let a volatile block dominate. (The paper's getRecBipart reaches the same idea by recursively bisecting the sorted order.)
# cut the dendrogram into its 5 natural clusters, then allocate across and within them
cl = sch.fcluster(link, t=5, criterion="maxclust")
clusters = {}
for lab, cid in zip(corr.index, cl): clusters.setdefault(cid, []).append(lab)
clusters = list(clusters.values())
clusters.sort(key=lambda c: min(sortIx.index(a) for a in c)) # order clusters left->right
for c in clusters: c.sort(key=lambda a: sortIx.index(a)) # order assets within a cluster
cvar = [getClusterVar(cov, c) for c in clusters] # each cluster's variance
inv = np.array([1.0/v for v in cvar]); cw = inv / inv.sum() # cluster weights: risk parity across clusters
asset_w = {}
for c, w in zip(clusters, cw):
ivp = getIVP(cov.loc[c, c]) # inverse-variance within the cluster
for a, iv in zip(c, ivp): asset_w[a] = w * iv
print(f"{len(clusters)} clusters from the dendrogram (source + its copies):")
for k, (c, w, v) in enumerate(zip(clusters, cw, cvar), 1):
print(f" C{k} = {c} (variance {v:.3f}) -> {w*100:.0f}% of capital, within: " + ", ".join(f"{a} {asset_w[a]*100:.1f}%" for a in c))
# 2-level tree: 100% -> 5 clusters -> assets
fig, ax = plt.subplots(figsize=(12, 6)); xpos = {a: i for i, a in enumerate(sortIx)}
cx = [np.mean([xpos[a] for a in c]) for c in clusters]; rootx = np.mean(list(xpos.values()))
colors = plt.cm.Set2(np.linspace(0, 1, len(clusters)))
ax.scatter([rootx], [2], s=130, color="#2A4E7E", zorder=3, edgecolor="white")
ax.text(rootx, 2.3, "100% of capital", ha="center", fontsize=10, weight="bold", color="#2A4E7E")
for k, (c, w, xc, colr) in enumerate(zip(clusters, cw, cx, colors), 1):
ax.plot([rootx, xc], [2, 1], "-", color="#B0B0B0", lw=1.4, zorder=1)
ax.text((rootx+xc)/2, 1.52, f"{w*100:.0f}%", fontsize=9, ha="center", va="center",
bbox=dict(boxstyle="round,pad=0.14", fc="white", ec="#B0B0B0", lw=.5), zorder=2)
ax.scatter([xc], [1], s=900, color=colr, zorder=3, edgecolor="#444", lw=1)
ax.text(xc, 1, f"C{k}\n{w*100:.0f}%", ha="center", va="center", fontsize=8.5, weight="bold", color="#222", zorder=4)
for a in c:
xa = xpos[a]; ax.plot([xc, xa], [1, 0], "-", color=colr, lw=1.4, alpha=.7, zorder=1)
ax.scatter([xa], [0], s=470, color="#B2182B", zorder=3, edgecolor="white", lw=1)
ax.text(xa, 0, f"{a}", ha="center", va="center", fontsize=10, weight="bold", color="white", zorder=4)
ax.text(xa, -0.34, f"{asset_w[a]*100:.1f}%", ha="center", va="top", fontsize=9, weight="bold", color="#B2182B", zorder=4)
ax.text(0.5, -0.13, "top branch = cluster's share of capital (risk parity across the 5 clusters); bottom = each asset's final weight",
transform=ax.transAxes, ha="center", fontsize=8.5, color="#555", style="italic")
ax.axis("off"); ax.set_xlim(-0.8, len(sortIx)-0.2); ax.set_ylim(-0.85, 2.65); plt.tight_layout(); plt.show()
print(f"\nall 10 assets keep a positive weight (sum {sum(asset_w.values()):.2f}); weights range from ~{min(asset_w.values())*100:.0f}% to ~{max(asset_w.values())*100:.0f}%.")5 clusters from the dendrogram (source + its copies): C1 = [8, 3, 9] (variance 0.072) -> 4% of capital, within: 8 1.1%, 3 1.3%, 9 1.7% C2 = [2, 10] (variance 0.006) -> 50% of capital, within: 2 22.5%, 10 27.8% C3 = [1, 7] (variance 0.036) -> 8% of capital, within: 1 4.5%, 7 3.5% C4 = [5] (variance 0.012) -> 24% of capital, within: 5 23.5% C5 = [4, 6] (variance 0.020) -> 14% of capital, within: 4 7.3%, 6 6.9%
all 10 assets keep a positive weight (sum 1.00); weights range from ~1% to ~28%.
How we get it: hrpMC runs 2,000 simulated markets. For each one: generateData_mc makes shocked data, then at every rebalance it estimates cov/corr on a 260-day window and pulls weights from getIVP / getHRP / getMV, applies them to the next 22 unseen days, and records the terminal return. We report the variance of those 2,000 terminal returns per method.
Even though min-variance is CLA's objective, HRP has the lowest out-of-sample variance, because it never inverts the covariance. (2,000 iters here; paper used 10,000.)
def hrpMC(numIters=2000, nObs=520, size0=5, size1=5, mu0=0, sigma0=1e-2, sigma1F=.25, sLength=260, rebal=22):
methods = {"IVP": lambda c, r: getIVP(c), "HRP": lambda c, r: getHRP(c, r).values, "CLA/MinVar": lambda c, r: getMV(c)}
terms = {k: [] for k in methods}; paths = {k: [] for k in methods}; pointers = range(sLength, nObs, rebal)
for _ in range(int(numIters)):
xx, _ = generateData_mc(nObs, sLength, size0, size1, mu0, sigma0, sigma1F)
r = {k: [] for k in methods}
for p in pointers:
xin = xx[p-sLength:p]; cov_ = np.cov(xin, rowvar=0); corr_ = np.corrcoef(xin, rowvar=0); xout = xx[p:p+rebal]
for k, f in methods.items(): r[k].append(xout @ f(cov_, corr_))
for k in methods:
eq = np.concatenate([[1.0], np.cumprod(1 + np.concatenate(r[k]))]) # keep the whole equity path, starts at $1
terms[k].append(eq[-1] - 1); paths[k].append(eq)
return pd.DataFrame(terms), {k: np.array(v) for k, v in paths.items()}
np.random.seed(0); random.seed(0)
S, PATHS = hrpMC(numIters=2000); var = S.var()
paper = pd.Series({"IVP": 0.0928, "HRP": 0.0671, "CLA/MinVar": 0.1157}) # paper's Exhibit (CLA, 10,000 sims)
pd.DataFrame({"std": S.std(), "variance": var, "paper variance (10,000 sims)": paper,
"extra variance vs HRP": var/var["HRP"] - 1}).round(4)| std | variance | paper variance (10,000 sims) | extra variance vs HRP | |
|---|---|---|---|---|
| IVP | 0.3137 | 0.0984 | 0.0928 | 0.3702 |
| HRP | 0.2680 | 0.0718 | 0.0671 | 0.0000 |
| CLA/MinVar | 0.3496 | 0.1222 | 0.1157 | 0.7021 |
How: the same 2,000 Monte-Carlo runs, but now we keep each run's equity path (not just its final value) and draw 80 of them per method. Same result as the table, but you can see it: HRP's paths bundle tightly (low variance); CLA's spray apart (high variance), even though minimizing variance is CLA's whole objective. This is the payoff picture.
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5), sharey=True)
col = {"HRP": "#B2182B", "IVP": "#2A4E7E", "CLA/MinVar": "#1B7837"}
for ax, m in zip(axes, ["HRP", "IVP", "CLA/MinVar"]):
for p in PATHS[m][:80]: ax.plot(p, color=col[m], lw=.5, alpha=.22)
ax.axhline(1, color="k", lw=.6, ls=":"); ax.set_xlabel("out-of-sample day")
ax.set_title(f"{m} (variance = {S[m].var():.3f})", weight="bold")
axes[0].set_ylabel("growth of $1")
fig.suptitle("Fan of out-of-sample paths, HRP bundles tight, CLA sprays wide", weight="bold", y=1.03)
plt.tight_layout(); plt.show()
The variance the paper reports is itself a sample statistic: too few simulations and it's noisy, enough and it settles. Drag the slider (100 → 2,000) and watch the HRP picture converge. Because 1 simulation = 1 equity path, the fan literally fills in as you add runs (drawn up to 800 to keep the file sane; the variance / σ readout uses all n). Notice the outline barely moves past a few hundred runs, the terminal-return distribution is fixed, so more sims just estimate the same spread more precisely. All 2,000 runs are computed above (PATHS["HRP"]); the slider only reveals more of them.
(Same Plotly slider style as the Greeks interactive, it stays live in the saved notebook, no kernel required.)
import numpy as np, plotly.graph_objects as go, plotly.io as pio
pio.renderers.default = "notebook"
# Reuse the 2,000 HRP runs already computed in the Monte Carlo above, nothing is recomputed.
POOL_PATHS = PATHS["HRP"] # (2000, n_days+1): each row an equity curve, starts at $1
POOL_TERMS = S["HRP"].to_numpy() # each run's terminal return
DAYS = np.arange(POOL_PATHS.shape[1])
OPTIONS = [100, 200, 400, 600, 800, 1000, 1200, 1400, 1700, 2000] # the 10 slider stops
# ---- gamma-interactive palette + gold slider ----
BG, INK, INK_MUTE, GRID, GOLD = "#000000", "#FFFFFF", "#B8CAD9", "#3A3A3A", "#FFC93A"
HRP_NEON = "#FF3B57"
CAP = 800 # max individual paths drawn (keeps the notebook a sane size)
STEP = 3 # plot every 3rd day; float32 to keep the embedded arrays small
# 1 simulation = 1 path, so the fan literally densifies as you add sims (drawn count = min(n, CAP)).
# All n paths are packed into ONE NaN-separated trace per frame, cheap, unlike n separate traces.
def fan_trace(n):
d = DAYS[::STEP]; k = min(n, CAP)
xs = np.empty((k, d.size + 1), np.float32); ys = np.empty((k, d.size + 1), np.float32)
xs[:, :-1] = d; xs[:, -1] = np.nan # NaN breaks the line between paths
ys[:, :-1] = POOL_PATHS[:k, ::STEP]; ys[:, -1] = np.nan
return go.Scatter(x=xs.ravel(), y=ys.ravel(), mode="lines", hoverinfo="skip",
line=dict(color="rgba(255,59,87,0.10)", width=1), showlegend=False, name="paths")
def frame_traces(n):
mean = go.Scatter(x=DAYS, y=POOL_PATHS[:n].mean(0).astype(np.float32), mode="lines",
line=dict(color=HRP_NEON, width=3.5), name="mean path")
return [fan_trace(n), mean] # constant 2 traces per frame -> clean animation
def var_note(n):
v, sd = POOL_TERMS[:n].var(), POOL_TERMS[:n].std()
return [dict(x=0.02, y=0.98, xref="paper", yref="paper", xanchor="left", yanchor="top",
showarrow=False, align="left",
text=(f"<b>{n:,} simulations</b><br>variance = {v:.4f}<br>std = {sd:.4f}"),
font=dict(color=INK, size=14, family="Consolas, monospace"),
bgcolor="rgba(255,59,87,0.12)", bordercolor=HRP_NEON, borderwidth=1, borderpad=8)]
DEF = len(OPTIONS) - 1 # start on 2,000
frames = [go.Frame(data=frame_traces(n), name=str(n),
layout=go.Layout(annotations=var_note(n))) for n in OPTIONS]
fig = go.Figure(data=frame_traces(OPTIONS[DEF]), frames=frames)
fig.add_hline(y=1.0, line=dict(color=INK_MUTE, width=1, dash="dot"))
fig.update_layout(
paper_bgcolor=BG, plot_bgcolor=BG,
title=dict(text=("<b style='font-size:24px;color:white'>HRP, Monte Carlo out-of-sample</b><br>"
"<span style='font-size:13px;color:#B8CAD9;font-family:Consolas,monospace'>"
"drag the slider — watch the fan & variance settle as simulations pile up</span>"),
x=0.5, xanchor="center", y=0.97),
xaxis=dict(title="out-of-sample day", color=INK_MUTE, gridcolor=GRID, zerolinecolor=GRID),
yaxis=dict(title="growth of $1", color=INK_MUTE, gridcolor=GRID, zerolinecolor=GRID),
annotations=var_note(OPTIONS[DEF]),
legend=dict(x=0.98, xanchor="right", y=0.06, font=dict(color=INK_MUTE)),
sliders=[dict(
active=DEF, pad={"t": 50, "b": 10},
currentvalue={"prefix": "simulations = ",
"font": {"size": 16, "color": GOLD, "family": "Consolas, monospace"}},
len=0.92, x=0.04, y=-0.06, ticklen=4, tickcolor=GOLD,
bgcolor="rgba(255,201,58,0.18)", activebgcolor=GOLD,
font={"color": INK_MUTE, "size": 11}, bordercolor="rgba(255,255,255,0.05)",
steps=[dict(method="animate",
args=[[str(n)], {"mode": "immediate",
"frame": {"duration": 0, "redraw": True},
"transition": {"duration": 0}}],
label=f"{n:,}") for n in OPTIONS],
)],
margin=dict(l=70, r=30, t=110, b=90), height=560,
font=dict(family="IBM Plex Sans, Segoe UI, sans-serif"),
)
fig.show()How: one_run() replays a single simulated market and stores the weight array from getIVP / getHRP / getMV at every rebalance; plt.stackplot draws each method's weights over time. Each panel's title also reports its turnover (avg sum|Δw| per rebalance): CLA's bands lurch from period to period while HRP's drift smoothly, that is HRP trading less.
def one_run(nObs=520, size0=5, size1=5, mu0=0, sigma0=1e-2, sigma1F=.25, sLength=260, rebal=22, seed=3):
np.random.seed(seed); random.seed(seed)
xx, _ = generateData_mc(nObs, sLength, size0, size1, mu0, sigma0, sigma1F)
ptr = list(range(sLength, nObs, rebal))
methods = {"IVP": lambda c, r: getIVP(c), "HRP": lambda c, r: getHRP(c, r).values, "CLA/MinVar": lambda c, r: getMV(c)}
W = {m: [] for m in methods}
for p in ptr:
xin = xx[p-sLength:p]; cov_ = np.cov(xin, rowvar=0); corr_ = np.corrcoef(xin, rowvar=0)
for m, f in methods.items(): W[m].append(f(cov_, corr_))
return {m: np.array(W[m]) for m in methods}, ptr
W, ptr = one_run()
turn1 = {m: np.abs(np.diff(W[m], axis=0)).sum(1).mean() for m in W} # avg L1 weight move per rebalance
fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)
for ax, m in zip(axes, ["IVP", "HRP", "CLA/MinVar"]):
ax.stackplot(range(len(ptr)), W[m].T, labels=[str(i+1) for i in range(10)], colors=plt.cm.tab10(np.linspace(0,1,10)))
ax.set_ylabel("weight"); ax.set_ylim(0, 1)
ax.set_title(f"Exhibit 8, allocations over time: {m} (avg turnover {turn1[m]:.3f})", weight="bold")
axes[-1].set_xlabel("rebalance #"); axes[0].legend(ncol=10, fontsize=7, loc="upper center")
plt.tight_layout(); plt.show()
One run is noisy, so here is the average over 60 simulated markets: at each rebalance we measure how much each method's weights move (sum|Δw|). HRP shifts far less than min-variance, it never inverts Σ, so estimation noise does not whip its weights around, which means lower trading costs for the same, or lower, out-of-sample risk.
def avg_turnover(nsim=60, nObs=520, size0=5, size1=5, mu0=0, sigma0=1e-2, sigma1F=.25, sLength=260, rebal=22):
methods = {"IVP": lambda c, r: getIVP(c), "HRP": lambda c, r: getHRP(c, r).values, "CLA/MinVar": lambda c, r: getMV(c)}
t = {k: [] for k in methods}; np.random.seed(1); random.seed(1)
for _ in range(nsim):
xx, _ = generateData_mc(nObs, sLength, size0, size1, mu0, sigma0, sigma1F); prev = {k: None for k in methods}
for p in range(sLength, nObs, rebal):
xin = xx[p-sLength:p]; cov_ = np.cov(xin, rowvar=0); corr_ = np.corrcoef(xin, rowvar=0)
for k, f in methods.items():
w = f(cov_, corr_)
if prev[k] is not None: t[k].append(np.abs(w - prev[k]).sum())
prev[k] = w
return {k: float(np.mean(v)) for k, v in t.items()}
tv = avg_turnover()
order = ["IVP", "HRP", "CLA/MinVar"]; col = {"IVP": "#2A4E7E", "HRP": "#B2182B", "CLA/MinVar": "#1B7837"}
fig, ax = plt.subplots(figsize=(6.5, 4.2))
bars = ax.bar(order, [tv[k] for k in order], color=[col[k] for k in order], width=.6)
for b, k in zip(bars, order): ax.text(b.get_x()+b.get_width()/2, b.get_height()+.003, f"{tv[k]:.3f}", ha="center", fontsize=11, weight="bold")
ax.set_ylabel("avg turnover per rebalance (sum of |weight change|)")
ax.set_title("HRP rebalances less than min-variance (lower turnover)", fontsize=10, weight="bold")
plt.tight_layout(); plt.show()
print("avg turnover: IVP {:.3f} HRP {:.3f} CLA/MinVar {:.3f}".format(tv["IVP"], tv["HRP"], tv["CLA/MinVar"]))
print("HRP trades {:.0f}% less than CLA/MinVar".format((1 - tv["HRP"]/tv["CLA/MinVar"]) * 100))
avg turnover: IVP 0.073 HRP 0.122 CLA/MinVar 0.186 HRP trades 34% less than CLA/MinVar
How: formBlockMatrix(10, 20, 0.6) builds a 200-asset block matrix; we shuffle the order (left, no visible structure), then getQuasiDiag on its single-linkage tree reorders it back into clean blocks (right).
big = pd.DataFrame(formBlockMatrix(10, 20, 0.6))
np.random.seed(5); cols_b = big.columns.tolist(); np.random.shuffle(cols_b)
bigS = big[cols_b].loc[cols_b]
order = bigS.index[getQuasiDiag(sch.linkage(correlDist(bigS).values, "single"))].tolist()
after = bigS.loc[order, order]
fig, (a1, a2) = plt.subplots(1, 2, figsize=(12, 5.5))
a1.imshow(bigS.values, cmap="viridis"); a1.set_title("Exhibit 9a, before clustering (shuffled)", weight="bold")
im = a2.imshow(after.values, cmap="viridis"); a2.set_title("Exhibit 9b, after clustering", weight="bold")
fig.colorbar(im, ax=a2, fraction=.046); plt.tight_layout(); plt.show()
The core exhibits, on his own simulated data, each traced to the function that made it:
- Ex 1 the problem, corr_with_m_correlated + np.linalg.eigvalsh show the condition number exploding as correlated assets are added; inverting Σ then blows up ($\|\Sigma^{-1}\|=1/\lambda_{\min}$), so leverage and input-error sensitivity scale with κ (Markowitz's curse). HRP never inverts Σ.
- Why Markowitz concentrates, on a 30-asset correlated universe getMV zeroes the redundant names and piles into a few (top-5 concentration climbs toward 100% as redundancy grows), while getHRP and getIVP stay broad.
- Ex 8 to 9, one_run shows HRP's weights stay smooth; getQuasiDiag turns a shuffled 200-asset matrix into 10 clean blocks.
- Monte Carlo, hrpMC (2,000 runs): HRP has the lowest out-of-sample variance, because it never inverts the covariance and weights only move within clusters.