Risk parity, mínima varianza, HRP y HERC.

Carteras que solo miran el riesgo

La lección 4.1.1 dejó una conclusión incómoda: la optimización de media-varianza es exacta con parámetros que nunca se conocen, y el parámetro peor estimado es el rendimiento esperado. En el universo sintético de este nivel, la media muestral de 20 años de un activo tiene un error típico de su volatilidad dividida por √20: para una bolsa con un 20 % de volatilidad, unos 4,5 puntos al año. Con un error así, el optimizador elige sobre todo ruido.

Hay una familia de carteras que renuncia a estimar μ y solo usa la matriz de covarianzas Σ, que se estima mucho mejor (lección 4.1.2). Cada una responde a una pregunta distinta:

Todas tienen una interpretación de media-varianza: cada una es la cartera óptima bajo un supuesto sobre μ. La mínima varianza lo es si todos los activos tienen el mismo rendimiento esperado; la máxima diversificación, si el rendimiento esperado es proporcional a la volatilidad (todos con el mismo Sharpe); la paridad de riesgo, si además todas las correlaciones son iguales (Maillard, Roncalli y Teïletche, 2010). Elegir una de ellas es elegir un supuesto.

El universo es el de universo4_precios.parquet: 38 índices ficticios de rentabilidad total en USD, de 2006 a 2025, generados con un modelo factorial publicado. Su covarianza verdadera, por régimen y total, está en universo4_cov_verdad.csv. Son datos sintéticos: sirven para medir cuánto se aleja cada método de la verdad, no para decidir qué funciona en los mercados reales. Se quitan el monetario y la letra del Tesoro a 1-3 años, que son casi efectivo (más abajo se ve qué pasa si se deja la letra).

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 la víspera
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)          # anual, verdadera
clase = activos.loc[cols, "clase"]
print(len(cols), "activos")
clase.value_counts()

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

Los rendimientos son simples y el de una cartera es la suma de pesos por rendimientos. Con rendimientos logarítmicos la suma no es exacta y el Sharpe de una cartera sale más bajo de lo que es. La mitad del universo es renta variable: bolsas por regiones, con y sin cobertura de divisa, sectores y estilos. Esa redundancia, muchos activos que son casi el mismo riesgo, es la que pone a prueba a los métodos de esta lección.

Contribuciones al riesgo

La volatilidad de una cartera, $\sigma_p = \sqrt{w'\Sigma w}$, es una función homogénea de grado 1 de los pesos: si se duplican todos, se duplica. Por el teorema de Euler, se descompone exactamente en la suma de lo que aporta cada activo:

$RC_i$ es la fracción del riesgo que aporta el activo $i$. Puede ser negativa: un activo con correlación negativa con el resto de la cartera reduce su volatilidad. La 60/40 clásica es el ejemplo de siempre:

def contribuciones(w, S):
    """Fracción del riesgo de la cartera que aporta cada activo."""
    w = pd.Series(w, index=S.index)
    sw = S @ w                         # Σw
    return w * sw / (w @ sw)           # suman 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"volatilidad verdadera: {vol(w6040, S):.2%}")
contribuciones(w6040, S)[["RV_EEUU", "RF_T710"]].round(3)

volatilidad verdadera: 11.87%
activo
RV_EEUU    1.01
RF_T710   -0.01
dtype: float64

El 60 % del capital está en acciones y el 101 % del riesgo también. El bono a 7-10 años, con correlación negativa con las acciones en el régimen dominante del universo, aporta un −1 %: reduce algo la volatilidad, pero no carga con riesgo propio. La figura de la cabecera lo resume para cinco carteras: el cuadrado hueco es el peso en capital de la renta variable (con inmobiliario y alternativos, que se comportan como ella) y el lleno, su parte del riesgo.

La cartera 1/N parece diversificada porque reparte el capital, pero no el riesgo:

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"volatilidad verdadera: {vol(w_1n, S):.2%}")
por_clase(w_1n, S).round(3)

volatilidad verdadera: 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

La renta fija tiene el 25 % del capital y el 2,3 % del riesgo. La renta variable, el 50 % del capital y el 71,8 % del riesgo; con inmobiliario y alternativos, el 89 %. Repartir el dinero por igual entre activos de volatilidades muy distintas es concentrar el riesgo en los volátiles.

Paridad de riesgo

