SHAP, partial dependence and the stability of feature importances.

Explaining the model, not the world

Lesson 5.2.5 measured the global importance of each feature with four methods and reached a conclusion that is the starting point here: an explanation describes what the model does, not how the market works. If the model learned noise, the explanation shows noise; if it did not learn something, the explanation cannot show it. This lesson goes three steps further:

1. Local and interaction explanations: why the model predicted this for this stock this week, and what part of each effect depends on another feature.
2. Shape of the effects: partial dependence and ICE curves, what they measure that SHAP does not and vice versa.
3. Stability: whether the explanation would change with another seed, another sample or another period. An explanation that does not survive that is not an explanation.

flowchart TB
  Q[Question] --> L["Why this row?"]
  Q --> F["What shape does it have?"]
  Q --> E["Is it stable?"]
  L --> S[Local SHAP]
  F --> P[PD and ICE]
  E --> B[Bootstrap]

The test bench is the level 5 panel: 400 stocks, 21 weekly features in normal rank (without the leak revision_analistas), two targets (the real label and the true alpha, the predictable part without noise) and the market regime of each week, which the generator publishes in panel5_coef_verdad.csv. The models are trained on 2012-2019 and explained with rows from 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):                 # in basis points
    return ((d[col] - sem[col].transform("mean")) * 1e4).to_numpy()
