Risk parity, minimum variance, HRP and HERC.

Portfolios that only look at risk

Lesson 4.1.1 left an uncomfortable conclusion: mean-variance optimization is exact with parameters that are never known, and the worst estimated parameter is the expected return. In the synthetic universe of this level, the 20-year sample mean of an asset has a standard error equal to its volatility divided by √20: for a stock market with 20% volatility, about 4.5 points a year. With an error like that, the optimizer mostly picks noise.

There is a family of portfolios that gives up on estimating μ and only uses the covariance matrix Σ, which is estimated much better (lesson 4.1.2). Each one answers a different question:

They all have a mean-variance interpretation: each one is the optimal portfolio under an assumption about μ. Minimum variance is optimal if all assets have the same expected return; maximum diversification, if the expected return is proportional to volatility (all with the same Sharpe); risk parity, if in addition all correlations are equal (Maillard, Roncalli and Teïletche, 2010). Choosing one of them is choosing an assumption.

The universe is that of universo4_precios.parquet: 38 fictitious total return indices in USD, from 2006 to 2025, generated with a published factor model. Its true covariance, by regime and overall, is in universo4_cov_verdad.csv. It is synthetic data: it serves to measure how far each method is from the truth, not to decide what works in real markets. The money market and the 1-3 year Treasury bill are removed, since they are almost cash (further down you will see what happens if the bill is left in).

import numpy as np
import pandas as pd

precios = pd.read_parquet("universo4_precios.parquet")
precios = precios.set_index("fecha")
fact = pd.read_parquet("universo4_factores_diarios.parquet")
fact = fact.set_index("fecha")
activos = pd.read_csv("universo4_activos.csv", index_col="activo")

r = precios.pct_change()
r.iloc[0] = precios.iloc[0] / 100 - 1   # base 100 the day before
exceso = r.sub(fact["rf"], axis=0)

cv = pd.read_csv("universo4_cov_verdad.csv")
def cov_verdad(regimen, cols):
    m = cv[cv["regimen"] == regimen].set_index("activo")
    return m.loc[cols, cols]

cols = [c for c in precios if c not in ("MONET", "RF_T13")]
S = cov_verdad("total", cols)          # annual, true
clase = activos.loc[cols, "clase"]
print(len(cols), "assets")
clase.value_counts()

36 assets
clase
renta_variable     18
renta_fija          9
materias_primas     5
inmobiliario        2
alternativo         2
Name: count, dtype: int64

Returns are simple and the return of a portfolio is the sum of weights times returns. With logarithmic returns the sum is not exact and the Sharpe of a portfolio comes out lower than it is. Half of the universe is equities: stock markets by region, with and without currency hedging, sectors and styles. That redundancy, many assets that are almost the same risk, is what puts the methods of this lesson to the test.

Risk contributions