La cartera de contribuciones iguales (ERC, equal risk contribution) pide $RC_i = 1/N$ para todos. No tiene fórmula cerrada salvo casos particulares, pero sí una formulación convexa (Maillard, Roncalli y Teïletche, 2010): el mínimo de

cumple $y_i(\Sigma y)_i = b_i$; normalizando $w = y / \sum y$, las contribuciones son proporcionales a los presupuestos $b_i$. Con $b_i = 1/N$ es la ERC. El problema se resuelve bien por descenso coordenado (Griveau-Billion, Richard y Roncalli, 2013): fijados los demás, la condición del activo $i$ es una ecuación de segundo grado en $y_i$ con una única raíz positiva.

def erc(S, presupuesto=None, vueltas=200):
    """Paridad de riesgo por descenso coordenado: minimiza
    ½·y'Σy − Σ b_i·ln(y_i) y normaliza los pesos para que sumen 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"contribuciones: de {rc.min():.4f} a {rc.max():.4f}")
print(f"volatilidad verdadera: {vol(w_erc, S):.2%}")
por_clase(w_erc, S).round(3)

contribuciones: de 0.0278 a 0.0278
volatilidad verdadera: 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

Las 36 contribuciones valen 0,0278, es decir, 1/36. El reparto por clase es exactamente el número de activos de cada clase entre 36: 18 de renta variable dan el 50 % del riesgo, 9 de renta fija el 25 %, 5 de materias primas el 13,9 %. Es la primera lección práctica de la paridad de riesgo: iguala por activo, y el reparto por clase lo decide quien diseña el universo. Diez fondos sectoriales más subirían la parte de la renta variable sin que el mercado haya cambiado. Si lo que se quiere es un tercio del riesgo en acciones, otro en bonos y otro en materias primas, hay que decirlo con los presupuestos $b_i$ (ejercicio 1).

La segunda: la cartera tiene el 65 % del capital en bonos y una volatilidad del 6,01 %. Para llegar a un 10 % habría que apalancarla 1,66 veces, y es lo que hacen los fondos de paridad de riesgo con futuros. Asness, Frazzini y Pedersen (2012) explican su rendimiento histórico justamente por ahí: muchos inversores no pueden o no quieren apalancarse y pagan de más por los activos volátiles, así que una cartera de activos poco volátiles apalancada ha tenido mejor Sharpe. El apalancamiento trae sus propios riesgos (financiación, margen, correlaciones que suben juntas, lección 4.2.4), y el volatility targeting de la lección 4.2.2 es la otra mitad del método.

Antes de seguir, conviene comprobar que el código propio coincide con las librerías. skfolio y Riskfolio-Lib estiman la covarianza desde los rendimientos, así que la comparación se hace con una muestra, la de 2024-2025:

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

X = exceso.loc["2024":"2025", cols]   # dos años de muestra
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

Las diferencias son del orden de $10^{-6}$: las dos librerías resuelven la misma formulación con un optimizador cónico y se paran en la tolerancia del solver. Coinciden.

Mínima varianza y máxima diversificación

La mínima varianza solo con largos es un problema cuadrático: minimizar $w'\Sigma w$ con $\sum w = 1$ y $w \ge 0$ (la fórmula cerrada sin límites está en 4.1.1). La máxima diversificación de Choueifaty y Coignard (2008) maximiza el ratio de diversificación

la media ponderada de las volatilidades entre la volatilidad de la cartera. Vale 1 si todo está perfectamente correlado y crece cuanto más se compensan los activos. Maximizarlo equivale a buscar la mínima varianza sobre la matriz de correlación y dividir después cada peso por su volatilidad; así se resuelve abajo, porque la matriz de correlación está mejor condicionada que la de covarianza y el solver trabaja más cómodo.

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):
    """Máxima diversificación = mínima varianza sobre la matriz de
    correlación, y después se divide cada peso por su volatilidad."""
    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),
                  ("mín. varianza", w_mv), ("máx. divers.", w_md)]:
    print(f"{nombre:14s} vol {vol(w, S):6.2%}  "
          f"DR {ratio_div(w, S):4.2f}  N_ef {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_ef  36.0
ERC            vol  6.01%  DR 1.77  N_ef  14.6
mín. varianza  vol  3.11%  DR 1.68  N_ef   1.3
máx. divers.   vol  4.87%  DR 2.09  N_ef   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) es $w'\Sigma w$ escrito con el factor de Cholesky $\Sigma = LL'$: es la forma que cvxpy reconoce como convexa sin comprobar que la matriz es semidefinida, y la que admite parámetros si el problema se resuelve muchas veces. El número efectivo de activos, $N_{ef} = 1/\sum w_i^2$, mide la concentración: 36 en la 1/N, 14,6 en la ERC, 6,5 en la máxima diversificación y 1,3 en la mínima varianza, que pone un 87,9 % en el bono del Tesoro a 3-7 años. Es la cartera menos volátil (3,11 %) y la menos diversificada. La máxima diversificación consigue el mayor DR (2,09) con seis activos efectivos.

La concentración de la mínima varianza es un rasgo del método, no un error. Con la letra a 1-3 años en el universo:

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

Un 97,5 % en la letra: la cartera de mínima varianza de un universo con un activo casi sin riesgo es ese activo. Por eso se quitó al principio y por eso los métodos de riesgo se aplican sobre universos de activos con riesgo, a menudo apalancados después.

HRP paso a paso

La paridad de riesgo jerárquica (HRP, hierarchical risk parity) de López de Prado (2016) parte de una crítica a los optimizadores cuadráticos: invertir Σ amplifica sus errores, más cuanto peor condicionada está, y una matriz de 36 activos con sectores casi duplicados lo está mucho. HRP evita la inversión con cuatro pasos:

flowchart TB
  C[Covarianza Σ] --> D[Distancia d_ij]
  D --> A[Árbol jerárquico]
  A --> O[Orden de hojas]
  O --> B[Bisección]
  B --> W[Pesos HRP]

1. Distancia. De la correlación a una distancia métrica: $d_{ij} = \sqrt{(1 - \rho_{ij})/2}$, que vale 0 con correlación 1 y 1 con correlación −1.
2. Árbol. Un agrupamiento jerárquico (en el artículo, enlace simple: la distancia entre dos grupos es la de sus miembros más cercanos) une primero los activos más parecidos y después los grupos.
3. Cuasi-diagonalización. Se ordenan los activos según las hojas del árbol. Con ese orden, los activos parecidos quedan juntos y la covarianza reordenada tiene los valores grandes cerca de la diagonal.
4. Bisección recursiva. Se parte la lista ordenada por la mitad. Cada mitad recibe una varianza: la de su cartera de varianza inversa. El peso se reparte entre las dos mitades en proporción inversa a esas varianzas, $\alpha = 1 - V_{izq}/(V_{izq} + V_{der})$, y se repite dentro de cada mitad hasta llegar a activos sueltos.

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

def var_ivp(S, ix):
    """Varianza de la cartera de varianza inversa del grupo 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. distancia de correlación
    d = np.sqrt(np.clip((1 - corr) / 2, 0, 1))
    # 2. árbol jerárquico (vector condensado de distancias)
    Z = linkage(squareform(d, checks=False), metodo)
    # 3. cuasi-diagonalización: el orden de las hojas
    orden = leaves_list(Z)
    # 4. bisección recursiva
    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"volatilidad verdadera: {vol(w_hrp, S):.2%}, "
      f"N_ef {n_efectivo(w_hrp):.1f}")
