---
title: "5.10.1 SHAP, partial dependence and stability"
description: "SHAP values and interactions, partial dependence and ICE against planted effects, and how much importances change with seeds, windows and bootstrap."
url: "https://pythonparatrading.com/en/python-trading-course/level-5/shap-partial-dependence-stability/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-5/shap-partial-dependence-stability.png"
---

# 5.10.1 · SHAP, partial dependence and the stability of feature importances

> Explaining a model is easy to do and hard to read well. In the level's panel every effect is known, the interaction of reversal with volume and the sign change of momentum in rebounds, so you can check what each tool recovers, what it splits, what it cannot see and how much the explanation changes if the seed or the sample changes.

- Course: [Level 5: Machine learning and deep learning (Keras/PyTorch)](https://pythonparatrading.com/en/python-trading-course/level-5/)
- Module: 5.10 · Interpretability, robustness and model risk
- Duration: 75 min
- Prerequisites: [5.2.5 · Feature importance: MDI, MDA, SFI and SHAP](https://pythonparatrading.com/en/python-trading-course/level-5/feature-importance-mdi-mda-sfi-shap/), [5.5.3 · Gradient boosting: XGBoost, LightGBM and CatBoost](https://pythonparatrading.com/en/python-trading-course/level-5/xgboost-lightgbm-catboost-gradient-boosting/), [5.1.1 · Why financial machine learning is different](https://pythonparatrading.com/en/python-trading-course/level-5/why-financial-machine-learning-is-different/)

## By the end of this lesson

- Use the SHAP values of a tree to explain a specific prediction and the interaction values to separate the main effect of a feature from its interaction with another.
- Compute partial dependence and ICE curves with scikit-learn and read a planted interaction in the fan of curves.
- Distinguish what SHAP splits (the attribution) from what partial dependence measures (the conditioned effect), and recognize that a model cannot explain a state it does not see.
- Measure the stability of an importance with seeds, time windows and a bootstrap of weeks, and know when a ranking is noise.
- Recognize when explanations mislead with correlated features: substitution, arbitrary splitting within a group and extrapolation of partial dependence.

## 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.

```mermaid
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.

```python
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)
```
```output
(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.

```python
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}")
```
```output
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).

```python
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)
```
```output
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.

```python
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})")
```
```output
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$:

$$
\text{PD}_j(x) = \frac{1}{n}\sum_{i=1}^{n} \hat f\left(x,\, \mathbf{x}_{i,-j}\right)
$$

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).

```python
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)
```
```output
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
```

*Fig. — Centred ICE curves of ret_1s for the LightGBM trained on the noise-free target (20 stock-weeks from 2020-2025, grey) and their mean in each volumen_rel tercile, with the slope in bp per unit and the generator's slope.* Curves falling from left to right: the prediction drops as last week's return rises (reversal). The low-volume tercile mean, highlighted, is the steepest: −7.6 bp per unit against the generator's −8.4. Medium volume: −5.1 against −5.5. High volume: −2.5 against −2.4. The individual curves fan out: the slope depends on each row's volume.

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.

<aside class="callout" data-kind="aviso">

Partial dependence moves one feature while leaving the others fixed, even when they are correlated: with `ret_1s` at
+2 and `rsi_14` at −2 it evaluates the model on combinations that do not exist (a stock that rises sharply over the
week with an oversold RSI). There the model extrapolates and the curve can invent shapes. ALE curves (*accumulated
local effects*, Apley and Zhu, 2020) avoid the problem by accumulating local differences within each interval of the
feature, using only rows that are really there.

</aside>

## 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.

```python
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)
```
```output
                  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
```

*Fig. — Momentum slope in the 2020-2025 rebound weeks according to partial dependence (PD) and SHAP, for three LightGBMs trained on the noise-free target: with the true regime as a feature, with an observable market state and without state.* Horizontal bars. Training truth: −7.07. PD of the model with the true regime: −6.47, close to the truth. SHAP of the same model, highlighted: −0.94, because it splits the interaction with the regime feature. PD with an observable state: +1.26. PD without state: +3.88, the same effect as in calm.

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.

```python
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)
```
```output
           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.

```python
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)
```
```output
            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.

```python
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)
```
```output
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.

```python
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)
```
```output
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:

| Situation | What the explanation shows | What to do |
|---|---|---|
| Correlated features | arbitrary split among them, changes with the replica | importance by groups |
| Feature removed or added | its substitute inherits the importance | read each importance with its list of substitutes |
| Interaction | SHAP splits it between the two features | interaction values, conditioned ICE |
| Unobserved state | a mean effect faithful to the model and false in the regime | add the state or do not trust the effect in that regime |
| Partial dependence with correlations | impossible combinations, extrapolation | ALE, or PD only within the joint range of the data |
| Noisy label | rankings that change with the sample | replicas, rank ranges, noise controls |

## How to present an explanation

```mermaid
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?