y = menos_media("ret_sig_semana")     # the real label
a = menos_media("alfa_verdadero")     # the noise-free target
reg = d["regimen"].to_numpy()
tr = (d["fecha"] < "2020-01-01").to_numpy()   # 2012-2019
te = ~tr                                       # 2020-2025
print(X.shape, f"training {tr.sum()} · test {te.sum()}")
periodo = np.where(tr, "2012-2019", "2020-2025")
print(pd.crosstab(periodo, reg, rownames=["weeks"],
                  colnames=[None]) // 400)

(291600, 21) training 166400 · test 125200
           calma  estres  rebote
weeks
2012-2019    353      37      26
2020-2025    222      52      39

There are 416 training weeks and 313 test weeks; the rebounds (after each stress episode, when true momentum reverses, lesson 5.1.1) are 26 in training and 39 in the test. The model is a moderate LightGBM (300 trees of 15 leaves, at least 500 rows per leaf, row and column subsampling) like the one in lesson 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):
    """Mean rank IC over the weeks of `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)          # learns the predictable part
m_y = lgbm(FEATS, y)          # what can be done in practice
Xte = X.loc[te, FEATS].to_numpy()
for nombre, m in [("noise-free", m_a), ("label", m_y)]:
    p = m.predict(Xte)
    print(f"{nombre:16s} IC with the alpha {ic_semanal(p, a, te):.3f}"
          f" · with the label {ic_semanal(p, y, te):.4f}")
print(f"oracle           IC with the label "
      f"{ic_semanal(a[te], y, te):.4f}")

noise-free       IC with the alpha 0.783 · with the label 0.0183
label            IC with the alpha 0.294 · with the label 0.0077
oracle           IC with the label 0.0262

With the noise-free target, the model reaches an IC of 0.783 with the 2020-2025 alpha: it has learned a good part of the planted structure. With the real label it stays at 0.0077, against the oracle's ceiling (0.0262). And the noise-free target model, evaluated against the label, gives 0.0183: what a model could extract if the signal were not buried. To see what each tool recovers, the first model is used most; to see how far an explanation with real data can be trusted, both.

Local SHAP: why this prediction

SHAP values split the difference between the prediction for a row and the base value (the mean prediction on the model's data) among the features. The property that makes it useful is additivity: base plus the sum of values equals the prediction, with no remainder. shap.TreeExplainer satisfies it exactly for tree models (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("max. |base + SHAP sum − prediction|:",
      f"{np.abs(suma - m_a.predict(Xm)).max():.1e}")
i = 0                                      # one stock-week
fila = pd.DataFrame({"value (normal rank)": Xm[i],
                     "shap (bp)": sv[i]}, index=FEATS)
print(d.loc[muestra[i], ["fecha", "ticker"]].to_dict(),
      f"base {ex_a.expected_value:.2f} bp · prediction "
      f"{m_a.predict(Xm[[i]])[0]:.2f} bp · true alpha "
      f"{a[muestra[i]]:.2f} bp")
fila.reindex(fila["shap (bp)"].abs().sort_values(
    ascending=False).index).head(6).round(2)

max. |base + SHAP sum − prediction|: 9.6e-14
{'fecha': Timestamp('2020-05-29 00:00:00'), 'ticker': 'MACE'} base 0.01 bp · prediction 15.65 bp · true alpha 11.34 bp
              value (normal rank)  shap (bp)
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

The maximum difference between base plus sum and the prediction is $10^{-13}$. For the stock MACE on 2020-05-29, the model predicts 15.65 bp of excess return for the following week (the true alpha was 11.34). Most of it comes from valor_bp (7.11 bp, with a normal rank of 2.08: a stock that is very cheap relative to book value), followed by low volatility (vol_63d at −0.78 contributes 2.82) and by ey (2.14), which has no effect of its own and carries part of the value effect through its correlation with valor_bp. It is the explanation a portfolio manager asks for (“why are we buying this?”) and it is exactly the model's, including its borrowing between correlated features.

Interaction values

The SHAP value of a feature mixes its own effect with its share of the interactions. SHAP interaction values (Lundberg et al., 2020) separate it into a matrix per row: the diagonal $\Phi_{jj}$ is the main effect of feature $j$ and each off-diagonal cell $\Phi_{jk}$ is half of its interaction with $k$, so that $\phi_j = \Phi_{jj} + \sum_{k \ne j} \Phi_{jk}$.

The generator plants reversal as $c_{rev}\cdot g(r)\cdot(1 - 0.5\, g(v))$, with $r$ = ret_1s and $v$ = volumen_rel. The main effect of ret_1s has slope $c_{rev}$ and the cross term, $-0.5\,c_{rev}$ per unit of $g(r)\,g(v)$, is split equally between $\Phi_{rv}$ and $\Phi_{vr}$: $-0.25\,c_{rev}$ in each.

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"mean true c_rev 2012-2019: {c_rev:.2f} bp")
print(f"main effect of ret_1s: slope "
      f"{np.polyfit(r, si[:, j, j], 1)[0]:.2f} (truth {c_rev:.2f})")
print(f"interaction with volumen_rel (each half): "
      f"{np.polyfit(r * z, si[:, j, k], 1)[0]:.2f} "
      f"(truth {-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
mean true c_rev 2012-2019: -5.46 bp
main effect of ret_1s: slope -4.74 (truth -5.46)
interaction with volumen_rel (each half): 0.93 (truth 1.37)

The pair with the most interaction, by far, is ret_1s × volumen_rel (1.47 bp on average in absolute value, counting both cells). The fourth, sue × dias_desde_informe, is also planted (the SUE decays with the days since the report). The other two are not: ret_1s × vol_63d and mom_12_1 × valor_bp are interactions the model found in training, perhaps because reversal and momentum change with the regime and volatility or value move with it, perhaps noise. A large interaction in the model is a hypothesis, not a discovery.

The sizes can be compared with the truth. The main effect of ret_1s has a slope of −4.74 against the mean c_rev of 2012-2019 (−5.46), and each half of the interaction, 0.93 against 1.37. The model learned reversal at 87% of its size and the interaction at 68%: regularization (leaves of at least 500 rows, 300 trees with a learning rate of 0.05) shrinks more what is harder to see.

Partial dependence and ICE

Partial dependence (Friedman, 2001) answers another question: how the mean prediction changes when one feature is moved with the others left as they are. For feature $j$ at value $x$:

Each term is an ICE curve (individual conditional expectation, Goldstein et al., 2015): the prediction for row $i$ when only its feature $j$ is moved. Partial dependence is the mean of the ICE curves; if the curves are parallel, the effect does not depend on the other features, and if they fan out, there is interaction. sklearn.inspection.partial_dependence computes them with any model that has predict (here, by brute force: 21 grid points × 1,000 rows).

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 of ret_1s (bp), from −2 to +2:",
      np.round(pdep["average"][0][::5], 1))
pend = np.array([np.polyfit(rejilla, c, 1)[0]
                 for c in pdep["individual"][0]])    # one per ICE
vol = Xm[:, k]
tercio = pd.cut(vol, [-np.inf, -0.43, 0.43, np.inf],
                labels=["low", "mid", "high"])
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(ICE slope, volumen_rel) = "
      f"{np.corrcoef(pend, vol)[0, 1]:.2f}")
pd.DataFrame({"ICE slope (bp)": ice,
              "relative ICE": ice / ice["mid"],
              "relative truth": verdad / verdad["mid"]}).round(2)

PD of ret_1s (bp), from −2 to +2: [ 10.    5.1  -0.4  -6.1 -10.7]
corr(ICE slope, volumen_rel) = 0.96
      ICE slope (bp)  relative ICE  relative truth
low            -7.58          1.48            1.52
mid            -5.12          1.00            1.00
high           -2.46          0.48            0.43

The partial dependence of ret_1s falls from +10 bp to −10.7 bp between the 2nd and 98th percentiles: the mean reversal. The ICE curves are not parallel: the correlation between the slope of each curve and the volumen_rel of its row is 0.96. Grouped by volume terciles, the slopes are −7.58, −5.12 and −2.46 bp per unit; relative to the middle tercile, 1.48 and 0.48, against 1.52 and 0.43 from the generator. The figure shows them in bp alongside the truth with the mean c_rev of the training period.

This is the difference with lesson 5.2.5, where the SHAP slopes of the same effect by terciles gave 1.25 and 0.72 against the exact Shapley value, 1.27 and 0.72. SHAP attributes: it splits the interaction between the two features. Grouped ICE conditions: it says how the prediction responds to ret_1s when volume is low, with the whole interaction. Neither is wrong; they answer different questions.

Momentum by regime: what the model does not see

Momentum in the panel changes with the market regime: 5.0 bp per unit in calm, 1.7 in stress and −7.1 in rebound (2012-2019 means). A model can only reproduce that if it sees the regime. Three versions with the noise-free target: no state (the 21 features), with an observable state (the cross-sectional mean of ret_1m, ret_3m and dist_max_52s: how the market has done over the last month, the last quarter and how far it is below its highs) and with the true regime as a feature, something nobody has in reality. In each 2020-2025 regime the slope of momentum is measured with SHAP and with the partial dependence computed on rows of that regime.

for c in ["ret_1m", "ret_3m", "dist_max_52s"]:   # observable state
    X[f"mercado_{c}"] = sem[c].transform("mean")
X["rebote"] = (reg == "rebote").astype(float)    # true state
X["estres"] = (reg == "estres").astype(float)
modelos = {
    "no state": FEATS,
    "observable state": FEATS + ["mercado_ret_1m", "mercado_ret_3m",
                                 "mercado_dist_max_52s"],
    "true regime": 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 == "no state" 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["truth 2012-2019"] = list(verdad.loc[False]) * 2
t.loc["truth 2020-2025"] = list(verdad.loc[True]) * 2
t.round(2)

                  shap                  pd
                 calma estres rebote calma estres rebote
no state          3.64   3.67   3.61  3.84   3.87   3.88
observable state  3.84   1.96   2.32  4.27   0.40   1.26
true regime       4.00   2.83  -0.94  4.58   2.00  -6.47
truth 2012-2019   5.01   1.74  -7.07  5.01   1.74  -7.07
truth 2020-2025   4.23   2.22 -12.68  4.23   2.22 -12.68

Without state, both tools give the same slope in all three regimes (3.6-3.9 bp): the model learned an average momentum and applies it always. The explanation is faithful to the model and the model is wrong in rebounds, where the truth is −12.68 in the test. With the true regime, partial dependence on rebound rows gives −6.47, very close to the −7.07 of training (the 2020-2025 rebounds are more violent than the ones the model saw). SHAP, on the other hand, gives −0.94: the momentum × regime interaction is split, and half goes to the rebote feature. Reading the SHAP of momentum as “momentum barely matters in rebounds” would be a mistake.

The observable state is not enough (+1.26 in rebound with partial dependence). The market's return over the last month and the last quarter, or the drawdown from highs, do not separate a rebound from stress in time, and with 26 rebound weeks in training the model has very few examples. It is the same conclusion as lesson 5.10.2 seen from the side of explanation: what the model does not see, it neither learns nor explains, and a clean explanation of a model without the right state gives a reassuring and false picture.

Stability of importances

An importance is an estimate and has an error. Three sources of variation, from least to most relevant:

• Seed: the same sample with different randomness in training (row and column subsampling).
• Bootstrap of weeks: another sample from the same process, resampling whole weeks with replacement (the rows of a week share that week's market; resampling individual rows would underestimate the variation).
• Time window: two different years of data. Besides sampling, the true effects change (the generator makes them vary with a weekly AR(1) and with the regime).

For each replica a somewhat lighter LightGBM is trained (150 trees, rate 0.1) and the mean absolute SHAP is measured on the same 1,000 test rows. Stability is the mean rank correlation between replicas.

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 [("alpha", a), ("label", y)]:
    for s in range(4):                 # same sample, another seed
        imp[(diana, "seed", s)] = importancia(
            lgbm(FEATS, obj, semilla=s, n=150, lr=0.1))
    for i0, i1 in ventanas:            # two years of data each
        filas = (d["fecha"] >= f"{i0}-01-01") & (
            d["fecha"] <= f"{i1}-12-31")
        imp[(diana, "window", i0)] = importancia(
            lgbm(FEATS, obj, filas.to_numpy(), n=150, lr=0.1))
    for b in range(4):                 # weeks with replacement
        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 ["seed", "window", "bootstrap"]}
              for diana in ["alpha", "label"]}).round(2)

           alpha  label
seed        0.99   0.81
window      0.90   0.31
bootstrap   0.98   0.51

With the noise-free target, the explanation is stable under everything: 0.99 across seeds, 0.98 across bootstraps and 0.90 across windows, the lowest because the true effects change from one period to another. With the real label it falls apart: 0.81 across seeds, 0.51 across bootstraps and 0.31 across windows. Stability across seeds is misleading, because the model reads the same noise of the same sample in a similar way; the uncertainty that matters is that of the data.

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] != "seed"]
res = {}
for diana in ["alpha", "label"]:
    rr = rangos[[c for c in replicas if c[0] == diana]]
    res[(diana, "min")] = rr.min(axis=1)
    res[(diana, "max")] = rr.max(axis=1)
pd.DataFrame(res).loc[sel].astype(int)

            alpha     label
              min max   min max
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

The figure at the top shows the rank ranges in the eight data replicas (windows and bootstrap). With the noise-free target, ret_1s is always first and ruido_1 is always among the last four. With the real label, valor_bp and sue, with a planted effect, range from 2nd to 18th place; ret_1m, with no effect of its own, from 1st to 17th; ruido_4 reaches 6th and ruido_1 8th. Only ret_1s stays at the top (from 1st to 5th), and sentimiento, which has an effect but 32% of its values missing and filled with zero, never gets above 14th. A report saying “the third most important feature is X” with these data would be reporting one realization of the noise.

When explanations mislead: correlated features

The first source of instability is the split within groups of correlated features: the model can use any of them and the split changes from one replica to another even if the total does not. This is checked with the correlation groups of lesson 5.2.5.

grupos = {"reversal": ["ret_1s", "ret_1m", "rsi_14"],
          "medium term": ["ret_3m", "dist_max_52s"],
          "momentum": ["mom_12_1", "ruido_4"],
          "volatility": ["vol_21d", "vol_63d"],
          "value-quality": ["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] == ("label", "bootstrap")]
cg = spearmanr(imp_g[boot]).statistic
print(f"label, bootstrap: mean rho by features "
      f"{rho_medio(boot):.2f} · by groups "
      f"{cg[np.triu_indices(4, 1)].mean():.2f}")
cuota = imp.loc[grupos["reversal"], boot]
(cuota / cuota.sum()).round(2)

label, bootstrap: mean rho by features 0.51 · by groups 0.83
           label
       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

Summing the SHAP of each group, the rank correlation across bootstraps goes from 0.51 to 0.83. Within the reversal group, the share of ret_1s ranges from 39% to 51% depending on the replica and ret_1m reaches 39% in one of them: the group is stable, the split is not.

The most extreme form is substitution: if a feature is removed, its substitute takes its place.

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({"with mom_12_1": shap_medio(m_a, FEATS),
                     "without mom_12_1": shap_medio(m_sin, sin_mom)})
p_sin = m_sin.predict(X.loc[te, sin_mom].to_numpy())
print("IC with the alpha in 2020-2025:",
      f"with {ic_semanal(m_a.predict(Xte), a, te):.3f} · without "
      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 with the alpha in 2020-2025: with 0.783 · without 0.777
corr(ruido_4, mom_12_1) = 0.6 · corr(dist_max_52s, mom_12_1) = 0.64
              with mom_12_1  without 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

Without mom_12_1, the mean SHAP of ruido_4 (0.6 × momentum plus noise, no effect of its own) goes from 0.03 to 0.94 bp and that of dist_max_52s (correlation 0.64), from 0.32 to 1.74. The model's IC barely changes (0.783 against 0.777): the momentum information is still inside, spread among its neighbors. For someone looking only at the explanation, the second model “proves” that the distance to the 52-week high is one of the most important features in the market.

The same traps, summarized:

How to present an explanation

flowchart TB
  M[Validated model] --> R[Replicas]
  R --> G[Groups]
  G --> E{Stable?}
  E -- no --> N[Do not interpret]
  E -- yes --> F[Shape: PD, ICE]
  F --> H[Hypothesis]

• First the model, then the explanation. Explaining a model that does not pass honest validation (lessons 5.4.2 and 5.4.4) is explaining noise.
• Always replicas. Bootstrap of time blocks and windows; report the range of each importance or of its rank, not a single number.
• Groups before individual features, with noise controls inside the model (lesson 5.2.5).
• SHAP for attribution, partial dependence and ICE for shape, and interaction values when the question is whether two features reinforce each other.
• An explanation is a hypothesis. If it suggests a mechanism (reversal is stronger with low volume), the mechanism is tested separately, with data that were not used to find it.

Exercises

1. By default TreeExplainer computes SHAP values by following the tree paths (tree_path_dependent), which partly respect the correlations of the training data. With feature_perturbation="interventional" and a background of 100 rows, it breaks those correlations. Compare the mean SHAP of the two versions on 500 rows. Where do they differ?

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({"interventional": np.abs(sv_int).mean(0),
                    "tree path": np.abs(sv_cam).mean(0)}, index=FEATS)
dif["difference"] = dif["interventional"] - dif["tree path"]
dif.sort_values("tree path", ascending=False).head(11).round(2)

              interventional  tree path  difference
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. Explain the same stock-week (the MACE row of 2020-05-29) with four LightGBMs with different seeds, trained with the noise-free target and with the real label. Do the three features with the highest SHAP agree?

fila = Xm[[0]]
for diana, obj in [("alpha", a), ("label", 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"  seed {s}: pred {m.predict(fila)[0]:6.2f} bp · "
              + ", ".join(f"{c} {sv1[c]:+.1f}" for c in top))

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

3. An investment committee asks every week for “the three reasons why the model is buying”. With what you have seen in the lesson, what would you deliver and with which caveats?