por_clase(w_hrp, S).round(3)

volatilidad verdadera: 4.18%, N_ef 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

Con la covarianza verdadera, HRP queda entre la paridad de riesgo y la mínima varianza: volatilidad del 4,18 %, 8,2 activos efectivos y un 86,3 % del capital en renta fija, que aporta el 79,3 % del riesgo. No es una cartera de contribuciones iguales: la bisección reparte a cada mitad en proporción inversa a su varianza, de modo que la mitad de los bonos, mucho menos volátil, se lleva casi todo el peso en el primer corte.

El árbol y el corte se ven en la figura siguiente. El orden de las hojas es este:

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

A la izquierda, el oro, la deuda europea sin cubrir y las materias primas, que se unen tarde; después los bonos del Tesoro, de crédito y emergentes; y a la derecha, una cadena de bolsas, sectores, inmobiliario y alternativos que se van uniendo de uno en uno. Esa cadena es típica del enlace simple (chaining): basta un vecino cercano para unirse a un grupo. La línea discontinua marca la primera bisección: 18 activos a cada lado, y el corte cae entre INFRA e INM_GLB, dentro de la cadena de renta variable, no en una rama del árbol. La bisección por mitades usa el orden, no la estructura del árbol. Es la crítica que motivó HERC, más abajo.