**Solution.**

```python
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)
```
```output
              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
```

The order barely changes and the differences are small for the features with an effect (`ret_1s` 4.63 against 4.37,
`mom_12_1` 3.22 against 3.06). The largest relative differences are in features with no effect of their own that are
correlated with others: `dist_max_52s` rises from 0.33 to 0.46 (39%) and `ret_1m` from 0.25 to 0.33. The
interventional version replaces the “absent” features with values from background rows chosen at random, breaking
their correlation with the present ones: it evaluates the model on unrealistic combinations and attributes somewhat
more to the correlated neighbors. It answers “what does the model do if I change this feature and nothing else”; the
tree-path version is closer to “what does it do with the data as they come”. With nearly independent features the two
agree; with highly correlated groups it is worth looking at both or, better, grouping.

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?

**Solution.**

```python
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))
```
```output
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
```

With the noise-free target, the four seeds tell the same story: a prediction between 14.3 and 15.5 bp and the same
three reasons (`valor_bp` +7.1 to +7.3, `vol_63d`, `ey`). With the real label, the same row receives predictions of
24.2, 2.9, −0.2 and −0.1 bp, and the reasons change with the seed: `valor_bp` appears three times, but alongside
`ret_1s`, `dist_max_52s`, `rsi_14`, `dias_desde_informe` and even `ruido_5`, pure noise, with +5.5 bp. Each explanation
is exact for its model; the problem is that models trained with such a noisy label differ a lot from one another on
each row even though their mean IC is similar. A local explanation of a model with little signal has to be averaged
over several models (or presented as a range) before telling it to anyone.

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?

**Solution.**

Per row, the model's SHAP values aggregated by groups of correlated features (value-quality, reversal, momentum,
volatility...), not by individual features, because the split within a group is arbitrary and changes with the seed and
the sample. The three reasons would be the stock's three groups with the highest SHAP, with their sign and their size in
basis points, and the total prediction, which is their sum plus the base. The caveats: the explanation describes the
model, not the market; with a signal-to-noise ratio like that of equities, the model mainly explains how it splits a
small effect, and the explanations of a row with the real label change a lot from one seed to another (exercise 2), so
it is advisable to average several models or show the range. If a group dominates disproportionately, it has to be
audited before being presented (it may be a leak, lesson 5.2.4). And periodically, global stability: importances by
groups with a bootstrap of weeks and by windows, and the shape of the main effects with partial dependence and ICE, to
detect whether the model has changed its criteria (lesson 5.10.2).

## Self-assessment

1. For one row, the explainer's base value is 0.01 bp, the sum of its SHAP values is 15.64 bp and the model's prediction is 15.65 bp. Which SHAP property does this show?

   - A. That the model is overfitted
   - B. Additivity (efficiency): the base value plus the sum of the row's SHAP values is exactly the prediction
   - C. That SHAP measures the causal importance of each feature
   - D. That the row is an outlier

   *Answer B: Additivity (efficiency): the base value plus the sum of the row's SHAP values is exactly the prediction.* Shapley values split the difference between the prediction and the base value without leaving anything out. In the cell the maximum difference between the two is 1e-13: rounding error. Additivity says how the prediction is split, not whether the features cause anything.

2. The planted interaction between `ret_1s` and `volumen_rel` is worth −0.5·c_rev per unit of the product (2.73 bp with the mean c_rev of 2012-2019). What slope is expected in each of the pair's two SHAP interaction cells, (`ret_1s`, `volumen_rel`) and (`volumen_rel`, `ret_1s`)?

   - A. 2.73 in each
   - B. Zero, because SHAP does not measure interactions
   - C. 5.46, twice as much
   - D. Half, 1.37: the interaction matrix is symmetric and splits the cross term equally between the two cells

   *Answer D: Half, 1.37: the interaction matrix is symmetric and splits the cross term equally between the two cells.* The interaction values of Lundberg et al. (2020) put half of the joint effect in each cell of the pair. The LightGBM gives 0.93 against 1.37: it learned the interaction at about two thirds of its size, as expected from a regularized model.

3. With the true regime as a feature, momentum in a rebound has a partial dependence slope of −6.47 and a SHAP slope of −0.94 (truth −7.07). What is the correct reading?

   - A. Both are correct and answer different questions: partial dependence on rebound rows measures the effect of momentum given the regime; SHAP splits the momentum × regime interaction between the two features, and half goes to the regime feature
   - B. SHAP is computed wrongly
   - C. Partial dependence is biased by the correlation
   - D. The model did not learn the sign change

   *Answer A: Both are correct and answer different questions: partial dependence on rebound rows measures the effect of momentum given the regime; SHAP splits the momentum × regime interaction between the two features, and half goes to the regime feature.* For a term b·x·z, the SHAP of x carries roughly half of b·x·z (plus a term with the mean of z). If you want to know how the prediction responds to momentum in a regime, the tool is partial dependence or ICE conditioned on that regime.