The volatility of a portfolio, $\sigma_p = \sqrt{w'\Sigma w}$, is a homogeneous function of degree 1 of the weights: if all of them double, it doubles. By Euler's theorem, it decomposes exactly into the sum of what each asset contributes:

$RC_i$ is the fraction of the risk contributed by asset $i$. It can be negative: an asset negatively correlated with the rest of the portfolio reduces its volatility. The classic 60/40 is the usual example:

def contribuciones(w, S):
    """Fraction of the portfolio's risk contributed by each asset."""
    w = pd.Series(w, index=S.index)
    sw = S @ w                         # Σw
    return w * sw / (w @ sw)           # they sum to 1

def vol(w, S):
    w = np.asarray(w)
    return float(np.sqrt(w @ np.asarray(S) @ w))

w6040 = pd.Series(0.0, index=cols)
w6040[["RV_EEUU", "RF_T710"]] = [0.6, 0.4]
print(f"true volatility: {vol(w6040, S):.2%}")
contribuciones(w6040, S)[["RV_EEUU", "RF_T710"]].round(3)

true volatility: 11.87%
activo
RV_EEUU    1.01
RF_T710   -0.01
dtype: float64

60% of the capital is in stocks and so is 101% of the risk. The 7-10 year bond, negatively correlated with stocks in the dominant regime of the universe, contributes −1%: it reduces volatility somewhat, but carries no risk of its own. The figure at the top summarizes it for five portfolios: the hollow square is the capital weight of equities (with real estate and alternatives, which behave like them) and the filled one, their share of the risk.

The 1/N portfolio looks diversified because it splits the capital, but not the risk:

def por_clase(w, S):
    w = pd.Series(w, index=S.index)
    return pd.DataFrame({
        "capital": w.groupby(clase).sum(),
        "riesgo": contribuciones(w, S).groupby(clase).sum(),
    })

w_1n = pd.Series(1 / len(cols), index=cols)
print(f"true volatility: {vol(w_1n, S):.2%}")
por_clase(w_1n, S).round(3)

true volatility: 13.41%
                 capital  riesgo
clase
alternativo        0.056   0.089
inmobiliario       0.056   0.084
materias_primas    0.139   0.087
renta_fija         0.250   0.023
renta_variable     0.500   0.718

Fixed income has 25% of the capital and 2.3% of the risk. Equities, 50% of the capital and 71.8% of the risk; with real estate and alternatives, 89%. Splitting the money equally among assets with very different volatilities means concentrating the risk in the volatile ones.

Risk parity

The equal risk contribution portfolio (ERC) asks for $RC_i = 1/N$ for all. It has no closed form except in special cases, but it does have a convex formulation (Maillard, Roncalli and Teïletche, 2010): the minimum of

satisfies $y_i(\Sigma y)_i = b_i$; normalizing $w = y / \sum y$, the contributions are proportional to the budgets $b_i$. With $b_i = 1/N$ it is the ERC. The problem is solved well by coordinate descent (Griveau-Billion, Richard and Roncalli, 2013): with the others fixed, the condition for asset $i$ is a quadratic equation in $y_i$ with a single positive root.

def erc(S, presupuesto=None, vueltas=200):
    """Risk parity by coordinate descent: minimizes
    ½·y'Σy − Σ b_i·ln(y_i) and normalizes the weights to sum to 1."""
    S = np.asarray(S)
    n = len(S)
    b = np.full(n, 1 / n) if presupuesto is None else presupuesto
    y = 1 / np.sqrt(np.diag(S))
    for _ in range(vueltas):
        for i in range(n):
            a = S[i] @ y - S[i, i] * y[i]
            disc = a * a + 4 * S[i, i] * b[i]
            y[i] = (-a + np.sqrt(disc)) / (2 * S[i, i])
    return pd.Series(y / y.sum(), index=cols)

w_erc = erc(S)
rc = contribuciones(w_erc, S)
print(f"contributions: from {rc.min():.4f} to {rc.max():.4f}")
print(f"true volatility: {vol(w_erc, S):.2%}")
por_clase(w_erc, S).round(3)

contributions: from 0.0278 to 0.0278
true volatility: 6.01%
                 capital  riesgo
clase
alternativo        0.019   0.056
inmobiliario       0.019   0.056
materias_primas    0.113   0.139
renta_fija         0.650   0.250
renta_variable     0.198   0.500

The 36 contributions are 0.0278, that is, 1/36. The split by class is exactly the number of assets in each class divided by 36: 18 equity assets give 50% of the risk, 9 fixed income 25%, 5 commodities 13.9%. It is the first practical lesson of risk parity: it equalizes per asset, and the split by class is decided by whoever designs the universe. Ten more sector funds would raise the equity share without the market having changed. If what you want is a third of the risk in stocks, another in bonds and another in commodities, you have to say so with the budgets $b_i$ (exercise 1).

The second: the portfolio has 65% of its capital in bonds and a volatility of 6.01%. To reach 10% it would have to be leveraged 1.66 times, and that is what risk parity funds do with futures. Asness, Frazzini and Pedersen (2012) explain its historical return precisely that way: many investors cannot or will not use leverage and overpay for volatile assets, so a leveraged portfolio of low-volatility assets has had a better Sharpe. Leverage brings its own risks (funding, margin, correlations that rise together, lesson 4.2.4), and the volatility targeting of lesson 4.2.2 is the other half of the method.

Before going on, it is worth checking that the custom code matches the libraries. skfolio and Riskfolio-Lib estimate the covariance from returns, so the comparison is made with a sample, that of 2024-2025:

import riskfolio as rp
from skfolio import RiskMeasure
from skfolio.optimization import RiskBudgeting

X = exceso.loc["2024":"2025", cols]   # two years of sample
S_m = X.cov()
propia = erc(S_m)

skf = RiskBudgeting(risk_measure=RiskMeasure.VARIANCE).fit(X)
port = rp.Portfolio(returns=X)
port.assets_stats(method_mu="hist", method_cov="hist")
rkf = port.rp_optimization(model="Classic", rm="MV", b=None)

print(f"skfolio:  {np.abs(skf.weights_ - propia).max():.1e}")
print(f"riskfolio: {np.abs(rkf['weights'] - propia).max():.1e}")

skfolio:  8.1e-06
riskfolio: 2.4e-06

The differences are of the order of $10^{-6}$: both libraries solve the same formulation with a conic optimizer and stop at the solver's tolerance. They match.

Minimum variance and maximum diversification

Long-only minimum variance is a quadratic problem: minimize $w'\Sigma w$ with $\sum w = 1$ and $w \ge 0$ (the unconstrained closed form is in 4.1.1). The maximum diversification of Choueifaty and Coignard (2008) maximizes the diversification ratio

the weighted average of the volatilities divided by the volatility of the portfolio. It equals 1 if everything is perfectly correlated and grows the more the assets offset each other. Maximizing it is equivalent to finding the minimum variance on the correlation matrix and then dividing each weight by its volatility; that is how it is solved below, because the correlation matrix is better conditioned than the covariance matrix and the solver works more comfortably.

import cvxpy as cp

def min_var(S):
    n = len(S)
    w = cp.Variable(n)
    L = np.linalg.cholesky(np.asarray(S))
    riesgo = cp.sum_squares(L.T @ w)        # = w'Σw
    cp.Problem(cp.Minimize(riesgo),
               [cp.sum(w) == 1, w >= 0]).solve(solver="CLARABEL")
    w = np.clip(w.value, 0, None)
    return pd.Series(w / w.sum(), index=S.index)

def max_div(S):
    """Maximum diversification = minimum variance on the correlation
    matrix, then each weight is divided by its volatility."""
    sd = np.sqrt(np.diag(S))
    y = min_var(S / np.outer(sd, sd))
    w = y / sd
    return w / w.sum()

def ratio_div(w, S):
    sd = np.sqrt(np.diag(S))
    return float(w @ sd / vol(w, S))

def n_efectivo(w):
    return float(1 / (np.asarray(w) ** 2).sum())

w_mv, w_md = min_var(S), max_div(S)
for nombre, w in [("1/N", w_1n), ("ERC", w_erc),
                  ("min. variance", w_mv), ("max. divers.", w_md)]:
    print(f"{nombre:14s} vol {vol(w, S):6.2%}  "
          f"DR {ratio_div(w, S):4.2f}  N_eff {n_efectivo(w):5.1f}")
w_mv[w_mv > 0.01].sort_values(ascending=False).round(3)

1/N            vol 13.41%  DR 1.32  N_eff  36.0
ERC            vol  6.01%  DR 1.77  N_eff  14.6
min. variance  vol  3.11%  DR 1.68  N_eff   1.3
max. divers.   vol  4.87%  DR 2.09  N_eff   6.5
activo
RF_T37    0.879
RF_CHY    0.046
SC_FIN    0.026
MP_AGR    0.025
ES_VAL    0.010
dtype: float64

cp.sum_squares(L.T @ w) is $w'\Sigma w$ written with the Cholesky factor $\Sigma = LL'$: it is the form cvxpy recognizes as convex without checking that the matrix is semidefinite, and the one that accepts parameters if the problem is solved many times. The effective number of assets, $N_{eff} = 1/\sum w_i^2$, measures concentration: 36 in the 1/N, 14.6 in the ERC, 6.5 in maximum diversification and 1.3 in minimum variance, which puts 87.9% in the 3-7 year Treasury bond. It is the least volatile portfolio (3.11%) and the least diversified. Maximum diversification achieves the highest DR (2.09) with six effective assets.

The concentration of minimum variance is a feature of the method, not a mistake. With the 1-3 year bill in the universe:

cols13 = cols + ["RF_T13"]
w13 = min_var(cov_verdad("total", cols13))
w13[w13 > 0.01].round(3)

activo
SC_FIN    0.020
RF_T13    0.975
dtype: float64

97.5% in the bill: the minimum variance portfolio of a universe with an almost riskless asset is that asset. That is why it was removed at the start, and why risk-based methods are applied to universes of risky assets, often leveraged afterwards.

HRP step by step

The hierarchical risk parity (HRP) of López de Prado (2016) starts from a criticism of quadratic optimizers: inverting Σ amplifies its errors, the more so the worse conditioned it is, and a matrix of 36 assets with nearly duplicated sectors is very badly conditioned. HRP avoids the inversion with four steps:

flowchart TB
  C[Covariance Σ] --> D[Distance d_ij]
  D --> A[Hierarchical tree]
  A --> O[Leaf order]
  O --> B[Bisection]
  B --> W[HRP weights]

1. Distance. From the correlation to a metric distance: $d_{ij} = \sqrt{(1 - \rho_{ij})/2}$, which is 0 with correlation 1 and 1 with correlation −1.
2. Tree. A hierarchical clustering (in the paper, single linkage: the distance between two clusters is that of their closest members) joins the most similar assets first and then the clusters.
3. Quasi-diagonalization. The assets are sorted according to the leaves of the tree. With that order, similar assets end up together and the reordered covariance has the large values near the diagonal.
4. Recursive bisection. The sorted list is split in half. Each half gets a variance: that of its inverse variance portfolio. The weight is split between the two halves in inverse proportion to those variances, $\alpha = 1 - V_{L}/(V_{L} + V_{R})$, and the process is repeated within each half until single assets are reached.

from scipy.cluster.hierarchy import leaves_list, linkage
from scipy.spatial.distance import squareform

def var_ivp(S, ix):
    """Variance of the inverse variance portfolio of group ix."""
    c = S[np.ix_(ix, ix)]
    w = 1 / np.diag(c)
    w /= w.sum()
    return w @ c @ w

def hrp(S, metodo="single"):
    S = np.asarray(S)
    sd = np.sqrt(np.diag(S))
    corr = S / np.outer(sd, sd)
    # 1. correlation distance
    d = np.sqrt(np.clip((1 - corr) / 2, 0, 1))
    # 2. hierarchical tree (condensed distance vector)
    Z = linkage(squareform(d, checks=False), metodo)
    # 3. quasi-diagonalization: the leaf order
    orden = leaves_list(Z)
    # 4. recursive bisection
    w = np.ones(len(S))
    grupos = [orden]
    while grupos:
        siguientes = []
        for g in grupos:
            if len(g) < 2:
                continue
            izq, der = g[: len(g) // 2], g[len(g) // 2 :]
            v_izq, v_der = var_ivp(S, izq), var_ivp(S, der)
            alfa = 1 - v_izq / (v_izq + v_der)
            w[izq] *= alfa
            w[der] *= 1 - alfa
            siguientes += [izq, der]
        grupos = siguientes
    return w

w_hrp = pd.Series(hrp(S), index=cols)
print(f"true volatility: {vol(w_hrp, S):.2%}, "
      f"N_eff {n_efectivo(w_hrp):.1f}")
por_clase(w_hrp, S).round(3)

true volatility: 4.18%, N_eff 8.2
                 capital  riesgo
clase
alternativo        0.006   0.013
inmobiliario       0.005   0.012
materias_primas    0.062   0.078
renta_fija         0.863   0.793
renta_variable     0.064   0.105

With the true covariance, HRP lies between risk parity and minimum variance: a volatility of 4.18%, 8.2 effective assets and 86.3% of the capital in fixed income, which contributes 79.3% of the risk. It is not an equal contribution portfolio: the bisection allocates to each half in inverse proportion to its variance, so the bond half, much less volatile, takes almost all the weight at the first cut.

The tree and the cut can be seen in the next figure. The leaf order is this:

sd = np.sqrt(np.diag(S))
d = np.sqrt(np.clip((1 - S / np.outer(sd, sd)) / 2, 0, 1))
Z = linkage(squareform(d.values, checks=False), "single")
print(" ".join(np.array(cols)[leaves_list(Z)]))

MP_ORO RF_EUR MP_AGR MP_IND MP_AMP MP_ENE RF_INF RF_GLB_C RF_T20 RF_T37 RF_T710 RF_CIG RF_EMG RF_CHY SC_UTI RV_JAP RV_JAP_C INFRA INM_GLB INM_EEUU SC_ENE SC_CBA RV_EMG RV_PAC SC_TEC SC_SAL RV_EUR RV_EUR_C SC_FIN ES_VAL RV_EEUU_P ALT_PRI ES_MOM SC_IND RV_EEUU ES_MIN

On the left, gold, unhedged European debt and commodities, which join late; then Treasury, credit and emerging market bonds; and on the right, a chain of stock markets, sectors, real estate and alternatives that join one by one. That chain is typical of single linkage (chaining): a single close neighbor is enough to join a cluster. The dashed line marks the first bisection: 18 assets on each side, and the cut falls between INFRA and INM_GLB, inside the equity chain, not at a branch of the tree. Bisection by halves uses the order, not the structure of the tree. That is the criticism that motivated HERC, below.

The same HRP in skfolio and Riskfolio-Lib

Both libraries implement HRP with more options (risk measures other than variance, other linkages, optimal leaf ordering). With the paper's configuration, single linkage and no leaf reordering, they should give the same as the custom code:

from skfolio.cluster import HierarchicalClustering, LinkageMethod
from skfolio.optimization import HierarchicalRiskParity
from skfolio.seriation import HierarchicalSeriation

propia = hrp(S_m)

hc = HierarchicalClustering(linkage_method=LinkageMethod.SINGLE)
ser = HierarchicalSeriation(hierarchical_clustering_estimator=hc,
                            optimal_ordering=False)
skf = HierarchicalRiskParity(seriation_estimator=ser).fit(X)

port = rp.HCPortfolio(returns=X)
rkf = port.optimization(model="HRP", codependence="pearson",
                        rm="MV", linkage="single",
                        leaf_order=False)

print(f"skfolio:   {np.abs(skf.weights_ - propia).max():.1e}")
print(f"riskfolio: {np.abs(rkf['weights'] - propia).max():.1e}")

# with each library's default values
skf_def = HierarchicalRiskParity().fit(X).weights_
rkf_def = rp.HCPortfolio(returns=X).optimization(model="HRP")
print(f"default, skfolio:   "
      f"{np.abs(skf_def - propia).max():.3f}")
print(f"default, riskfolio: "
      f"{np.abs(rkf_def['weights'] - propia).max():.3f}")

skfolio:   8.3e-17
riskfolio: 1.7e-16
default, skfolio:   0.123
default, riskfolio: 0.111

They match up to machine rounding ($10^{-16}$). With the default values, on the other hand, the weight of some asset changes by 12.3 points (skfolio) and 11.1 points (Riskfolio-Lib). skfolio uses Ward linkage and optimal leaf ordering by default; Riskfolio-Lib, single linkage with optimal ordering (leaf_order=True). The leaf order decides where each bisection falls, so two “identical” HRPs can give different portfolios. Before comparing one result with another, you have to read which configuration each one used.

The advantage nobody disputes is another one: HRP does not need the covariance to be invertible. With fewer sessions than assets, the sample covariance is singular:

X20 = exceso.loc["2025-12", cols].tail(20)
S20 = X20.cov().values * 252
print("rank of the covariance:", np.linalg.matrix_rank(S20))
w20 = hrp(S20)
print(f"HRP: weights from {w20.min():.4f} to {w20.max():.4f}")
np.linalg.cholesky(S20)

rank of the covariance: 19
HRP: weights from 0.0008 to 0.1968
numpy.linalg.LinAlgError: Matrix is not positive definite

With 20 sessions the matrix has rank 19 (one degree of freedom is lost by subtracting the mean). The Cholesky factorization that min_var used fails and the inverse does not exist; HRP gives reasonable weights, because it only uses the diagonal and quadratic forms of blocks. That it works does not mean that 20 sessions are enough to estimate anything; it means that the method does not break.

HERC: equalizing risk across clusters

Raffinot (2017, 2018) proposed two changes to HRP. The first: cut the tree at its branches and not at the middle of the list. The second: choose a number of clusters $k$ (with a statistical criterion, such as the gap statistic) and split the risk equally across clusters, top-down through the tree, and within each cluster naively (inverse volatility or equal weights). This is the hierarchical equal risk contribution (HERC).

from skfolio.optimization import HierarchicalEqualRiskContribution

herc = HierarchicalEqualRiskContribution().fit(X)
grupos = herc.hierarchical_clustering_estimator_.n_clusters_
w_herc = pd.Series(herc.weights_, index=cols)
print("clusters:", grupos)

rkf_herc = rp.HCPortfolio(returns=X).optimization(
    model="HERC", rm="MV", linkage="ward")["weights"]
pd.DataFrame({
    "herc_skfolio": w_herc.groupby(clase).sum(),
    "herc_riskfolio": rkf_herc.groupby(clase).sum(),
    "hrp": pd.Series(propia, cols).groupby(clase).sum(),
}).round(3)

clusters: 7
                 herc_skfolio  herc_riskfolio    hrp
clase
alternativo             0.069           0.051  0.007
inmobiliario            0.056           0.042  0.008
materias_primas         0.003           0.014  0.054
renta_fija              0.159           0.364  0.855
renta_variable          0.713           0.529  0.076

With the 2024-2025 sample, skfolio finds 7 clusters and puts 71.3% of the capital in equities; Riskfolio-Lib, with the same Ward linkage and its own criterion for the number of clusters, 52.9%; HRP, 7.6%. HERC splits the risk across clusters: if the tree separates the stock markets into several clusters and leaves the bonds in one, equities get several slices of risk and bonds only one. The result depends on how many clusters there are and on how each library forms them, and that changes from one sample to another. It is its conceptual strength (it uses the structure of the tree) and its practical weakness (the portfolio jumps when $k$ changes).

Out of sample with the truth known

López de Prado (2016) compared HRP with the minimum variance of the critical line algorithm (CLA) and with inverse variance in a Monte Carlo: 10 normal series with random jumps and a random correlation structure, a window of 260 observations and rebalancing every 22. The out-of-sample variance was 0.1157 for CLA, 0.0928 for inverse variance and 0.0671 for HRP: minimum variance, whose goal is precisely minimum variance, turned out to be the worst. The question here is whether that happens in this universe.

With the true covariance it can be measured without a backtest. Sixty sessions (or 252) of normal returns are simulated with the covariance of the calm regime, the sample covariance is estimated, the weights of each method are computed and the true volatility of those weights is measured, $\sqrt{w'\Sigma w}$ with the true Σ. This is repeated 60 times.

S_calma = cov_verdad("calma", cols).values
n = len(cols)
optimo = vol(min_var(pd.DataFrame(S_calma, cols, cols)), S_calma)
L = np.linalg.cholesky(S_calma / 252)
rng = np.random.default_rng(414)

def todos(Sh, Xd):
    Sd = pd.DataFrame(Sh, cols, cols)
    iv = 1 / np.sqrt(np.diag(Sh))
    return {
        "1/N": np.full(n, 1 / n),
        "inv_vol": iv / iv.sum(),
        "erc": erc(Sh, vueltas=60).values,
        "min_var": min_var(Sd).values,
        "max_div": max_div(Sd).values,
        "hrp": hrp(Sh),
        "herc": HierarchicalEqualRiskContribution().fit(Xd).weights_,
    }

filas = []
for T in (60, 252):
    for k in range(60):
        Xd = rng.standard_normal((T, n)) @ L.T
        for m, w in todos(np.cov(Xd.T) * 252, Xd).items():
            filas.append((T, m, vol(w, S_calma)))
mc = pd.DataFrame(filas, columns=["T", "metodo", "vol"])
print(f"true optimum: {optimo:.2%}")
tabla = mc.groupby(["metodo", "T"])["vol"].agg(
    ["mean", lambda v: v.quantile(0.05), lambda v: v.quantile(0.95)])
tabla.columns = ["media", "p5", "p95"]
(100 * tabla.unstack("T")).round(2)

true optimum: 2.84%
        media          p5          p95
T         60    252   60    252    60     252
metodo
1/N      9.99  9.99  9.99  9.99   9.99   9.99
erc      5.23  5.29  4.75  5.04   5.70   5.57
herc     6.47  7.94  3.82  4.05  10.26  10.72
hrp      4.02  4.07  3.78  3.86   4.27   4.22
inv_vol  6.64  6.73  6.23  6.49   7.05   6.95
max_div  4.81  4.50  3.48  3.67   7.02   5.48
min_var  2.95  2.87  2.87  2.85   3.08   2.92

The result is not the paper's. Minimum variance estimated with only 60 sessions has a true volatility of 2.95%, 0.11 points from the true optimum (2.84%), and with 252 sessions, 2.87%. HRP ends up at 4.02%, with a narrow range (3.78-4.27%): it is stable, but it does not minimize variance, because it splits the weight by halves and does not look for the minimum. HERC is the most unstable: between 3.8% and 10.3% depending on the sample, because the number of clusters changes.

The difference from the paper lies in the universe. Here the covariance comes from a ten-factor model with independent idiosyncratic terms, the portfolio is long only (the $w \ge 0$ constraint is in itself a form of regularization) and there are no asset-specific jumps to punish a concentrated portfolio. In the paper's experiment the random jumps do punish concentration. Neither result is “the truth” about HRP: each one is true for its data generator. What can be stated in general is that HRP does not break with badly conditioned or singular matrices and that its portfolio is more diversified than the minimum variance one.

The second experiment is a monthly rebalancing from 2008 to 2025: on the first day of each month, each method computes its weights with the previous two years; between rebalances the portfolio is held and the weights drift with prices. Besides the realized metrics, the last column gives the average true Sharpe of the weights chosen: $w'\mu / \sqrt{w' \Sigma w}$ with the true μ and Σ of the universe. It separates what each method chose from the luck it had afterwards.

V = 504                                   # two-year window
fechas = exceso.index
primeras = pd.Series(fechas, fechas).groupby(
    fechas.to_period("M")).first()
rebal = [d for d in primeras if fechas.get_loc(d) >= V]

metodos = ["1/N", "inv_vol", "erc", "min_var", "max_div", "hrp",
           "herc"]
pesos = {m: {} for m in metodos}
for d in rebal:
    i = fechas.get_loc(d)
    Xd = exceso[cols].iloc[i - V : i]
    for m, w in todos(Xd.cov().values * 252, Xd).items():
        pesos[m][d] = w

def backtest(P):
    """Buy and hold between rebalances (the weights drift)."""
    trozos, rot, previo = [], [], None
    for j, d in enumerate(P.index):
        fin = P.index[j + 1] if j + 1 < len(P) else None
        b = r.loc[d:fin, cols]
        b = b.iloc[:-1] if fin is not None else b
        w = P.loc[d].values
        if previo is not None:
            rot.append(np.abs(w - previo).sum())
        g = (1 + b).cumprod().values
        valor = g @ w
        trozos.append(pd.Series(
            valor / np.r_[1, valor[:-1]] - 1, b.index))
        previo = g[-1] * w / valor[-1]
    total = pd.concat(trozos)
    return total - fact["rf"].loc[total.index], np.mean(rot) * 12

mu_v = pd.read_csv("universo4_verdad.csv",
                   index_col="activo").loc[cols, "mu_verdadero"]
res = {}
for m in pesos:
    P = pd.DataFrame(pesos[m], index=cols).T
    ex_m, rot_anual = backtest(P)
    curva = (1 + ex_m + fact["rf"].loc[ex_m.index]).cumprod()
    res[m] = {
        "exceso_%": 252 * ex_m.mean() * 100,
        "vol_%": np.sqrt(252) * ex_m.std() * 100,
        "sharpe": np.sqrt(252) * ex_m.mean() / ex_m.std(),
        "max_caida_%": 100 * (curva / curva.cummax() - 1).min(),
        "rotacion_%": 100 * rot_anual,
        "n_ef": np.mean([n_efectivo(w) for w in P.values]),
        "sharpe_verdad": np.mean([w @ mu_v / vol(w, S)
                                  for w in P.values]),
    }
print(rebal[0].date(), "to", fechas[-1].date())
pd.DataFrame(res).T.round(2)

2008-02-01 to 2025-12-31
         exceso_%  vol_%  sharpe  max_caida_%  rotacion_%   n_ef  sharpe_verdad
1/N          3.93  13.79    0.28       -46.19       37.40  36.00           0.39
inv_vol      2.65   8.51    0.31       -26.91       34.89  22.76           0.47
erc          1.85   6.25    0.30       -18.67       37.42  15.43           0.50
min_var      0.41   3.21    0.13       -10.77       38.52   1.31           0.45
max_div      0.94   6.06    0.16       -17.68      165.83   5.51           0.52
hrp          0.91   4.45    0.20       -18.24      178.78   8.93           0.48
herc         6.87  10.33    0.66       -21.72      399.08  16.74           0.43

Five readings of the table:

• The order of the volatilities is that of the Monte Carlo: 1/N 13.79%, inverse volatility 8.51%, ERC 6.25%, maximum diversification 6.06%, HRP 4.45% and minimum variance 3.21%, with maximum drawdowns ranging from −46% for the 1/N to −11% for minimum variance. The methods do what they say: less risk the more they concentrate in what is not very volatile.
• The realized Sharpes (0.13 to 0.66) do not say which method is better. With 18 years, the standard error of an annual Sharpe is about $1/\sqrt{18} \approx 0.24$. The true Sharpes of what each method chose range from 0.39 (1/N) to 0.52 (maximum diversification), and the realized ones of the methods with more bonds fall below because the universe's bonds returned less in the sample than their expectation (RF_T710: 0.5% average excess return in 2006-2025 against a true 1.6%).
• HERC has the best realized Sharpe (0.66) and one of the worst true ones (0.43). It weighted equities heavily in years when they paid well. It is exactly the kind of result the backtest rewards and the truth does not back up.
• Turnover separates the methods more than the Sharpe. 1/N, inverse volatility, ERC and minimum variance turn over between 35% and 39% a year; maximum diversification and HRP, 166% and 179%; HERC, 399%. Turnover costs money (lesson 4.1.5) and it is the symptom of an unstable method.
• Minimum variance was not the unstable one. Its turnover (39%) is that of the ERC. The commonplace that quadratic optimizers jump and HRP does not is reversed on this data.

Where does HRP's turnover come from? From the bisection by halves. A small change in the correlations can move an asset from one end of the tree's chain to the other; then which assets fall on each side of a cut changes and the weight of a whole cluster jumps:

P_hrp = pd.DataFrame(pesos["hrp"], index=cols).T
P_erc = pd.DataFrame(pesos["erc"], index=cols).T
salto = P_hrp.diff().abs().sum(axis=1)
d = salto.idxmax()
antes = P_hrp.index[P_hrp.index.get_loc(d) - 1]
cambio = (P_hrp.loc[d] - P_hrp.loc[antes])
top = cambio.abs().sort_values(ascending=False).index[:4]
print(d.date(), f"HRP turnover {salto[d]:.0%}, ERC "
      f"{P_erc.diff().abs().sum(axis=1)[d]:.0%}")
pd.DataFrame({"antes": P_hrp.loc[antes, top],
              "después": P_hrp.loc[d, top]}).round(3)

2015-12-01 HRP turnover 69%, ERC 3%
          antes  después
RF_T37    0.211    0.429
RF_GLB_C  0.242    0.068
RF_INF    0.183    0.095
RF_CIG    0.036    0.067

On December 1, 2015, with one more month of data in a two-year window, HRP moved 69% of the portfolio and risk parity 3%. The 3-7 year bond went from 21.1% to 42.9% and hedged global debt, from 24.2% to 6.8%. The figure follows the weight of the 3-7 year bond in the two portfolios over the whole period:

Risk parity moves slowly, because each weight depends continuously on the covariance. HRP depends on a discrete object, the leaf order, and it jumps. The variants that smooth the tree (another linkage, optimal leaf ordering) move the jumps around but do not remove them, because they all still cut a sorted list (exercise 2 measures how much the weights change with the tree).

Which method and when

flowchart TB
  A{Reliable μ?} -- yes --> B[Mean-variance]
  A -- no --> C{Goal?}
  C -- less risk --> D[Minimum variance]
  C -- spread risk --> E{Redundant assets?}
  E -- no --> F[ERC]
  E -- yes --> G[ERC by groups]
  G --> H[or HRP]

• Minimum variance when the goal is literally the lowest volatility and concentrated portfolios are acceptable. With per-asset limits (lesson 4.1.1) it is forced to diversify.
• ERC when you want to spread the risk and you control the universe: few assets, different from each other, or budgets per class. It is stable and easy to explain. It needs leverage to earn like a portfolio with stocks.
• HRP when there are many redundant assets, few observations or a covariance that cannot be inverted. Its turnover has to be watched and its configuration (linkage, leaf order) set explicitly.
• HERC as an idea (risk by clusters of the tree), with care: the number of clusters changes with the sample and the portfolio with it. Fixing $k$ or setting the clusters by hand (by asset class) gives something more stable: it is an ERC by groups.
• None of them removes correlation risk. In a crisis correlations rise together (in the universe's stress regime they all come from a single risk-aversion component), and a portfolio calibrated on two calm years has more risk than it measures. That is lesson 4.2.4.

Typical mistakes with these methods:

• Comparing methods by the Sharpe of a backtest. With 18 years, differences of 0.2 cannot be told apart from noise.
• Forgetting turnover: a method with 400% annual turnover pays costs that a backtest without costs does not see.
• Reading “risk parity” as “risk split across classes”: it is split across assets.
• Using a library's default values without saying so: the same method name gives different weights.

Exercises

1. Risk budgets by class. Using the erc function, build a portfolio in which equities (with real estate and alternatives), fixed income and commodities contribute a third of the risk each, split equally within each group. Check the contributions by class and the volatility.

grupo = clase.replace({"inmobiliario": "renta_variable",
                       "alternativo": "renta_variable"})
tam = grupo.map(grupo.value_counts())
b = (1 / 3 / tam).values                 # they sum to 1
w_b = erc(S, presupuesto=b)
print(f"true volatility: {vol(w_b, S):.2%}")
pd.DataFrame({
    "capital": w_b.groupby(grupo).sum(),
    "riesgo": contribuciones(w_b, S).groupby(grupo).sum(),
}).round(3)

true volatility: 5.48%
                 capital  riesgo
clase
materias_primas    0.196   0.333
renta_fija         0.661   0.333
renta_variable     0.142   0.333

2. HRP variants. With the 2024-2025 sample covariance (S_m), compute the maximum weight difference between the lesson's HRP and: (a) the tree of the paper's code, passing the square distance matrix to linkage; (b) Ward linkage. Which one changes more?

def hrp_arbol(S, Z):
    """Recursive bisection on the order of a given tree."""
    S = np.asarray(S)
    w, grupos = np.ones(len(S)), [leaves_list(Z)]
    while grupos:
        sig = []
        for g in grupos:
            if len(g) < 2:
                continue
            a, b = g[: len(g) // 2], g[len(g) // 2 :]
            va, vb = var_ivp(S, a), var_ivp(S, b)
            w[a] *= 1 - va / (va + vb)
            w[b] *= va / (va + vb)
            sig += [a, b]
        grupos = sig
    return w

sd_m = np.sqrt(np.diag(S_m))
d_m = np.sqrt(np.clip((1 - S_m.values / np.outer(sd_m, sd_m)) / 2,
                      0, 1))
import warnings
from scipy.cluster.hierarchy import ClusterWarning

base = hrp(S_m)
with warnings.catch_warnings():     # scipy warns about the square matrix
    warnings.simplefilter("ignore", ClusterWarning)
    articulo = hrp_arbol(S_m, linkage(d_m, "single"))
ward = hrp(S_m, metodo="ward")
print(f"paper's tree:  {np.abs(articulo - base).max():.3f}")
print(f"Ward linkage:  {np.abs(ward - base).max():.3f}")

paper's tree:  0.118
Ward linkage:  0.112

3. Without code: a manager proposes moving from a 60/40 to a risk parity portfolio “because it diversifies better” and keeping the same volatility. Explain with the lesson's numbers what changes in the capital, what is needed to match the volatility and what new risk appears.