El mismo HRP en skfolio y Riskfolio-Lib

Las dos librerías implementan HRP con más opciones (medidas de riesgo distintas de la varianza, otros enlaces, orden óptimo de las hojas). Con la configuración del artículo, enlace simple y sin reordenar las hojas, deben dar lo mismo que el código propio:

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}")

# con los valores por defecto de cada librería
skf_def = HierarchicalRiskParity().fit(X).weights_
rkf_def = rp.HCPortfolio(returns=X).optimization(model="HRP")
print(f"por defecto, skfolio:   "
      f"{np.abs(skf_def - propia).max():.3f}")
print(f"por defecto, riskfolio: "
      f"{np.abs(rkf_def['weights'] - propia).max():.3f}")

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

Coinciden hasta el redondeo de la máquina ($10^{-16}$). Con los valores por defecto, en cambio, el peso de algún activo cambia en 12,3 puntos (skfolio) y 11,1 puntos (Riskfolio-Lib). skfolio usa por defecto el enlace de Ward y el orden óptimo de las hojas; Riskfolio-Lib, el enlace simple con orden óptimo (leaf_order=True). El orden de las hojas decide dónde cae cada bisección, así que dos HRP «iguales» pueden dar carteras distintas. Antes de comparar un resultado con otro, hay que leer qué configuración usó cada uno.

La ventaja que nadie discute es otra: HRP no necesita que la covarianza sea invertible. Con menos sesiones que activos, la covarianza muestral es singular:

X20 = exceso.loc["2025-12", cols].tail(20)
S20 = X20.cov().values * 252
print("rango de la covarianza:", np.linalg.matrix_rank(S20))
w20 = hrp(S20)
print(f"HRP: pesos de {w20.min():.4f} a {w20.max():.4f}")
np.linalg.cholesky(S20)

rango de la covarianza: 19
HRP: pesos de 0.0008 a 0.1968
numpy.linalg.LinAlgError: Matrix is not positive definite

Con 20 sesiones la matriz tiene rango 19 (se pierde un grado de libertad al restar la media). La factorización de Cholesky que usaba min_var falla y la inversa no existe; HRP da unos pesos razonables, porque solo usa la diagonal y formas cuadráticas de bloques. Que funcione no quiere decir que 20 sesiones basten para estimar nada; quiere decir que el método no se rompe.

HERC: igualar el riesgo entre grupos