4. With the real label, the mean rank correlation between the SHAP importances of four seeds is 0.81, that of four bootstraps of weeks 0.51, and that of four two-year windows 0.31. Which one measures the uncertainty that matters?

   - A. The seeds one, because it is the highest
   - B. None: SHAP is exact for trees
   - C. The bootstrap and the windows: seeds only change the randomness of training on the same data, and the doubt lies in the data (another sample, another period)
   - D. Only the windows

   *Answer C: The bootstrap and the windows: seeds only change the randomness of training on the same data, and the doubt lies in the data (another sample, another period).* An importance that is stable across seeds can be a stable reading of the same noise. The bootstrap of weeks imitates another sample from the same process; the windows add that the effects change over time. With the noise-free target all three are above 0.90.

5. When `mom_12_1` is removed, the mean SHAP of `ruido_4` (no effect of its own, correlation 0.60 with momentum) goes from 0.03 to 0.94 bp and the model's IC barely changes (0.783 against 0.777). What does this teach?

   - A. That `ruido_4` has an effect that momentum was hiding
   - B. That a high explanation does not prove that the feature causes anything: the model uses whatever correlated feature it has at hand, and importance moves between substitutes depending on which features are included
   - C. That momentum should always be removed
   - D. That SHAP fails with continuous features

   *Answer B: That a high explanation does not prove that the feature causes anything: the model uses whatever correlated feature it has at hand, and importance moves between substitutes depending on which features are included.* It is the substitution effect, now in SHAP. `dist_max_52s` (correlation 0.64) rises from 0.32 to 1.74. The importance of a feature can only be read together with the list of its substitutes; that is why features are grouped.

## Key points

- SHAP is additive: the base value plus the sum of a row's values is the prediction. That makes it possible to explain a specific decision, but it explains the model, not the world.
- SHAP interaction values find the planted reversal × volume interaction as the strongest pair (1.47 bp) and measure its size: 0.93 per cell against 1.37 for the truth. The following pairs mix a planted interaction (SUE × days since the report) with others that are not planted.
- Partial dependence gives the mean effect and the ICE curves its heterogeneity. Grouped by volume terciles, the ICE slopes recover the full interaction (1.48 and 0.48 against 1.52 and 0.43 for the truth), which SHAP splits.
- A model can only explain what it sees. Without the regime, momentum has the same slope in calm and in rebound (3.9 bp); with the true regime, partial dependence in rebound gives −6.5 against −7.1; with an observable market state, +1.3: rebounds are not recognized in time.
- With the real label the importances are unstable: a rank correlation of 0.51 across bootstraps of weeks and of 0.31 across windows, and `valor_bp` ranges from 2nd to 18th place. By groups of correlated features it rises to 0.83.
- With correlated features importance moves between substitutes: without momentum, `ruido_4` and `dist_max_52s` take its place and the model predicts almost the same. Group, measure replicas and report the range, not the rank.

## Resources

### Documentation

- [Lundberg and Lee (2017): A Unified Approach to Interpreting Model Predictions, arXiv:1705.07874](https://arxiv.org/abs/1705.07874) · arxiv.org
- [Lundberg et al. (2020): From local explanations to global understanding with explainable AI for trees, Nature Machine Intelligence 2, 56-67](https://doi.org/10.1038/s42256-019-0138-9) · doi.org
- [Friedman (2001): Greedy function approximation: a gradient boosting machine, Annals of Statistics 29(5) (partial dependence)](https://doi.org/10.1214/aos/1013203451) · doi.org
- [Goldstein, Kapelner, Bleich and Pitkin (2015): Peeking Inside the Black Box (ICE curves), JCGS 24(1)](https://doi.org/10.1080/10618600.2014.907095) · doi.org
- [Apley and Zhu (2020): Visualizing the Effects of Predictor Variables in Black Box Supervised Learning Models (ALE), JRSS B 82(4)](https://doi.org/10.1111/rssb.12377) · doi.org
- [scikit-learn: partial dependence and individual conditional expectation](https://scikit-learn.org/stable/modules/partial_dependence.html) · scikit-learn.org

### Data

- [Level 5 data: true coefficients and regime of each week (CSV)](https://pythonparatrading.com/curso/datos/nivel-5/panel5_coef_verdad.csv) · Download · panel5_coef_verdad.csv

Previous lesson: [5.9.3 · Why reinforcement learning overfits in directional trading](https://pythonparatrading.com/en/python-trading-course/level-5/why-rl-overfits-directional-trading/)

Next lesson: [5.10.2 · Data drift and concept drift: detecting them and what to do](https://pythonparatrading.com/en/python-trading-course/level-5/data-drift-concept-drift/)