SHAP, dependencia parcial y estabilidad de las importancias.

Explicar el modelo, no el mundo

La lección 5.2.5 midió la importancia global de cada feature con cuatro métodos y llegó a una conclusión que aquí es el punto de partida: una explicación describe lo que hace el modelo, no cómo funciona el mercado. Si el modelo aprendió ruido, la explicación enseña ruido; si no aprendió algo, la explicación no lo puede enseñar. Esta lección va tres pasos más allá:

1. Explicaciones locales y de interacción: por qué el modelo predijo esto para esta acción esta semana, y qué parte de cada efecto depende de otra feature.
2. Forma de los efectos: la dependencia parcial y las curvas ICE, qué miden que SHAP no mide y al revés.
3. Estabilidad: si la explicación cambiaría con otra semilla, otra muestra u otro periodo. Una explicación que no sobrevive a eso no es una explicación.

flowchart TB
  Q[Pregunta] --> L["¿Por qué esta fila?"]
  Q --> F["¿Qué forma tiene?"]
  Q --> E["¿Es estable?"]
  L --> S[SHAP local]
  F --> P[PD e ICE]
  E --> B[Bootstrap]

El banco de pruebas es el panel del nivel 5: 400 acciones, 21 features semanales en rango normal (sin la fuga revision_analistas), dos dianas (la etiqueta real y el alfa verdadero, la parte predecible sin ruido) y el régimen de mercado de cada semana, que el generador publica en panel5_coef_verdad.csv. Los modelos se entrenan con 2012-2019 y se explican con filas de 2020-2025.

import numpy as np
import pandas as pd
from scipy import stats

f = pd.read_parquet("panel5_features.parquet")
e = pd.read_parquet("panel5_etiquetas.parquet")
v = pd.read_parquet("panel5_verdad.parquet")
cv = pd.read_csv("panel5_coef_verdad.csv", parse_dates=["fecha"])
d = (f.merge(e[["fecha", "ticker", "ret_sig_semana"]],
             on=["fecha", "ticker"])
      .merge(v[["fecha", "ticker", "alfa_verdadero"]],
             on=["fecha", "ticker"])
      .merge(cv[["fecha", "regimen"]], on="fecha")
      .dropna(subset=["ret_sig_semana"])
      .reset_index(drop=True))
FEATS = [c for c in f.columns
         if c not in ("fecha", "ticker", "sector",
                      "revision_analistas")]

sem = d.groupby("fecha")
def gauss(col):
    n = sem[col].transform("count")
    return stats.norm.ppf((sem[col].rank() - 0.5) / n)

X = pd.DataFrame({c: gauss(c) for c in FEATS}).fillna(0.0)
def menos_media(col):                 # en puntos básicos
    return ((d[col] - sem[col].transform("mean")) * 1e4).to_numpy()