Raffinot (2017, 2018) propuso dos cambios a HRP. El primero: cortar el árbol por sus ramas y no por la mitad de la lista. El segundo: elegir un número de grupos $k$ (con un criterio estadístico, como el gap) y repartir el riesgo por igual entre grupos, de arriba abajo por el árbol, y dentro de cada grupo de forma ingenua (volatilidad inversa o a partes iguales). Es la contribución igual al riesgo jerárquica (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("grupos:", 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)

grupos: 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

Con la muestra de 2024-2025, skfolio encuentra 7 grupos y pone un 71,3 % del capital en renta variable; Riskfolio-Lib, con el mismo enlace de Ward y su propio criterio para el número de grupos, un 52,9 %; HRP, un 7,6 %. HERC reparte el riesgo entre grupos: si el árbol separa las bolsas en varios grupos y deja los bonos en uno, la renta variable recibe varias porciones de riesgo y los bonos una sola. El resultado depende de cuántos grupos haya y de cómo los forme cada librería, y eso cambia de una muestra a otra. Es su fortaleza conceptual (usa la estructura del árbol) y su debilidad práctica (la cartera salta cuando cambia $k$).

Fuera de muestra con la verdad conocida

López de Prado (2016) comparó HRP con la mínima varianza del algoritmo de línea crítica (CLA) y con la varianza inversa en un Monte Carlo: 10 series normales con saltos aleatorios y una estructura de correlación aleatoria, una ventana de 260 observaciones y rebalanceo cada 22. La varianza fuera de muestra fue 0,1157 para CLA, 0,0928 para la varianza inversa y 0,0671 para HRP: la mínima varianza, cuyo objetivo es precisamente la varianza mínima, resultó la peor. La pregunta aquí es si eso pasa en este universo.

Con la covarianza verdadera se puede medir sin backtest. Se simulan 60 sesiones (o 252) de rendimientos normales con la covarianza del régimen de calma, se estima la covarianza muestral, se calculan los pesos de cada método y se mide la volatilidad verdadera de esos pesos, $\sqrt{w'\Sigma w}$ con la Σ de verdad. Se repite 60 veces.

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"óptimo verdadero: {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)

óptimo verdadero: 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

El resultado no es el del artículo. La mínima varianza estimada con solo 60 sesiones tiene una volatilidad verdadera del 2,95 %, a 0,11 puntos del óptimo verdadero (2,84 %), y con 252 sesiones del 2,87 %. HRP queda en el 4,02 %, con un rango estrecho (3,78-4,27 %): es estable, pero no minimiza la varianza, porque reparte el peso por mitades y no busca el mínimo. HERC es el más inestable: entre el 3,8 % y el 10,3 % según la muestra, porque el número de grupos cambia.

La diferencia con el artículo está en el universo. Aquí la covarianza sale de un modelo de diez factores con idiosincráticos independientes, la cartera es solo de largos (la restricción $w \ge 0$ es en sí misma una forma de regularización) y no hay saltos propios de un activo que castiguen a una cartera concentrada. En el experimento del artículo los saltos aleatorios sí castigan la concentración. Ninguno de los dos resultados es «la verdad» sobre HRP: cada uno es cierto para su generador de datos. Lo que se puede afirmar en general es que HRP no se rompe con matrices mal condicionadas o singulares y que su cartera es más diversificada que la de mínima varianza.

El segundo experimento es un rebalanceo mensual de 2008 a 2025: el primer día de cada mes, cada método calcula sus pesos con los dos años anteriores; entre rebalanceos la cartera se mantiene y los pesos derivan con los precios. Además de las métricas realizadas, la última columna da el Sharpe verdadero medio de los pesos elegidos: $w'\mu / \sqrt{w' \Sigma w}$ con el μ y la Σ verdaderos del universo. Separa lo que cada método eligió de la suerte que tuvo después.

V = 504                                   # dos años de ventana
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):
    """Compra y mantiene entre rebalanceos (los pesos derivan)."""
    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(), "a", fechas[-1].date())
pd.DataFrame(res).T.round(2)

2008-02-01 a 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

Cinco lecturas de la tabla:

• El orden de las volatilidades es el del Monte Carlo: 1/N 13,79 %, volatilidad inversa 8,51 %, ERC 6,25 %, máxima diversificación 6,06 %, HRP 4,45 % y mínima varianza 3,21 %, con caídas máximas que van del −46 % de la 1/N al −11 % de la mínima varianza. Los métodos hacen lo que dicen: menos riesgo cuanto más se concentran en lo poco volátil.
• Los Sharpes realizados (0,13 a 0,66) no dicen qué método es mejor. Con 18 años, el error típico de un Sharpe anual es de unos $1/\sqrt{18} \approx 0,24$. Los Sharpes verdaderos de lo que eligió cada método van de 0,39 (1/N) a 0,52 (máxima diversificación), y los realizados de los métodos con más bonos quedan por debajo porque los bonos del universo rindieron en la muestra menos que su esperanza (RF_T710: un 0,5 % de exceso medio en 2006-2025 frente a un 1,6 % verdadero).
• HERC tiene el mejor Sharpe realizado (0,66) y uno de los peores verdaderos (0,43). Pesó mucho la renta variable en años en que pagó bien. Es exactamente el tipo de resultado que el backtest premia y la verdad no respalda.
• La rotación separa los métodos más que el Sharpe. 1/N, volatilidad inversa, ERC y mínima varianza rotan entre un 35 % y un 39 % al año; máxima diversificación y HRP, un 166 % y un 179 %; HERC, un 399 %. La rotación cuesta (lección 4.1.5) y es el síntoma de un método inestable.
• La mínima varianza no fue la inestable. Su rotación (39 %) es la de la ERC. El tópico de que los optimizadores cuadráticos saltan y HRP no, en estos datos, se invierte.

¿De dónde sale la rotación de HRP? De la bisección por mitades. Un cambio pequeño en las correlaciones puede mover un activo de un extremo a otro de la cadena del árbol; entonces cambia qué activos caen a cada lado de un corte y el peso de un grupo entero salta:

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"rotación HRP {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 rotación HRP 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