y = menos_media("ret_sig_semana")     # la etiqueta real
a = menos_media("alfa_verdadero")     # la diana sin ruido
reg = d["regimen"].to_numpy()
tr = (d["fecha"] < "2020-01-01").to_numpy()   # 2012-2019
te = ~tr                                       # 2020-2025
print(X.shape, f"entrenamiento {tr.sum()} · prueba {te.sum()}")
periodo = np.where(tr, "2012-2019", "2020-2025")
print(pd.crosstab(periodo, reg, rownames=["semanas"],
                  colnames=[None]) // 400)

(291600, 21) entrenamiento 166400 · prueba 125200
           calma  estres  rebote
semanas
2012-2019    353      37      26
2020-2025    222      52      39

Hay 416 semanas de entrenamiento y 313 de prueba; los rebotes (tras cada episodio de estrés, cuando el momentum verdadero se invierte, lección 5.1.1) son 26 en el entrenamiento y 39 en la prueba. El modelo es un LightGBM moderado (300 árboles de 15 hojas, al menos 500 filas por hoja, submuestras de filas y columnas) como el de la lección 5.5.3.

import lightgbm as lgb

def lgbm(cols, obj, filas=tr, semilla=0, n=300, lr=0.05, w=None):
    return lgb.LGBMRegressor(
        n_estimators=n, learning_rate=lr, num_leaves=15,
        min_child_samples=500, subsample=0.5, subsample_freq=1,
        colsample_bytree=0.8, random_state=semilla, n_jobs=4,
        deterministic=True, force_row_wise=True, verbose=-1,
    ).fit(X.loc[filas, cols].to_numpy(), obj[filas],
          sample_weight=w)

def ic_semanal(pred, obj, filas):
    """IC de rangos medio de las semanas de `filas`."""
    r = pd.DataFrame({"p": pred, "y": obj[filas],
                      "s": d["fecha"].to_numpy()[filas]})
    c = r.groupby("s")[["p", "y"]].corr("spearman")
    return c.xs("p", level=1)["y"].mean()

m_a = lgbm(FEATS, a)          # aprende la parte predecible
m_y = lgbm(FEATS, y)          # lo que se puede hacer en la práctica
Xte = X.loc[te, FEATS].to_numpy()
for nombre, m in [("diana sin ruido", m_a), ("etiqueta", m_y)]:
    p = m.predict(Xte)
    print(f"{nombre:16s} IC con el alfa {ic_semanal(p, a, te):.3f}"
          f" · con la etiqueta {ic_semanal(p, y, te):.4f}")
print(f"oráculo          IC con la etiqueta "
      f"{ic_semanal(a[te], y, te):.4f}")

diana sin ruido  IC con el alfa 0.783 · con la etiqueta 0.0183
etiqueta         IC con el alfa 0.294 · con la etiqueta 0.0077
oráculo          IC con la etiqueta 0.0262

Con la diana sin ruido, el modelo alcanza un IC de 0,783 con el alfa de 2020-2025: ha aprendido buena parte de la estructura plantada. Con la etiqueta real se queda en 0,0077, frente al techo del oráculo (0,0262). Y el modelo de la diana sin ruido, evaluado contra la etiqueta, da 0,0183: lo que un modelo podría sacar si la señal no estuviera enterrada. Para ver qué recupera cada herramienta se usa sobre todo el primero; para ver cuánto se puede fiar uno de una explicación con datos reales, los dos.

SHAP local: por qué esta predicción

Los valores SHAP reparten la diferencia entre la predicción de una fila y el valor base (la predicción media en los datos del modelo) entre las features. La propiedad que lo hace útil es la aditividad: base más suma de valores igual a predicción, sin resto. shap.TreeExplainer la cumple de forma exacta en modelos de árboles (Lundberg et al., 2020).

import shap

rng = np.random.default_rng(1)
muestra = rng.choice(np.flatnonzero(te), 1000, replace=False)
Xm = X.loc[muestra, FEATS].to_numpy()
ex_a = shap.TreeExplainer(m_a)
sv = ex_a.shap_values(Xm)
suma = ex_a.expected_value + sv.sum(axis=1)
print("máx. |base + suma SHAP − predicción|:",
      f"{np.abs(suma - m_a.predict(Xm)).max():.1e}")
i = 0                                      # una acción-semana
fila = pd.DataFrame({"valor (rango normal)": Xm[i],
                     "shap (pb)": sv[i]}, index=FEATS)
print(d.loc[muestra[i], ["fecha", "ticker"]].to_dict(),
      f"base {ex_a.expected_value:.2f} pb · predicción "
      f"{m_a.predict(Xm[[i]])[0]:.2f} pb · alfa verdadero "
      f"{a[muestra[i]]:.2f} pb")
fila.reindex(fila["shap (pb)"].abs().sort_values(
    ascending=False).index).head(6).round(2)

máx. |base + suma SHAP − predicción|: 9.6e-14
{'fecha': Timestamp('2020-05-29 00:00:00'), 'ticker': 'MACE'} base 0.01 pb · predicción 15.65 pb · alfa verdadero 11.34 pb
              valor (rango normal)  shap (pb)
valor_bp                      2.08       7.11
vol_63d                      -0.78       2.82
ey                            1.68       2.14
mom_12_1                      0.32       1.25
ret_1s                       -0.12       1.16
dist_max_52s                  1.68      -0.90

La diferencia máxima entre base más suma y predicción es de $10^{-13}$. Para la acción MACE el 2020-05-29, el modelo predice 15,65 pb de rendimiento extra para la semana siguiente (el alfa verdadero era 11,34). La mayor parte viene de valor_bp (7,11 pb, con un rango normal de 2,08: una acción muy barata por patrimonio), seguida de una volatilidad baja (vol_63d en −0,78 aporta 2,82) y de ey (2,14), que no tiene efecto propio y carga con parte del value por su correlación con valor_bp. Es la explicación que pide un gestor («¿por qué compramos esto?») y es exactamente la del modelo, incluidos sus préstamos entre features correlacionadas.

Valores de interacción

El valor SHAP de una feature mezcla su efecto propio con su parte de las interacciones. Los valores de interacción SHAP (Lundberg et al., 2020) lo separan en una matriz por fila: la diagonal $\Phi_{jj}$ es el efecto principal de la feature $j$ y cada casilla fuera de la diagonal $\Phi_{jk}$ es la mitad de su interacción con $k$, de modo que $\phi_j = \Phi_{jj} + \sum_{k \ne j} \Phi_{jk}$.

El generador planta la reversión como $c_{rev}\cdot g(r)\cdot(1 - 0{,}5\, g(v))$, con $r$ = ret_1s y $v$ = volumen_rel. El efecto principal de ret_1s tiene pendiente $c_{rev}$ y el término cruzado, $-0{,}5\,c_{rev}$ por unidad de $g(r)\,g(v)$, que se reparte a partes iguales entre $\Phi_{rv}$ y $\Phi_{vr}$: $-0{,}25\,c_{rev}$ en cada una.

si = ex_a.shap_interaction_values(Xm[:600])
inter = np.abs(si).mean(axis=0)
np.fill_diagonal(inter, 0)
pares = pd.Series(
    {f"{FEATS[j]} × {FEATS[k]}": 2 * inter[j, k]
     for j in range(len(FEATS)) for k in range(j + 1, len(FEATS))})
print(pares.sort_values(ascending=False).head(4).round(2)
      .to_string())

j, k = FEATS.index("ret_1s"), FEATS.index("volumen_rel")
r, z = Xm[:600, j], Xm[:600, k]
c_rev = cv.loc[cv["fecha"] < "2020-01-01", "c_rev"].mean() * 1e4
print(f"c_rev verdadero medio 2012-2019: {c_rev:.2f} pb")
print(f"efecto principal de ret_1s: pendiente "
      f"{np.polyfit(r, si[:, j, j], 1)[0]:.2f} (verdad {c_rev:.2f})")
print(f"interacción con volumen_rel (cada mitad): "
      f"{np.polyfit(r * z, si[:, j, k], 1)[0]:.2f} "
      f"(verdad {-0.25 * c_rev:.2f})")

ret_1s × volumen_rel        1.47
ret_1s × vol_63d            0.39
mom_12_1 × valor_bp         0.35
sue × dias_desde_informe    0.32
c_rev verdadero medio 2012-2019: -5.46 pb
efecto principal de ret_1s: pendiente -4.74 (verdad -5.46)
interacción con volumen_rel (cada mitad): 0.93 (verdad 1.37)

El par con más interacción, con diferencia, es ret_1s × volumen_rel (1,47 pb de media en valor absoluto, contando las dos casillas). El cuarto, sue × dias_desde_informe, también está plantado (el SUE decae con los días desde el informe). Los otros dos no: ret_1s × vol_63d y mom_12_1 × valor_bp son interacciones que el modelo encontró en el entrenamiento, quizá porque la reversión y el momentum cambian con el régimen y la volatilidad o el value se mueven con él, quizá ruido. Una interacción grande en el modelo es una hipótesis, no un descubrimiento.

Los tamaños se pueden comparar con la verdad. El efecto principal de ret_1s tiene pendiente −4,74 frente al c_rev medio de 2012-2019 (−5,46) y cada mitad de la interacción, 0,93 frente a 1,37. El modelo aprendió la reversión con un 87 % de su tamaño y la interacción con un 68 %: la regularización (hojas de al menos 500 filas, 300 árboles con tasa de aprendizaje 0,05) encoge más lo que es más difícil de ver.

Dependencia parcial e ICE

La dependencia parcial (Friedman, 2001) responde a otra pregunta: cómo cambia la predicción media al mover una feature con las demás como están. Para la feature $j$ en el valor $x$:

Cada sumando es una curva ICE (individual conditional expectation, Goldstein et al., 2015): la predicción de la fila $i$ cuando solo se mueve su feature $j$. La dependencia parcial es la media de las curvas ICE; si las curvas son paralelas, el efecto no depende de las demás features, y si se abren en abanico, hay interacción. sklearn.inspection.partial_dependence las calcula con cualquier modelo con predict (aquí, por fuerza bruta: 21 puntos de la rejilla × 1.000 filas).

from sklearn.inspection import partial_dependence

pdep = partial_dependence(m_a, Xm, [j], kind="both",
                          grid_resolution=21,
                          percentiles=(0.02, 0.98))
rejilla = pdep["grid_values"][0]
print("PD de ret_1s (pb), de −2 a +2:",
      np.round(pdep["average"][0][::5], 1))
pend = np.array([np.polyfit(rejilla, c, 1)[0]
                 for c in pdep["individual"][0]])    # una por ICE
vol = Xm[:, k]
tercio = pd.cut(vol, [-np.inf, -0.43, 0.43, np.inf],
                labels=["bajo", "medio", "alto"])
ice = pd.Series(pend).groupby(tercio, observed=True).mean()
media_v = pd.Series(vol).groupby(tercio, observed=True).mean()
verdad = 1 - 0.5 * media_v
print(f"corr(pendiente ICE, volumen_rel) = "
      f"{np.corrcoef(pend, vol)[0, 1]:.2f}")
pd.DataFrame({"pendiente ICE (pb)": ice,
              "ICE relativa": ice / ice["medio"],
              "verdad relativa": verdad / verdad["medio"]}).round(2)

PD de ret_1s (pb), de −2 a +2: [ 10.    5.1  -0.4  -6.1 -10.7]
corr(pendiente ICE, volumen_rel) = 0.96
       pendiente ICE (pb)  ICE relativa  verdad relativa
bajo                -7.58          1.48             1.52
medio               -5.12          1.00             1.00
alto                -2.46          0.48             0.43

La dependencia parcial de ret_1s baja de +10 pb a −10,7 pb entre los percentiles 2 y 98: la reversión media. Las curvas ICE no son paralelas: la correlación entre la pendiente de cada curva y el volumen_rel de su fila es 0,96. Agrupadas por tercios de volumen, las pendientes son −7,58, −5,12 y −2,46 pb por unidad; relativas al tercio medio, 1,48 y 0,48, frente a 1,52 y 0,43 del generador. La figura las pone en pb junto a la verdad con el c_rev medio del entrenamiento.

Es la diferencia con la lección 5.2.5, donde las pendientes SHAP del mismo efecto por tercios daban 1,25 y 0,72 frente al valor exacto de Shapley, 1,27 y 0,72. SHAP atribuye: reparte la interacción entre las dos features. El ICE agrupado condiciona: dice cómo responde la predicción a ret_1s cuando el volumen está bajo, con la interacción entera. Ninguno está mal; responden a preguntas distintas.

El momentum según el régimen: lo que el modelo no ve

El momentum del panel cambia con el régimen de mercado: 5,0 pb por unidad en calma, 1,7 en estrés y −7,1 en rebote (medias de 2012-2019). Un modelo solo puede reproducir eso si ve el régimen. Tres versiones con la diana sin ruido: sin estado (las 21 features), con un estado observable (la media transversal de ret_1m, ret_3m y dist_max_52s: cómo le ha ido al mercado el último mes, el último trimestre y cuánto cae desde máximos) y con el régimen verdadero como feature, algo que en la realidad nadie tiene. En cada régimen de 2020-2025 se mide la pendiente del momentum con SHAP y con la dependencia parcial calculada sobre filas de ese régimen.

for c in ["ret_1m", "ret_3m", "dist_max_52s"]:   # estado observable
    X[f"mercado_{c}"] = sem[c].transform("mean")
X["rebote"] = (reg == "rebote").astype(float)    # estado verdadero
X["estres"] = (reg == "estres").astype(float)
modelos = {
    "sin estado": FEATS,
    "estado observable": FEATS + ["mercado_ret_1m", "mercado_ret_3m",
                                  "mercado_dist_max_52s"],
    "régimen verdadero": FEATS + ["rebote", "estres"],
}
rng = np.random.default_rng(2)
jm = FEATS.index("mom_12_1")
filas_t = {}
for nombre, cols in modelos.items():
    m = m_a if nombre == "sin estado" else lgbm(cols, a)
    for rg in ["calma", "estres", "rebote"]:
        s = rng.choice(np.flatnonzero(te & (reg == rg)), 800,
                       replace=False)
        Xs = X.loc[s, cols].to_numpy()
        sv_s = shap.TreeExplainer(m).shap_values(Xs)
        pdr = partial_dependence(m, Xs, [jm], grid_resolution=11,
                                 percentiles=(0.05, 0.95))
        filas_t[(nombre, rg)] = {
            "shap": np.polyfit(Xs[:, jm], sv_s[:, jm], 1)[0],
            "pd": np.polyfit(pdr["grid_values"][0],
                             pdr["average"][0], 1)[0]}
t = pd.DataFrame(filas_t).T.unstack().loc[list(modelos)]
verdad = cv.groupby([cv["fecha"] >= "2020-01-01", "regimen"]
                    )["c_mom"].mean().unstack() * 1e4
t.loc["verdad 2012-2019"] = list(verdad.loc[False]) * 2
t.loc["verdad 2020-2025"] = list(verdad.loc[True]) * 2
t.round(2)

                   shap                  pd
                  calma estres rebote calma estres rebote
sin estado         3.64   3.67   3.61  3.84   3.87   3.88
estado observable  3.84   1.96   2.32  4.27   0.40   1.26
régimen verdadero  4.00   2.83  -0.94  4.58   2.00  -6.47
verdad 2012-2019   5.01   1.74  -7.07  5.01   1.74  -7.07
verdad 2020-2025   4.23   2.22 -12.68  4.23   2.22 -12.68

Sin estado, las dos herramientas dan la misma pendiente en los tres regímenes (3,6-3,9 pb): el modelo aprendió un momentum medio y lo aplica siempre. La explicación es fiel al modelo y el modelo se equivoca en los rebotes, donde la verdad es −12,68 en la prueba. Con el régimen verdadero, la dependencia parcial en filas de rebote da −6,47, muy cerca de los −7,07 del entrenamiento (los rebotes de 2020-2025 son más violentos que los que vio el modelo). SHAP, en cambio, da −0,94: la interacción momentum × régimen se reparte, y la mitad va a la feature rebote. Leer el SHAP del momentum como «el momentum casi no importa en los rebotes» sería un error.

El estado observable no basta (+1,26 en rebote con la dependencia parcial). La rentabilidad del mercado en el último mes y en el último trimestre, o la caída desde máximos, no separan un rebote de un estrés a tiempo, y con 26 semanas de rebote en el entrenamiento el modelo tiene muy pocos ejemplos. Es la misma conclusión de la lección 5.10.2 vista desde la explicación: lo que el modelo no ve, ni lo aprende ni lo explica, y una explicación limpia de un modelo sin el estado correcto da una imagen tranquilizadora y falsa.

Estabilidad de las importancias

Una importancia es una estimación y tiene error. Tres fuentes de variación, de menos a más relevante:

• Semilla: la misma muestra con otro azar en el entrenamiento (submuestras de filas y columnas).
• Bootstrap de semanas: otra muestra del mismo proceso, remuestreando semanas enteras con reemplazo (las filas de una semana comparten el mercado de esa semana; remuestrear filas sueltas subestimaría la variación).
• Ventana temporal: dos años de datos distintos. Además del muestreo, los efectos verdaderos cambian (el generador los hace variar con un AR(1) semanal y con el régimen).

Para cada réplica se entrena un LightGBM algo más ligero (150 árboles, tasa 0,1) y se mide el SHAP medio en valor absoluto en las mismas 1.000 filas de prueba. La estabilidad es la correlación de rangos media entre réplicas.

from scipy.stats import spearmanr

def importancia(m):
    return np.abs(shap.TreeExplainer(m).shap_values(Xm)).mean(0)

fechas = d["fecha"].to_numpy()
semanas = np.unique(fechas[tr])
ventanas = [("2012", "2013"), ("2014", "2015"),
            ("2016", "2017"), ("2018", "2019")]
imp = {}
for diana, obj in [("alfa", a), ("etiqueta", y)]:
    for s in range(4):                 # misma muestra, otra semilla
        imp[(diana, "semilla", s)] = importancia(
            lgbm(FEATS, obj, semilla=s, n=150, lr=0.1))
    for i0, i1 in ventanas:            # dos años de datos cada una
        filas = (d["fecha"] >= f"{i0}-01-01") & (
            d["fecha"] <= f"{i1}-12-31")
        imp[(diana, "ventana", i0)] = importancia(
            lgbm(FEATS, obj, filas.to_numpy(), n=150, lr=0.1))
    for b in range(4):                 # semanas con reemplazo
        veces = pd.Series(np.random.default_rng(b).choice(
            semanas, len(semanas))).value_counts()
        w = pd.Series(fechas[tr]).map(veces).fillna(0).to_numpy()
        imp[(diana, "bootstrap", b)] = importancia(
            lgbm(FEATS, obj, n=150, lr=0.1, w=w))
imp = pd.DataFrame(imp, index=FEATS)

def rho_medio(cols):
    c = spearmanr(imp[cols]).statistic
    return c[np.triu_indices(len(cols), 1)].mean()

pd.DataFrame({diana: {p: rho_medio([c for c in imp
                                     if c[:2] == (diana, p)])
                      for p in ["semilla", "ventana", "bootstrap"]}
              for diana in ["alfa", "etiqueta"]}).round(2)

           alfa  etiqueta
semilla    0.99      0.81
ventana    0.90      0.31
bootstrap  0.98      0.51

Con la diana sin ruido, la explicación es estable con todo: 0,99 entre semillas, 0,98 entre bootstrap y 0,90 entre ventanas, la más baja porque los efectos verdaderos cambian de un periodo a otro. Con la etiqueta real se desmorona: 0,81 entre semillas, 0,51 entre bootstrap y 0,31 entre ventanas. La estabilidad entre semillas es engañosa, porque el modelo lee el mismo ruido de la misma muestra de forma parecida; la incertidumbre que importa es la de los datos.

rangos = imp.rank(ascending=False)
sel = ["ret_1s", "mom_12_1", "vol_63d", "valor_bp", "sue",
       "sentimiento", "ret_1m", "ruido_1", "ruido_4"]
replicas = [c for c in imp if c[1] != "semilla"]
res = {}
for diana in ["alfa", "etiqueta"]:
    rr = rangos[[c for c in replicas if c[0] == diana]]
    res[(diana, "mín")] = rr.min(axis=1)
    res[(diana, "máx")] = rr.max(axis=1)
pd.DataFrame(res).loc[sel].astype(int)

            alfa     etiqueta
             mín máx      mín máx
ret_1s         1   1        1   5
mom_12_1       2   4        1   8
vol_63d        2   5        1   9
valor_bp       3   6        2  18
sue            3   5        2  18
sentimiento    4   6       14  21
ret_1m        11  20        1  17
ruido_1       18  21        8  20
ruido_4       15  21        6  19

La figura de la cabecera pone los rangos de puestos en las ocho réplicas de datos (ventanas y bootstrap). Con la diana sin ruido, ret_1s es siempre la primera y ruido_1 está siempre entre las cuatro últimas. Con la etiqueta real, valor_bp y sue, con efecto plantado, van del 2.º al 18.º puesto; ret_1m, sin efecto propio, del 1.º al 17.º; ruido_4 llega al 6.º y ruido_1 al 8.º. Solo ret_1s se mantiene arriba (del 1.º al 5.º), y sentimiento, con efecto pero con un 32 % de valores vacíos rellenados con cero, no pasa del 14.º. Un informe que dijera «la tercera feature más importante es X» con estos datos estaría informando de una realización del ruido.

Cuando las explicaciones engañan: features correlacionadas

La primera fuente de inestabilidad es el reparto dentro de grupos de features correlacionadas: el modelo puede usar cualquiera de ellas y el reparto cambia de una réplica a otra aunque el total no cambie. Se comprueba con los grupos por correlación de la lección 5.2.5.

grupos = {"reversión": ["ret_1s", "ret_1m", "rsi_14"],
          "medio plazo": ["ret_3m", "dist_max_52s"],
          "momentum": ["mom_12_1", "ruido_4"],
          "volatilidad": ["vol_21d", "vol_63d"],
          "value-calidad": ["valor_bp", "ey", "roe"]}
etiqueta = {c: g for g, cs in grupos.items() for c in cs}
imp_g = imp.groupby([etiqueta.get(c, c) for c in FEATS]).sum()
boot = [c for c in imp if c[:2] == ("etiqueta", "bootstrap")]
cg = spearmanr(imp_g[boot]).statistic
print(f"etiqueta, bootstrap: rho medio por features "
      f"{rho_medio(boot):.2f} · por grupos "
      f"{cg[np.triu_indices(4, 1)].mean():.2f}")
cuota = imp.loc[grupos["reversión"], boot]
(cuota / cuota.sum()).round(2)

etiqueta, bootstrap: rho medio por features 0.51 · por grupos 0.83
        etiqueta
       bootstrap
               0     1     2     3
ret_1s      0.44  0.45  0.51  0.39
ret_1m      0.27  0.27  0.27  0.39
rsi_14      0.29  0.27  0.22  0.22

Sumando el SHAP de cada grupo, la correlación de rangos entre bootstrap pasa de 0,51 a 0,83. Dentro del grupo de la reversión, la cuota de ret_1s va del 39 % al 51 % según la réplica y ret_1m llega al 39 % en una de ellas: el grupo es estable, el reparto no.

La forma más extrema es la sustitución: si se quita una feature, su sustituta ocupa el sitio.

sin_mom = [c for c in FEATS if c != "mom_12_1"]
m_sin = lgbm(sin_mom, a)
def shap_medio(m, cols):
    s = shap.TreeExplainer(m).shap_values(X.loc[muestra, cols]
                                          .to_numpy())
    return pd.Series(np.abs(s).mean(0), index=cols)

comp = pd.DataFrame({"con mom_12_1": shap_medio(m_a, FEATS),
                     "sin mom_12_1": shap_medio(m_sin, sin_mom)})
p_sin = m_sin.predict(X.loc[te, sin_mom].to_numpy())
print("IC con el alfa en 2020-2025:",
      f"con {ic_semanal(m_a.predict(Xte), a, te):.3f} · sin "
      f"{ic_semanal(p_sin, a, te):.3f}")
print("corr(ruido_4, mom_12_1) =",
      round(X["ruido_4"].corr(X["mom_12_1"]), 2),
      "· corr(dist_max_52s, mom_12_1) =",
      round(X["dist_max_52s"].corr(X["mom_12_1"]), 2))
comp.loc[["mom_12_1", "dist_max_52s", "ruido_4", "ret_3m",
          "ret_1m"]].round(2)

IC con el alfa en 2020-2025: con 0.783 · sin 0.777
corr(ruido_4, mom_12_1) = 0.6 · corr(dist_max_52s, mom_12_1) = 0.64
              con mom_12_1  sin mom_12_1
mom_12_1              3.08           NaN
dist_max_52s          0.32          1.74
ruido_4               0.03          0.94
ret_3m                0.09          0.17
ret_1m                0.25          0.34

Sin mom_12_1, el SHAP medio de ruido_4 (0,6 × momentum más ruido, sin efecto propio) pasa de 0,03 a 0,94 pb y el de dist_max_52s (correlación 0,64), de 0,32 a 1,74. El IC del modelo apenas cambia (0,783 frente a 0,777): la información del momentum sigue dentro, repartida entre sus vecinas. Para quien mire solo la explicación, el segundo modelo «demuestra» que la distancia al máximo de 52 semanas es una de las features más importantes del mercado.

Las mismas trampas, resumidas:

Cómo presentar una explicación

flowchart TB
  M[Modelo validado] --> R[Réplicas]
  R --> G[Grupos]
  G --> E{¿Estable?}
  E -- no --> N[No interpretar]
  E -- sí --> F[Forma: PD, ICE]
  F --> H[Hipótesis]

• Primero el modelo, después la explicación. Explicar un modelo que no supera la validación honesta (lecciones 5.4.2 y 5.4.4) es explicar ruido.
• Réplicas siempre. Bootstrap de bloques de tiempo y ventanas; se informa del rango de cada importancia o de su puesto, no de un número.
• Grupos antes que features sueltas, con controles de ruido dentro del modelo (lección 5.2.5).
• SHAP para atribuir, dependencia parcial e ICE para la forma, y valores de interacción cuando la pregunta es si dos features se refuerzan.
• Una explicación es una hipótesis. Si sugiere un mecanismo (la reversión es más fuerte con poco volumen), el mecanismo se contrasta aparte, con datos que no se usaron para encontrarlo.

Ejercicios

1. TreeExplainer calcula por defecto los valores SHAP siguiendo los caminos del árbol (tree_path_dependent), que respetan en parte las correlaciones de los datos de entrenamiento. Con feature_perturbation="interventional" y un fondo de 100 filas, rompe esas correlaciones. Compara el SHAP medio de las dos versiones en 500 filas. ¿Dónde difieren?

fondo = X.loc[np.random.default_rng(3).choice(
    np.flatnonzero(tr), 100, replace=False), FEATS].to_numpy()
sv_int = shap.TreeExplainer(
    m_a, data=fondo, feature_perturbation="interventional"
).shap_values(Xm[:500])
sv_cam = ex_a.shap_values(Xm[:500])
dif = pd.DataFrame({"intervencional": np.abs(sv_int).mean(0),
                    "caminos": np.abs(sv_cam).mean(0)}, index=FEATS)
dif["diferencia"] = dif["intervencional"] - dif["caminos"]
dif.sort_values("caminos", ascending=False).head(11).round(2)

              intervencional  caminos  diferencia
ret_1s                  4.63     4.37        0.26
mom_12_1                3.22     3.06        0.16
vol_63d                 2.56     2.52        0.04
sue                     2.42     2.45       -0.04
valor_bp                2.35     2.29        0.06
sentimiento             1.84     1.82        0.02
ey                      0.88     0.85        0.03
volumen_rel             0.82     0.77        0.05
roe                     0.49     0.43        0.05
dist_max_52s            0.46     0.33        0.13
ret_1m                  0.33     0.25        0.08

2. Explica la misma acción-semana (la fila MACE del 2020-05-29) con cuatro LightGBM de semillas distintas, entrenados con la diana sin ruido y con la etiqueta real. ¿Coinciden las tres features con más SHAP?

fila = Xm[[0]]
for diana, obj in [("alfa", a), ("etiqueta", y)]:
    print(diana)
    for s in range(4):
        m = lgbm(FEATS, obj, semilla=s, n=150, lr=0.1)
        sv1 = pd.Series(shap.TreeExplainer(m).shap_values(fila)[0],
                        index=FEATS)
        top = sv1.abs().sort_values(ascending=False).index[:3]
        print(f"  semilla {s}: pred {m.predict(fila)[0]:6.2f} pb · "
              + ", ".join(f"{c} {sv1[c]:+.1f}" for c in top))

alfa
  semilla 0: pred  15.54 pb · valor_bp +7.1, vol_63d +2.6, ey +2.3
  semilla 1: pred  15.22 pb · valor_bp +7.3, vol_63d +2.9, ey +2.2
  semilla 2: pred  14.32 pb · valor_bp +7.1, vol_63d +2.6, ey +2.1
  semilla 3: pred  15.04 pb · valor_bp +7.1, vol_63d +2.8, ey +2.0
etiqueta
  semilla 0: pred  24.23 pb · valor_bp +15.5, ret_1s +6.1, volumen_rel +4.5
  semilla 1: pred   2.87 pb · valor_bp +9.1, dist_max_52s -5.1, dias_desde_informe +2.3
  semilla 2: pred  -0.15 pb · valor_bp +9.4, rsi_14 -6.8, ruido_5 +5.5
  semilla 3: pred  -0.08 pb · dias_desde_informe +6.5, vol_63d -4.8, volumen_rel -4.3

3. Un comité de inversión pide «las tres razones por las que el modelo compra» cada semana. Con lo visto en la lección, ¿qué le entregarías y con qué advertencias?