El 1 de diciembre de 2015, con un mes más de datos en una ventana de dos años, HRP movió el 69 % de la cartera y la paridad de riesgo el 3 %. El bono a 3-7 años pasó del 21,1 % al 42,9 % y el de deuda global cubierta, del 24,2 % al 6,8 %. La figura sigue el peso del bono a 3-7 años en las dos carteras durante todo el periodo:

La paridad de riesgo se mueve despacio, porque cada peso depende de forma continua de la covarianza. HRP depende de un objeto discreto, el orden de las hojas, y salta. Las variantes que suavizan el árbol (otro enlace, orden óptimo de las hojas) cambian los saltos de sitio pero no los eliminan, porque todas siguen cortando una lista ordenada (el ejercicio 2 mide cuánto cambian los pesos con el árbol).

Qué método y cuándo

flowchart TB
  A{¿μ fiable?} -- sí --> B[Media-varianza]
  A -- no --> C{¿Objetivo?}
  C -- menos riesgo --> D[Mínima varianza]
  C -- repartir riesgo --> E{¿Redundantes?}
  E -- no --> F[ERC]
  E -- sí --> G[ERC por grupos]
  G --> H[o HRP]

• Mínima varianza cuando el objetivo es literalmente la menor volatilidad y se aceptan carteras concentradas. Con límites por activo (lección 4.1.1) se diversifica a la fuerza.
• ERC cuando se quiere repartir el riesgo y se controla el universo: pocos activos, distintos entre sí, o presupuestos por clase. Es estable y fácil de explicar. Necesita apalancamiento para rentar como una cartera con acciones.
• HRP cuando hay muchos activos redundantes, pocas observaciones o una covarianza que no se puede invertir. Hay que vigilar su rotación y fijar la configuración (enlace, orden de hojas) de forma explícita.
• HERC como idea (riesgo por grupos del árbol), con cuidado: el número de grupos cambia con la muestra y la cartera con él. Fijar $k$ o los grupos a mano (por clase de activo) da algo más estable: es una ERC por grupos.
• Ninguno elimina el riesgo de correlación. En una crisis las correlaciones suben juntas (en el régimen de estrés del universo salen todas de un único componente de aversión al riesgo), y una cartera calibrada con dos años de calma tiene más riesgo del que mide. Eso es la lección 4.2.4.

Errores típicos con estos métodos:

• Comparar métodos por el Sharpe de un backtest. Con 18 años, diferencias de 0,2 no se distinguen del ruido.
• Olvidar la rotación: un método con un 400 % de rotación anual paga costes que el backtest sin costes no ve.
• Leer «paridad de riesgo» como «riesgo repartido entre clases»: lo está entre activos.
• Usar los valores por defecto de una librería sin decirlo: el mismo nombre de método da pesos distintos.

Ejercicios

1. Presupuestos de riesgo por clase. Construye con la función erc una cartera en la que renta variable (con inmobiliario y alternativos), renta fija y materias primas aporten un tercio del riesgo cada una, repartido por igual dentro de cada grupo. Comprueba las contribuciones por clase y la volatilidad.

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

volatilidad verdadera: 5.48%
                 capital  riesgo
clase
materias_primas    0.196   0.333
renta_fija         0.661   0.333
renta_variable     0.142   0.333

2. Variantes de HRP. Con la covarianza muestral de 2024-2025 (S_m), calcula la diferencia máxima de pesos entre el HRP de la lección y: (a) el árbol del código del artículo, pasando a linkage la matriz cuadrada de distancias; (b) el enlace de Ward. ¿Cuál cambia más?

def hrp_arbol(S, Z):
    """Bisección recursiva sobre el orden de un árbol dado."""
    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 avisa de la matriz cuadrada
    warnings.simplefilter("ignore", ClusterWarning)
    articulo = hrp_arbol(S_m, linkage(d_m, "single"))
ward = hrp(S_m, metodo="ward")
print(f"árbol del artículo: {np.abs(articulo - base).max():.3f}")
print(f"enlace de Ward:     {np.abs(ward - base).max():.3f}")

árbol del artículo: 0.118
enlace de Ward:     0.112

3. Sin código: un gestor propone pasar de una 60/40 a una paridad de riesgo «porque diversifica mejor» y mantener la misma volatilidad. Explica con los números de la lección qué cambia en el capital, qué hace falta para igualar la volatilidad y qué riesgo nuevo aparece.