Feature selection

Every other module here builds one more feature. This one answers the question that follows: you have 200+ indicators — which of them are actually different from each other?

This module is not an expression API

selection is the deliberate exception to the expression convention. Choosing features is a cross-feature question that needs the whole materialized matrix at once, so these functions take a pl.DataFrame and return NumPy arrays and plain Python objects. Nothing here is on the .ta namespace, and nothing here is lazy or streaming-safe.

Why you need this

polars_ta exposes over 200 indicators. They are nowhere near 200 independent measurements: rsi(14), rsi(21), cmo, stochrsi and roc are one underlying quantity read five ways. Handing all of them to a model is not just wasteful — collinearity actively destroys feature-importance scores. Two near-identical features split the credit between them and both look useless, while a mediocre-but-unique third feature outranks them. This is the substitution effect, and it breaks MDI and MDA alike.

So the order matters: remove redundancy first, score relevance second. This module is the first half. It never looks at a target, which is precisely why it is safe to run first — with no label involved, there is no label to leak.

Which features do I keep?

The short answer, and the one thing to take from this page:

result = selection.select_features(feats, FEATURES, k=6)
report = result.report()

report() gives one row per feature with a decision column. That column is the answer.

decision	Meaning	What to do
keep	This feature is its cluster's representative	Model on it
drop	Its representative correlates ≥ 0.7 with it	Discard — nothing is lost
review	Its representative correlates < 0.7 with it	Do not drop blindly — the partition is too coarse

Real output, 16 indicators at k=4:

feature	cluster	kept	representative	similarity	decision
rsi_21	0	✅	rsi_21	1.00	keep
rsi_14	0		rsi_21	0.80	drop
mfi	0		rsi_21	0.77	drop
cmo_14	0		rsi_21	0.75	drop
roc_10	0		rsi_21	0.66	review
obv	0		rsi_21	0.62	review
stochrsi	1	✅	stochrsi	1.00	keep
adx_14	2	✅	adx_14	1.00	keep
hurst	2		adx_14	0.76	drop
vol_21	3	✅	vol_21	1.00	keep
ulcer	3		vol_21	0.69	review

How to choose k: raise it until nothing says review

This is the concrete rule the rest of the page builds up to.

for k in (4, 6, 8):
    n = (select_features(feats, FEATURES, k=k).report()["decision"] == "review").sum()
    print(k, n)
# 4 → 4 flagged
# 6 → 0 flagged   ← use this
# 8 → 0 flagged

k=6 is the smallest cut where every drop is safe — and it agrees with effective_rank (5.66), which is a reassuring independent check. Note that signal_rank said 4, which was too aggressive: three features would have been thrown away that nothing left behind could reconstruct.

Two numbers, one decision

effective_rank proposes k; the review count verifies it. If they disagree, trust the review count — it is measured against the features you actually have, not against a random-matrix assumption your data may violate.

The pipeline

Stage	Function	What it does
0	nonstationary_features	Flag integrated features before you correlate anything
0	screen	Materialize the matrix; drop sparse, warm-up and constant columns
1	corr_matrix	Spearman (default) or Pearson correlation
1	vif	Multicollinearity a pairwise screen cannot see
1	mutual_information	Non-linear dependence; works on the sparse candle patterns
1	variation_of_information	The same, as a true metric distance for clustering
2	denoise	Marchenko-Pastur eigenvalue clipping, optional detoning
3	cluster_features	Hierarchical clustering, k chosen by silhouette
3	cluster_by_distance	The same engine on any distance — e.g. variation of information
4	select_representatives	One feature per cluster
1-4	select_features	Runs screen through select_representatives and returns a SelectionResult

select_features deliberately does not run nonstationary_features for you — fixing a flagged feature means choosing a transform, which is your call, not a default. Run it first and act on what it reports.

Diagnostics: signal_rank and effective_rank tell you how many dimensions your feature set really has.

Stage 0 — stationarity is not optional

Correlations between integrated series are spurious (Granger & Newbold, 1974): two independent random walks correlate strongly by construction. Denoise a matrix built on raw levels and you are carefully denoising nonsense.

Several indicators in this library are integrated by design — on_balance_volume, adi, volume_price_trend, cumulative_return, and every price-level moving average. nonstationary_features flags them with a lag-1 autocorrelation screen (a heuristic, not a unit-root test). Fix a flagged feature with quant.frac_diff — which keeps most of the memory — or quant.rolling_z_score, then re-run.

Stage 2 — what Marchenko-Pastur actually tells you

For a correlation matrix built from \(T\) observations of \(N\) independent variables with variance \(\sigma^2\), the eigenvalues concentrate on

\[\lambda_\pm = \sigma^2\left(1 \pm \sqrt{1/q}\right)^2, \qquad q = T/N.\]

Any empirical eigenvalue below \(\lambda_+\) is therefore consistent with pure noise. With 205 features over 5000 bars, \(q \approx 24\) and \(\lambda_+ \approx 1.45\) — and the average eigenvalue of any correlation matrix is exactly 1, so that threshold cuts deep.

denoise replaces the entire noise bulk with its average. This preserves the trace while destroying the spurious structure, and it lifts the smallest eigenvalues away from zero — the practical payoff, since it is the near-singular directions that wreck anything that inverts the matrix.

\(\sigma^2\) is estimated by fitting the theoretical MP density to a kernel density estimate of the observed eigenvalues (López de Prado, 2020, ch. 2). The tempting shortcut — iterate "everything below the edge is noise, its mean is the next \(\sigma^2\)" — is a trap: with real signal present it spirals downward until most of the bulk is misclassified as signal.

Detoning (detone=1) additionally strips the largest eigenvector. For technical indicators on a single asset the dominant mode is essentially "trend/level" and it swamps everything else, so detoning is what you want when you are after features that say something beyond direction. It is off by default because it is a modelling choice, not a correction.

Stage 3 — clustering, and why this distance

\[d_{ij} = \sqrt{\tfrac{1}{2}\left(1 - \rho_{ij}\right)}\]

Unlike \(1 - \rho\), this is a true metric (Mantegna, 1999) — it satisfies the triangle inequality, which is what makes hierarchical clustering on it meaningful rather than merely suggestive. Note that perfectly anti-correlated features land maximally far apart even though they carry the same information with a flipped sign; pass np.abs(corr) if that matters for your feature set.

The number of clusters is chosen by maximising the mean silhouette (Rousseeuw, 1987) over every cut of the dendrogram — the correlation-clustering step of López de Prado's ONC, minus the recursive re-clustering of low-quality clusters.

Part 2 — is it real, and does it predict?

Everything above is target-free and therefore leak-free. These two are not, and they are where the answers get hard.

Significance: block_permute / permutation_test

With 200 candidate features, spurious winners are guaranteed. The defence is a null distribution — but not from a plain shuffle, which would also destroy the autocorrelation every rolling indicator has by construction and hand you a null far too easy to beat. block_permute rebuilds each column from randomly placed contiguous blocks of itself, so its serial dependence survives while its alignment with every other column is destroyed.

This test found a real problem with signal_rank

Marchenko-Pastur assumes i.i.d. rows. Rolling indicators are ~0.99 autocorrelated, so n_obs grossly overstates the effective sample size and the edge sits too low. On the 16-indicator case study:

Statistic	Observed	Null	p
Eigenvalues above the edge	4	5.1 ± 1.2	0.99 — not significant
Top eigenvalue	6.07	1.54 ± 0.08	0.008
effective_rank	5.66	15.5	0.008

The shared structure is unambiguously real; the precise count at the edge is not. Prefer effective_rank as a headline, and treat signal_rank as a ceiling unless you have tested it.

Relevance: the feature-importance toolkit

Everything above is target-free. This is the half that brings in a label — and therefore the half where leakage, the substitution effect and multiple testing all live at once.

Start with the cross-validation, because every measure below sits on top of it. Purged K-fold with embargo (López de Prado, 2018, ch. 7) is not an optimization, it is a correctness fix: bars \(t\) and \(t+1\) share \(w-1\) observations of any \(w\)-window indicator, and a forward-return label spans \(h\) more bars. A plain K-fold therefore trains on data that overlaps its own test fold and reports memorization as skill. purged_kfold drops training samples whose label window reaches into the test fold, then embargoes a further \(h + \max(w)\) bars to handle the serial correlation in the features themselves.

Function	Family	Looks at	Cost	Answers
target_screen	Filter	One feature at a time	Trivial	Does this co-move with the label at all?
single_feature_importance	Wrapper (SFI)	One cluster, alone	\(k\) fits/fold	What is this worth by itself?
clustered_mdi	Embedded (MDI)	All features, in-sample	Free with the fit	What did the model use?
clustered_mda	Wrapper (MDA)	All features, out-of-sample	\(k\) scores/fold	What breaks if I remove this last?
clustered_shapley	Wrapper	Every coalition	\(2^k\) fits/fold	What is this worth on average, over all orderings?
null_importance	Calibration	Any of the above	\(n_{\text{draws}} \times\)	Is any of it beyond chance?

The model is always a callable you supply (fit/predict), so scikit-learn and LightGBM stay your dependency, not this library's; RidgeRegressor is the zero-dependency default. clustered_mdi is the one exception — it needs a model that reports feature_importances_, so it has no default at all.

Step 1 — target_screen: the free triage

Before fitting anything, ask whether a feature moves with the label. Two statistics, both computed out of sample on the purged test folds:

• The information coefficient, a Spearman rank correlation (Grinold & Kahn, 1999). Rank rather than Pearson because one outlier bar otherwise sets the number. On real forward returns a mean \(|IC|\) of 0.02-0.05 is a useful feature; above 0.1 you should suspect your label before celebrating.
• Mutual information, which catches what an IC structurally cannot. A U-shaped relationship — both volatility tails predict, the middle does not — scores an IC of zero and an MI well above it. High MI with near-zero IC means non-monotone, not useless: it needs a model that can represent a fold.

ic_t_stat measures consistency, not size

A negligible-but-reliably-signed IC gets a large t-stat. In this library's own tests a pure-noise feature scores \(|IC|\) 0.008 with a t-stat near 5 — a stable nothing. Sort by abs_ic; use the t-stat only to discard what near the top fails to reproduce across folds.

Step 2 — the three classical measures, and why you run all three

MDI, MDA and SFI (López de Prado, 2018 ch. 8; 2020 ch. 6) are not competing implementations of one idea. They answer three different questions, and their disagreement is the diagnostic:

• MDI — in-sample impurity reduction, summed within each cluster. Free with a tree fit and the only measure that sees interactions directly, since a deep split is conditional on everything above it. But it is in-sample, it is biased toward high-cardinality features (a continuous feature offers a tree far more split points than a 0/±100 candlestick pattern), and it has no null — a pure-noise feature still scores above zero.
• MDA — the drop in out-of-sample score when a cluster is shuffled. The strongest of the three, and the only one with a natural zero.
• SFI — fit each cluster alone. Structurally immune to substitution because nothing is left to substitute with, and equally blind to interaction.

SFI	MDA	Reading
high	high	Genuinely predictive. Keep.
high	low	Redundant — something else already carries it.
low	high	Only works in combination. Keep, but never alone.
low	low	Nothing there.

Similarly, high MDI with zero MDA means the model memorized that feature — the single most useful cross-check the pair provides.

Step 3 — clustered_shapley: the one that adds up

MDA asks what breaks if a cluster is removed last; SFI asks what it earns first. On interacting features those differ, with no principled reconciliation — which is exactly what the Shapley value provides, by averaging a cluster's marginal contribution over every order in which clusters could be added:

\[\phi_i = \sum_{S \subseteq N \setminus \{i\}} \frac{|S|!\,(|N|-|S|-1)!}{|N|!}\,\bigl[v(S \cup \{i\}) - v(S)\bigr]\]

SFI is the \(S = \emptyset\) term, MDA is close to the \(S = N \setminus \{i\}\) term, and \(\phi_i\) is the properly weighted average of those and everything between. The payoff is efficiency: the values sum exactly to the full model's score over an empty one, so "cluster 0 is 40% of the edge" becomes a defensible statement rather than a vibe. It is also the only measure here that reports a negative importance honestly.

Cost is \(2^k\) fits per fold

32 fits for \(k=5\), 4096 for \(k=12\), times n_splits. This is why the function takes clusters rather than features — per-feature Shapley on a 200-column matrix is not slow, it is infeasible. max_clusters (default 12) exists to stop you finding that out the hard way.

Step 4 — null_importance: the number that decides

An importance of 0.03 is not a finding until you know what 0.03 looks like when nothing is there — and it is rarely zero. With 200 candidates, a flexible model and a weak label, the best cluster always shows a positive score.

The fix (Altmann et al., 2010) breaks the one link that matters — feature to label — and leaves everything else intact. Features keep their correlation structure, the model keeps its capacity, the CV keeps its geometry; only the target is scrambled, and it is scrambled with block_permute rather than a plain shuffle, because a forward return is autocorrelated and destroying that too would make the null far too easy to beat.

null = selection.null_importance(
    scored, sel.names, "fwd", labels=sel.labels,
    n_draws=100, label_horizon=12, embargo=300,
)
null.report()                  # observed vs chance, with z-scores
null.significant(alpha=0.05)   # the ids surviving Bonferroni

significant() applies a Bonferroni correction by default, because testing 20 clusters at \(\alpha = 0.05\) produces a false positive about 64% of the time. Note that n_draws floors the smallest resolvable p-value at \(1/(n_\text{draws}+1)\) — 50 draws cannot resolve below 0.02, which is not enough to clear Bonferroni past a handful of clusters.

Worked example: the whole library against a real label

examples/plot_feature_importance.py runs all six on every indicator polars_ta exposes — 205 built automatically, screened to 183 — clustered into 8 groups and ranked against a 12-bar forward return on 5000 BTCUSDT 5m bars, with label_horizon=12 and embargo=300.

The univariate screen looks encouraging at first: the top features reach \(|IC| \approx 0.047\) with t-stats near 3.8, which on a research desk would pass for a signal. Then the model-based measures arrive — MDI from a RandomForestRegressor (120 trees, depth 4), the rest from the default ridge:

Cluster	n	Example member	SFI	MDA	MDI	Shapley
0	44	acc_dist_index	0.0095	−0.0022	0.191	0.0171
7	5	cdl_harami	0.0203	0.0010	0.000	0.0100
6	14	cdl_3_white_soldiers	0.0190	0.0077	0.001	0.0097
5	2	cdl_3_line_strike	0.0073	0.0012	0.000	0.0013
4	2	cdl_2_crows	−0.0095	−0.0009	0.003	−0.0037
3	7	approximate_entropy	−0.0116	−0.0047	0.074	−0.0049
2	63	adx_pos	−0.0399	−0.0141	0.262	−0.0166
1	46	adx	−0.0409	0.0001	0.469	−0.0208

The full-model out-of-sample IC is −0.0079. Negative. The Shapley values sum to exactly that, as efficiency requires.

MDI ranks the ranking exactly backwards

This is the memorization signature, on real data. The forest spends 47% of its total impurity reduction on cluster 1 — more than on all seven other clusters' worth of continuous features combined — and that cluster's out-of-sample MDA is 0.0001, its Shapley value −0.021. Clusters 1 and 2 together take 73% of the MDI budget while contributing nothing out of sample.

Meanwhile the three candlestick clusters (5, 6, 7) score essentially zero MDI — 0.001 between them — despite cluster 6 having the best MDA and the lowest p-value of any cluster in the study. That is MDI's cardinality bias stated in numbers: a continuous feature offers a tree hundreds of split points, a 0/±100 pattern offers two, and MDI rewards the former for the opportunity alone.

Read MDI alone and you would keep the 109 trend/momentum features and throw away the only cluster that came close to significance.

Not one cluster survives the null

With 50 shuffled-label draws, the best p-value is 0.118 (cluster 6) and nothing clears Bonferroni at 5% — or even an uncorrected 5%. The encouraging \(|IC| = 0.047\) from the univariate screen does not survive contact with a model, purged cross-validation, and a null.

This is the correct result, not a failed experiment. 200 technical indicators on 5000 bars of one asset, against a 12-bar forward return, do not predict it. A pipeline that reported otherwise would be measuring its own leakage. Every guard in this module — purging, embargoing, block permutation, Bonferroni — exists to make this the outcome rather than a confident, wrong number.

Two secondary readings worth taking from the same table:

• SFI and MDA disagree sharply on clusters 1 and 2 (SFI ≈ −0.04, MDA ≈ 0). Read via the table above, that is the low-SFI / low-MDA quadrant with noise on top: these 109 trend and momentum features are worthless alone and cost nothing when removed.
• Cluster 0 has the top Shapley value but a negative MDA. Its contribution is real only in combination, and only in some orderings — the precise case MDA alone would have scored as harmful and SFI alone as mediocre. It is also well within the null, so the honest conclusion is still "nothing here".

---

polars_ta.selection

Feature selection: find the handful of independent signals hiding in a wide indicator matrix.

Every other module in this library produces a pl.Expr — one column, computed lazily. This module is the exception, and deliberately so: choosing features is a cross-feature question that needs the whole materialized matrix at once, so these functions take a pl.DataFrame and return NumPy arrays and plain Python containers. Nothing here is registered on the .ta namespace.

The problem it solves: polars_ta exposes 200+ indicators, and they are nowhere near 200 independent measurements. rsi(14), rsi(21), cmo, stochrsi and roc are one underlying quantity read five ways. Feeding all of them to a model is not just wasteful — collinearity actively destroys every per-feature importance score, because two near-identical features split the credit between them and both look useless (the "substitution effect").

So the work splits in two, and the order is not negotiable.

Part 1 — redundancy (no target, therefore no possible leakage):

1. screen — materialize the matrix, drop sparse, warm-up and degenerate columns.
2. corr_matrix / mutual_information / vif — linear, rank, and non-linear dependence, plus multicollinearity that pairwise screening misses.
3. denoise — Marchenko-Pastur eigenvalue clipping: everything below the random-matrix edge is indistinguishable from noise, so it is flattened rather than trusted. Optionally detone the dominant market mode.
4. cluster_features / cluster_by_distance — hierarchical clustering on the correlation distance or on variation of information.
5. select_representatives — one feature per cluster.

select_features runs 1-5 and hands back a SelectionResult. Run nonstationary_features yourself first — correlations between integrated series are spurious, and select_features does not apply a fix for you because choosing the transform is a modelling decision, not a default.

Part 2 — is it real, and does it predict?

• permutation_test with block_permute — a null distribution for any statistic of the matrix. With 200 candidate features, spurious winners are guaranteed; this is what distinguishes structure from chance.
• target_screen — the cheap univariate triage: information coefficient and mutual information against the label, per feature, out of sample. No model, so no model to blame; run it before anything expensive.
• clustered_mda with purged_kfold — importance against an actual target. This does look at a label, so it is where leakage lives: the cross-validation purges and embargoes overlapping samples, and shuffling is done per cluster so the substitution effect cannot scramble the scores.
• clustered_mdi — the in-sample counterpart, for tree models that report feature_importances_. Fast and free with the fit, but in-sample and biased toward high-cardinality features, so it is a cross-check on MDA, not a replacement.
• single_feature_importance (SFI) — fit each feature alone. Structurally immune to substitution (nothing is left to substitute with) and equally blind to interaction, which makes the MDA/MDI/SFI disagreement itself informative.
• clustered_shapley — exact Shapley attribution of the out-of-sample score across clusters. Unlike MDA it is additive: the values sum to what the full model earns over an empty one, so "cluster A is worth 40% of the edge" is a statement you can actually make.
• null_importance — refit against a shuffled label many times to get a per-cluster null. This turns any of the scores above into a p-value, which is the only defensible way to draw the keep/drop line on 200 candidates.

The three importance families answer different questions and are meant to be read together (López de Prado, 2020, ch. 6): MDI says what the model used in-sample, MDA what it needs out of sample, SFI what each feature is worth alone. A feature that scores high on MDI and zero on MDA was memorized; high on SFI and zero on MDA is redundant with something else; high on MDA and low on SFI only works in combination.

References

• Marchenko & Pastur (1967), Distribution of eigenvalues for some sets of random matrices.
• López de Prado (2018), Advances in Financial Machine Learning, ch. 7 (purged K-fold cross-validation with embargo) and ch. 8 (MDI, MDA and SFI feature importance).
• López de Prado (2020), Machine Learning for Asset Managers, ch. 2 (denoising and detoning), ch. 3 (information-theoretic distance), ch. 4 (optimal clustering) and ch. 6 (clustered MDI and MDA).
• Mantegna (1999), Hierarchical structure in financial markets — the correlation-to-distance metric.
• Meila (2007), Comparing clusterings — variation of information as a metric.
• Rousseeuw (1987), Silhouettes — choosing the number of clusters.
• Hacine-Gharbi et al. (2012) — optimal histogram binning for mutual information.
• Altmann et al. (2010), Permutation importance: a corrected feature importance measure — the shuffled-label null behind null_importance.
• Shapley (1953), A value for n-person games; Štrumbelj & Kononenko (2014) for its use as a feature attribution.
• Grinold & Kahn (1999), Active Portfolio Management — the information coefficient as the unit of predictive skill.

ScreenResult dataclass

ScreenResult(values: ndarray, names: list[str], dropped: list[str] = list(), n_rows_dropped: int = 0, sparse: list[str] = list(), missing_rate: dict[str, float] = dict())

Output of :func:screen — a clean feature matrix and what was removed.

q property

q: float

Marchenko-Pastur aspect ratio T / N (observations per feature).

worst_missing

worst_missing(limit: int = 5) -> list[tuple[str, float]]

The limit columns with the highest non-finite rate, worst first.

Reach for this when n_rows_dropped is large: row removal is driven by the worst column, so a single sparse feature can cost the whole matrix most of its sample.

Source code in polars_ta/selection.py

def worst_missing(self, limit: int = 5) -> list[tuple[str, float]]:
    """The `limit` columns with the highest non-finite rate, worst first.

    Reach for this when `n_rows_dropped` is large: row removal is driven by
    the *worst* column, so a single sparse feature can cost the whole
    matrix most of its sample.
    """
    ranked = sorted(self.missing_rate.items(), key=lambda kv: -kv[1])
    return [(name, rate) for name, rate in ranked[:limit] if rate > 0.0]

SelectionResult dataclass

SelectionResult(selected: list[str], clusters: dict[int, list[str]], labels: ndarray, names: list[str], dropped: list[str], corr: ndarray, silhouette: float, signal_rank: int, effective_rank: float, n_obs: int)

Output of :func:select_features.

report

report(min_similarity: float = 0.7) -> DataFrame

Per-feature keep/drop table, with the evidence for each decision.

selected tells you what survived; this tells you why, and — more usefully — which of the drops you should not trust.

A cluster representative stands in for its members only as well as it correlates with them. Dropping a feature that correlates 0.97 with the kept one costs you nothing; dropping one that correlates 0.35 throws away real information, and means the partition was too coarse. Those two cases are indistinguishable in selected alone, which is why this exists.

Columns:

Column	Meaning
feature	The feature name
cluster	Which cluster it landed in
cluster_size	How many features share that cluster
kept	Whether it is the cluster's representative
representative	The feature kept in its place
similarity	abs correlation to that representative, in [0, 1]
decision	"keep", "drop", or "review"

The decision rule:

• keep — this is the representative; model on it.
• drop — similarity >= min_similarity, so the representative carries what this feature carried. Safe to discard.
• review — similarity < min_similarity. The representative is a poor stand-in. Either raise k so this feature gets its own cluster, or keep it alongside the representative.

A high review count means the clustering was cut too coarse, and is the signal to increase k — a more grounded way to choose it than silhouette alone.

report = result.report()
report.filter(pl.col("decision") == "review")   # the risky drops
report.filter(pl.col("kept"))                   # the final feature set

Note that similarity uses absolute correlation: a representative correlating -0.95 with a member carries that member's information with a flipped sign, which a linear model recovers for free.

Parameters:

Name	Type	Description	Default
min_similarity	float	Correlation below which a drop is flagged for review rather than treated as safe.	0.7

Returns:

Type	Description
DataFrame	One row per surviving feature, sorted by cluster then by descending
DataFrame	similarity, so each cluster reads as "the kept one, then what it
DataFrame	replaced, worst last".

Source code in polars_ta/selection.py

def report(self, min_similarity: float = 0.7) -> pl.DataFrame:
    """Per-feature keep/drop table, with the evidence for each decision.

    `selected` tells you *what* survived; this tells you **why**, and — more
    usefully — which of the drops you should not trust.

    A cluster representative stands in for its members only as well as it
    correlates with them. Dropping a feature that correlates 0.97 with the
    kept one costs you nothing; dropping one that correlates 0.35 throws
    away real information, and means the partition was too coarse. Those
    two cases are indistinguishable in `selected` alone, which is why this
    exists.

    Columns:

    | Column | Meaning |
    |---|---|
    | `feature` | The feature name |
    | `cluster` | Which cluster it landed in |
    | `cluster_size` | How many features share that cluster |
    | `kept` | Whether it is the cluster's representative |
    | `representative` | The feature kept in its place |
    | `similarity` | `abs` correlation to that representative, in `[0, 1]` |
    | `decision` | `"keep"`, `"drop"`, or `"review"` |

    The **decision rule**:

    - `keep` — this is the representative; model on it.
    - `drop` — `similarity >= min_similarity`, so the representative
      carries what this feature carried. Safe to discard.
    - `review` — `similarity < min_similarity`. The representative is a
      poor stand-in. Either raise `k` so this feature gets its own cluster,
      or keep it alongside the representative.

    A high `review` count means the clustering was cut too coarse, and is
    the signal to increase `k` — a more grounded way to choose it than
    silhouette alone.

    ```python
    report = result.report()
    report.filter(pl.col("decision") == "review")   # the risky drops
    report.filter(pl.col("kept"))                   # the final feature set
    ```

    Note that `similarity` uses **absolute** correlation: a representative
    correlating -0.95 with a member carries that member's information with
    a flipped sign, which a linear model recovers for free.

    Args:
        min_similarity: Correlation below which a drop is flagged for
            review rather than treated as safe.

    Returns:
        One row per surviving feature, sorted by cluster then by descending
        similarity, so each cluster reads as "the kept one, then what it
        replaced, worst last".
    """
    index = {name: position for position, name in enumerate(self.names)}
    representative_of = {
        cluster: next(name for name in self.selected if name in members)
        for cluster, members in self.clusters.items()
    }

    rows = []
    for name, label in zip(self.names, self.labels):
        cluster = int(label)
        representative = representative_of[cluster]
        similarity = abs(float(self.corr[index[name], index[representative]]))
        if name == representative:
            decision = "keep"
        elif similarity >= min_similarity:
            decision = "drop"
        else:
            decision = "review"
        rows.append(
            {
                "feature": name,
                "cluster": cluster,
                "cluster_size": len(self.clusters[cluster]),
                "kept": name == representative,
                "representative": representative,
                "similarity": similarity,
                "decision": decision,
            }
        )

    return pl.DataFrame(rows).sort(
        ["cluster", "similarity"], descending=[False, True]
    )

PermutationTest dataclass

PermutationTest(observed: float, null: ndarray, p_value: float, alternative: str)

Output of :func:permutation_test.

null_mean property

null_mean: float

Mean of the null distribution — what pure chance produces.

significant property

significant: bool

p_value < 0.05. A convenience, not a licence to stop thinking.

RidgeRegressor

RidgeRegressor(alpha: float = 1.0)

A minimal ridge regression, used as :func:clustered_mda's default model.

Exists so importance scoring works out of the box without pulling scikit-learn into a library whose dependencies are Polars and NumPy. It is deliberately plain — swap in any object exposing fit(X, y) and predict(X) (an sklearn estimator, a gradient-boosted tree) via model_factory when you want something with more capacity.

Ridge rather than ordinary least squares because a feature matrix drawn from this library is routinely singular to machine precision, where OLS has no unique solution at all.

Source code in polars_ta/selection.py

def __init__(self, alpha: float = 1.0) -> None:
    self.alpha = alpha
    self.coef_: np.ndarray | None = None
    self.intercept_: float = 0.0

fit

fit(x: ndarray, y: ndarray) -> RidgeRegressor

Fit on centered data, keeping the intercept out of the penalty.

Source code in polars_ta/selection.py

def fit(self, x: np.ndarray, y: np.ndarray) -> "RidgeRegressor":
    """Fit on centered data, keeping the intercept out of the penalty."""
    x_mean = x.mean(axis=0)
    y_mean = float(y.mean())
    xc, yc = x - x_mean, y - y_mean
    gram = xc.T @ xc + self.alpha * np.eye(x.shape[1])
    self.coef_ = np.linalg.solve(gram, xc.T @ yc)
    self.intercept_ = y_mean - float(x_mean @ self.coef_)
    return self

predict

predict(x: ndarray) -> ndarray

Predict for x.

Source code in polars_ta/selection.py

def predict(self, x: np.ndarray) -> np.ndarray:
    """Predict for `x`."""
    return x @ self.coef_ + self.intercept_

ImportanceResult dataclass

ImportanceResult(clusters: dict[int, list[str]], importance: dict[int, float], std: dict[int, float], baseline: float, n_splits: int)

Output of the clustered importance routines.

Returned by :func:clustered_mda, :func:clustered_mdi, :func:single_feature_importance and :func:clustered_shapley, so the four are directly comparable — same keys, same ordering, same ranked property. What the number means differs by routine (a score drop for MDA and Shapley, an impurity share for MDI, a standalone score for SFI), which is why they should be compared by rank rather than by level.

ranked property

ranked: list[int]

Cluster ids, most important first.

t_stats

t_stats() -> dict[int, float]

Per-cluster mean / (std / sqrt(n_splits)).

A crude significance check across folds. An importance whose sign is not stable fold-to-fold is not an importance — and with few folds this statistic is itself noisy, so read it as a smell test, not a p-value.

Source code in polars_ta/selection.py

def t_stats(self) -> dict[int, float]:
    """Per-cluster `mean / (std / sqrt(n_splits))`.

    A crude significance check across folds. An importance whose sign is
    not stable fold-to-fold is not an importance — and with few folds this
    statistic is itself noisy, so read it as a smell test, not a p-value.
    """
    out = {}
    for cluster, mean in self.importance.items():
        spread = self.std[cluster]
        out[cluster] = (
            0.0 if spread == 0.0 else mean / (spread / math.sqrt(self.n_splits))
        )
    return out

NullImportanceResult dataclass

NullImportanceResult(clusters: dict[int, list[str]], observed: dict[int, float], null: dict[int, ndarray], p_value: dict[int, float])

Output of :func:null_importance.

ranked property

ranked: list[int]

Cluster ids, smallest p-value first, ties broken by importance.

significant

significant(alpha: float = 0.05, correction: str = 'bonferroni') -> set[int]

Cluster ids whose importance beats the null at level alpha.

Parameters:

Name	Type	Description	Default
alpha	float	Family-wise error rate to control.	0.05
correction	str	"bonferroni" (default) divides alpha by the number of clusters tested — the honest default, because testing 20 clusters at 0.05 produces a false positive about 64% of the time. "none" compares each p-value to alpha directly.	'bonferroni'

Returns:

Type	Description
set[int]	The surviving cluster ids.

Raises:

Type	Description
ValueError	If correction is not "bonferroni" or "none".

Source code in polars_ta/selection.py

def significant(
    self, alpha: float = 0.05, correction: str = "bonferroni"
) -> set[int]:
    """Cluster ids whose importance beats the null at level `alpha`.

    Args:
        alpha: Family-wise error rate to control.
        correction: `"bonferroni"` (default) divides `alpha` by the number
            of clusters tested — the honest default, because testing 20
            clusters at 0.05 produces a false positive about 64% of the
            time. `"none"` compares each p-value to `alpha` directly.

    Returns:
        The surviving cluster ids.

    Raises:
        ValueError: If `correction` is not `"bonferroni"` or `"none"`.
    """
    if correction == "bonferroni":
        threshold = alpha / len(self.p_value)
    elif correction == "none":
        threshold = alpha
    else:
        raise ValueError(
            f"correction must be 'bonferroni' or 'none', got {correction!r}"
        )
    return {c for c, p in self.p_value.items() if p < threshold}

report

report() -> DataFrame

Per-cluster table: observed importance, null mean, z-score, p-value.

Columns are cluster, features, observed, null_mean, null_std, z_score and p_value, sorted by ascending p-value. The z_score is (observed - null_mean) / null_std — how many null standard deviations the real importance sits above chance, which stays informative when the p-value floors at 1/(n_draws+1).

Source code in polars_ta/selection.py

def report(self) -> pl.DataFrame:
    """Per-cluster table: observed importance, null mean, z-score, p-value.

    Columns are `cluster`, `features`, `observed`, `null_mean`, `null_std`,
    `z_score` and `p_value`, sorted by ascending p-value. The `z_score` is
    `(observed - null_mean) / null_std` — how many null standard deviations
    the real importance sits above chance, which stays informative when the
    p-value floors at `1/(n_draws+1)`.
    """
    rows = []
    for cluster in sorted(self.observed):
        draws = self.null[cluster]
        spread = float(draws.std())
        rows.append(
            {
                "cluster": cluster,
                "features": ", ".join(self.clusters[cluster]),
                "observed": self.observed[cluster],
                "null_mean": float(draws.mean()),
                "null_std": spread,
                "z_score": (
                    0.0
                    if spread == 0.0
                    else (self.observed[cluster] - float(draws.mean())) / spread
                ),
                "p_value": self.p_value[cluster],
            }
        )
    return pl.DataFrame(rows).sort("p_value")

screen

screen(df: DataFrame, features: Sequence[str] | None = None, max_missing: float | None = None) -> ScreenResult

Turn a frame of indicator columns into a clean, finite feature matrix.

Three kinds of garbage are removed, all of which silently corrupt a correlation matrix if left in:

• Sparse columns (opt-in, via max_missing). Row removal is driven by the worst column, so one feature that is frequently undefined costs every other feature its sample. This is not hypothetical: on the full indicator library, roll_spread is null wherever its serial-covariance premise fails — about 42% of bars, scattered throughout — and on its own it takes a 5,000-bar sample down to 2,610 rows for all 184 features. Setting max_missing=0.2 drops such a column instead of letting it truncate the matrix. worst_missing() on the result names the culprits.
• Warm-up rows. Indicators have wildly different warm-up lengths (an EMA(9) is live after 9 bars, momentum_12_1 after 273). Any row where any surviving feature is null, NaN or infinite is dropped, so every pairwise correlation is computed on the same aligned sample.
• Constant columns. A zero-variance feature has undefined correlation and puts NaN into the matrix, which then poisons the entire eigendecomposition. Also not hypothetical: on a 5,000-bar sample, 21 of the 61 candlestick patterns never fire even once.

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame containing the feature columns.	required
features	Sequence[str] | None	Column names to use. Defaults to every numeric column.	None
max_missing	float | None	Drop any column whose non-finite fraction exceeds this, before dropping rows. None (the default) keeps every column and lets the worst one set the sample size.	None

Returns:

Name	Type	Description
A	ScreenResult	class:ScreenResult.

Raises:

Type	Description
ValueError	If a requested column is missing or non-numeric, if max_missing is outside [0, 1], if fewer than two features survive, or if no rows survive warm-up removal.

Source code in polars_ta/selection.py

def screen(
    df: pl.DataFrame,
    features: Sequence[str] | None = None,
    max_missing: float | None = None,
) -> ScreenResult:
    """Turn a frame of indicator columns into a clean, finite feature matrix.

    Three kinds of garbage are removed, all of which silently corrupt a
    correlation matrix if left in:

    - **Sparse columns** (opt-in, via `max_missing`). Row removal is driven by
      the *worst* column, so one feature that is frequently undefined costs
      every other feature its sample. This is not hypothetical: on the full
      indicator library, `roll_spread` is null wherever its serial-covariance
      premise fails — about 42% of bars, scattered throughout — and on its own
      it takes a 5,000-bar sample down to 2,610 rows for all 184 features.
      Setting `max_missing=0.2` drops such a column instead of letting it
      truncate the matrix. `worst_missing()` on the result names the culprits.
    - **Warm-up rows.** Indicators have wildly different warm-up lengths (an
      EMA(9) is live after 9 bars, `momentum_12_1` after 273). Any row where
      *any* surviving feature is null, NaN or infinite is dropped, so every
      pairwise correlation is computed on the same aligned sample.
    - **Constant columns.** A zero-variance feature has undefined correlation
      and puts `NaN` into the matrix, which then poisons the entire
      eigendecomposition. Also not hypothetical: on a 5,000-bar sample, 21 of
      the 61 candlestick patterns never fire even once.

    Args:
        df: Frame containing the feature columns.
        features: Column names to use. Defaults to every numeric column.
        max_missing: Drop any column whose non-finite fraction exceeds this,
            *before* dropping rows. `None` (the default) keeps every column and
            lets the worst one set the sample size.

    Returns:
        A :class:`ScreenResult`.

    Raises:
        ValueError: If a requested column is missing or non-numeric, if
            `max_missing` is outside `[0, 1]`, if fewer than two features
            survive, or if no rows survive warm-up removal.
    """
    if features is None:
        names = df.select(cs.numeric()).columns
    else:
        names = list(features)
        missing = [c for c in names if c not in df.columns]
        if missing:
            raise ValueError(f"columns not in frame: {missing}")
        non_numeric = [c for c in names if not df.schema[c].is_numeric()]
        if non_numeric:
            raise ValueError(f"non-numeric feature columns: {non_numeric}")

    if len(names) < 2:
        raise ValueError(f"need at least 2 feature columns, got {len(names)}")
    if max_missing is not None and not 0.0 <= max_missing <= 1.0:
        raise ValueError(f"max_missing must be in [0, 1], got {max_missing}")

    values = df.select(names).to_numpy().astype(np.float64)
    finite = np.isfinite(values)
    rates = 1.0 - finite.mean(axis=0)
    missing_rate = {name: float(rate) for name, rate in zip(names, rates)}

    sparse: list[str] = []
    if max_missing is not None:
        dense = rates <= max_missing
        sparse = [n for n, keep_it in zip(names, dense) if not keep_it]
        names = [n for n, keep_it in zip(names, dense) if keep_it]
        values = values[:, dense]
        finite = finite[:, dense]
        if len(names) < 2:
            raise ValueError(
                f"only {len(names)} column(s) survived max_missing="
                f"{max_missing}; dropped as too sparse: {sparse}"
            )

    row_ok = finite.all(axis=1)
    n_rows_dropped = int((~row_ok).sum())
    values = values[row_ok]
    if values.shape[0] == 0:
        worst = max(missing_rate.items(), key=lambda kv: kv[1])
        raise ValueError(
            "no rows survived warm-up removal — every row has a null/NaN/inf "
            f"in at least one feature (worst: {worst[0]} at {worst[1]:.0%} "
            "non-finite). Pass max_missing to drop the sparse columns, shorten "
            "the longest indicator's window, or supply more history"
        )

    keep = values.std(axis=0) > 0.0
    dropped = [n for n, k in zip(names, keep) if not k]
    kept = [n for n, k in zip(names, keep) if k]
    if len(kept) < 2:
        raise ValueError(
            f"only {len(kept)} non-constant feature(s) left after screening; "
            f"dropped as constant: {dropped}"
        )

    return ScreenResult(
        values=values[:, keep],
        names=kept,
        dropped=dropped,
        n_rows_dropped=n_rows_dropped,
        sparse=sparse,
        missing_rate=missing_rate,
    )

nonstationary_features

nonstationary_features(df: DataFrame, features: Sequence[str] | None = None, threshold: float = 0.99) -> list[str]

Flag features whose level behaves like an integrated series.

Correlations between non-stationary series are spurious (Granger & Newbold, 1974) — two independent random walks correlate strongly by construction — so running the rest of this module on raw levels denoises nonsense. Several indicators in this library are integrated by design: on_balance_volume, adi, volume_price_trend, cumulative_return, and every price-level moving average.

The screen is a cheap heuristic, not a unit-root test: lag-1 autocorrelation of the level. A near-unit value means the series barely mean-reverts bar to bar. Fix a flagged feature with :func:polars_ta.quant.frac_diff (keeps most of the memory) or :func:polars_ta.quant.rolling_z_score, then re-run.

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame containing the feature columns.	required
features	Sequence[str] | None	Column names to test. Defaults to every numeric column.	None
threshold	float	Lag-1 autocorrelation at or above which a feature is flagged.	0.99

Returns:

Type	Description
list[str]	The flagged column names, in input order. Columns that are constant (so
list[str]	autocorrelation is undefined) are not flagged — :func:screen drops
list[str]	those separately.

Source code in polars_ta/selection.py

def nonstationary_features(
    df: pl.DataFrame,
    features: Sequence[str] | None = None,
    threshold: float = 0.99,
) -> list[str]:
    """Flag features whose *level* behaves like an integrated series.

    Correlations between non-stationary series are spurious (Granger & Newbold,
    1974) — two independent random walks correlate strongly by construction —
    so running the rest of this module on raw levels denoises nonsense. Several
    indicators in this library are integrated by design: `on_balance_volume`,
    `adi`, `volume_price_trend`, `cumulative_return`, and every price-level
    moving average.

    The screen is a cheap heuristic, **not** a unit-root test: lag-1
    autocorrelation of the level. A near-unit value means the series barely
    mean-reverts bar to bar. Fix a flagged feature with
    :func:`polars_ta.quant.frac_diff` (keeps most of the memory) or
    :func:`polars_ta.quant.rolling_z_score`, then re-run.

    Args:
        df: Frame containing the feature columns.
        features: Column names to test. Defaults to every numeric column.
        threshold: Lag-1 autocorrelation at or above which a feature is flagged.

    Returns:
        The flagged column names, in input order. Columns that are constant (so
        autocorrelation is undefined) are not flagged — :func:`screen` drops
        those separately.
    """
    names = df.select(cs.numeric()).columns if features is None else list(features)
    flagged = []
    for name in names:
        col = df[name].cast(pl.Float64).to_numpy()
        col = col[np.isfinite(col)]
        if col.size < 3:
            continue
        current, lagged = col[1:], col[:-1]
        if current.std() == 0.0 or lagged.std() == 0.0:
            continue
        if float(np.corrcoef(current, lagged)[0, 1]) >= threshold:
            flagged.append(name)
    return flagged

corr_matrix

corr_matrix(values: ndarray, method: str = 'spearman') -> ndarray

Correlation matrix of a (n_obs, n_features) matrix.

Parameters:

Name	Type	Description	Default
values	ndarray	Finite feature matrix, e.g. ScreenResult.values.	required
method	str	"spearman" (default) correlates ranks, "pearson" correlates levels. Spearman is the better default here: bounded oscillators and heavy-tailed microstructure features break Pearson's linearity assumption, and rank correlation is invariant to the monotone rescalings indicators are full of.	'spearman'

Returns:

Type	Description
ndarray	A symmetric (n_features, n_features) matrix with a unit diagonal.

Raises:

Type	Description
ValueError	If method is not "spearman" or "pearson".

Source code in polars_ta/selection.py

def corr_matrix(values: np.ndarray, method: str = "spearman") -> np.ndarray:
    """Correlation matrix of a ``(n_obs, n_features)`` matrix.

    Args:
        values: Finite feature matrix, e.g. `ScreenResult.values`.
        method: `"spearman"` (default) correlates *ranks*, `"pearson"`
            correlates levels. Spearman is the better default here: bounded
            oscillators and heavy-tailed microstructure features break
            Pearson's linearity assumption, and rank correlation is invariant
            to the monotone rescalings indicators are full of.

    Returns:
        A symmetric ``(n_features, n_features)`` matrix with a unit diagonal.

    Raises:
        ValueError: If `method` is not `"spearman"` or `"pearson"`.
    """
    if method == "spearman":
        values = np.apply_along_axis(_average_ranks, 0, values)
    elif method != "pearson":
        raise ValueError(f"method must be 'spearman' or 'pearson', got {method!r}")

    corr = np.corrcoef(values, rowvar=False)
    corr = np.clip(corr, -1.0, 1.0)
    corr = (corr + corr.T) / 2.0
    np.fill_diagonal(corr, 1.0)
    return corr

vif

vif(corr: ndarray, ridge: float = 1e-10) -> ndarray

Variance inflation factor per feature.

\(\mathrm{VIF}_i = 1/(1 - R_i^2)\), where \(R_i^2\) is from regressing feature \(i\) on all the others — equivalently the \(i\)-th diagonal entry of the inverted correlation matrix. It answers a different question from a pairwise correlation: a feature can correlate weakly with every other feature individually and still be an exact linear combination of several of them, which pairwise screening never sees. macd versus a pair of EMAs is exactly that case.

Conventional reading: VIF > 5 is worth a look, > 10 is severe. On a wide indicator matrix expect values in the thousands — which is the finding, not a bug.

Computed by eigendecomposition as \(\mathrm{VIF}_i = \sum_j V_{ij}^2 / (\lambda_j + \varepsilon)\) rather than by inverting the matrix. A wide indicator matrix is routinely singular to machine precision (the full library reaches a condition number of ~1e21), where inv returns noise and — worse — pinv quietly hides the problem: the pseudo-inverse discards the null space, so exactly-collinear features come back with a reassuring VIF of 1.0. The ridge keeps the blow-up visible and finite instead.

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal.	required
ridge	float	Added to each eigenvalue before inversion. Sets the ceiling on a reported VIF (roughly 1/ridge) and is what makes a perfectly collinear feature report a huge number rather than raising.	1e-10

Returns:

Type	Description
ndarray	One VIF per feature, aligned with corr's rows, each at least 1.

Source code in polars_ta/selection.py

def vif(corr: np.ndarray, ridge: float = 1e-10) -> np.ndarray:
    r"""Variance inflation factor per feature.

    $\mathrm{VIF}_i = 1/(1 - R_i^2)$, where $R_i^2$ is from regressing feature
    $i$ on all the others — equivalently the $i$-th diagonal entry of the
    inverted correlation matrix. It answers a different question from a
    pairwise correlation: a feature can correlate weakly with *every* other
    feature individually and still be an exact linear combination of several
    of them, which pairwise screening never sees. `macd` versus a pair of EMAs
    is exactly that case.

    Conventional reading: VIF > 5 is worth a look, > 10 is severe. On a wide
    indicator matrix expect values in the thousands — which is the finding, not
    a bug.

    Computed by eigendecomposition as
    $\mathrm{VIF}_i = \sum_j V_{ij}^2 / (\lambda_j + \varepsilon)$ rather than
    by inverting the matrix. A wide indicator matrix is routinely singular to
    machine precision (the full library reaches a condition number of ~1e21),
    where `inv` returns noise and — worse — `pinv` quietly *hides* the problem:
    the pseudo-inverse discards the null space, so exactly-collinear features
    come back with a reassuring VIF of 1.0. The ridge keeps the blow-up
    visible and finite instead.

    Args:
        corr: Symmetric correlation matrix with a unit diagonal.
        ridge: Added to each eigenvalue before inversion. Sets the ceiling on
            a reported VIF (roughly `1/ridge`) and is what makes a perfectly
            collinear feature report a huge number rather than raising.

    Returns:
        One VIF per feature, aligned with `corr`'s rows, each at least 1.
    """
    corr = _check_corr(corr)
    eigvals, eigvecs = np.linalg.eigh(corr)
    inverted = (eigvecs**2) / (np.clip(eigvals, 0.0, None) + ridge)
    return np.clip(inverted.sum(axis=1), 1.0, None)

mutual_information

mutual_information(values: ndarray, bins: int | None = None, normalize: bool = True) -> ndarray

Pairwise mutual information matrix.

\(I(X;Y) = H(X) + H(Y) - H(X,Y)\), estimated from a 2-D histogram. Unlike correlation it is not limited to monotone association, and it handles the sparse discrete features (candlestick patterns) that make a rank correlation near-useless.

Parameters:

Name	Type	Description	Default
values	ndarray	Finite (n_obs, n_features) feature matrix.	required
bins	int | None	Histogram bins per axis. Defaults to the Hacine-Gharbi optimum for the sample size.	None
normalize	bool	Divide by min(H(X), H(Y)) so the result lands in [0, 1] and is comparable across pairs, with 1 meaning one feature determines the other. Raw nats when False.	True

Returns:

Type	Description
ndarray	A symmetric (n_features, n_features) matrix. The diagonal is 1 when
ndarray	normalized (a feature determines itself), else its own entropy.

Source code in polars_ta/selection.py

def mutual_information(
    values: np.ndarray, bins: int | None = None, normalize: bool = True
) -> np.ndarray:
    """Pairwise mutual information matrix.

    $I(X;Y) = H(X) + H(Y) - H(X,Y)$, estimated from a 2-D histogram. Unlike
    correlation it is **not** limited to monotone association, and it handles
    the sparse discrete features (candlestick patterns) that make a rank
    correlation near-useless.

    Args:
        values: Finite `(n_obs, n_features)` feature matrix.
        bins: Histogram bins per axis. Defaults to the Hacine-Gharbi optimum
            for the sample size.
        normalize: Divide by `min(H(X), H(Y))` so the result lands in `[0, 1]`
            and is comparable across pairs, with 1 meaning one feature
            determines the other. Raw nats when `False`.

    Returns:
        A symmetric `(n_features, n_features)` matrix. The diagonal is 1 when
        normalized (a feature determines itself), else its own entropy.
    """
    values = np.asarray(values, dtype=np.float64)
    n_obs, n_features = values.shape
    bins = _bin_count(n_obs) if bins is None else bins

    out = np.zeros((n_features, n_features))
    for i in range(n_features):
        for j in range(i, n_features):
            h_x, h_y, h_xy = _entropies(values[:, i], values[:, j], bins)
            info = h_x + h_y - h_xy
            if normalize:
                floor = min(h_x, h_y)
                # A constant column has zero entropy and shares no information
                # with anything; call that 0 rather than 0/0.
                info = 0.0 if floor <= 0.0 else info / floor
            out[i, j] = out[j, i] = info
    if normalize:
        np.fill_diagonal(out, 1.0)
    return out

variation_of_information

variation_of_information(values: ndarray, bins: int | None = None) -> ndarray

Normalized variation of information — a metric distance.

\(VI(X;Y) = H(X,Y) - I(X;Y)\), divided by \(H(X,Y)\) to land in \([0, 1]\). It satisfies the triangle inequality (Meila, 2007), which is what makes it a legitimate input to :func:cluster_by_distance — and it is the non-linear counterpart to :func:distance_matrix: 0 when two features determine each other, 1 when they are independent.

Use it instead of the correlation distance when the feature set contains discrete or sparse columns, or when you suspect non-monotone relationships that Spearman scores as zero.

Parameters:

Name	Type	Description	Default
values	ndarray	Finite (n_obs, n_features) feature matrix.	required
bins	int | None	Histogram bins per axis; defaults to the Hacine-Gharbi optimum.	None

Returns:

Type	Description
ndarray	A symmetric distance matrix with a zero diagonal.

Source code in polars_ta/selection.py

def variation_of_information(values: np.ndarray, bins: int | None = None) -> np.ndarray:
    r"""Normalized variation of information — a **metric** distance.

    $VI(X;Y) = H(X,Y) - I(X;Y)$, divided by $H(X,Y)$ to land in $[0, 1]$. It
    satisfies the triangle inequality (Meila, 2007), which is what makes it a
    legitimate input to :func:`cluster_by_distance` — and it is the non-linear
    counterpart to :func:`distance_matrix`: 0 when two features determine each
    other, 1 when they are independent.

    Use it instead of the correlation distance when the feature set contains
    discrete or sparse columns, or when you suspect non-monotone relationships
    that Spearman scores as zero.

    Args:
        values: Finite `(n_obs, n_features)` feature matrix.
        bins: Histogram bins per axis; defaults to the Hacine-Gharbi optimum.

    Returns:
        A symmetric distance matrix with a zero diagonal.
    """
    values = np.asarray(values, dtype=np.float64)
    n_obs, n_features = values.shape
    bins = _bin_count(n_obs) if bins is None else bins

    out = np.zeros((n_features, n_features))
    for i in range(n_features):
        for j in range(i + 1, n_features):
            h_x, h_y, h_xy = _entropies(values[:, i], values[:, j], bins)
            info = h_x + h_y - h_xy
            # Two constant columns share a zero-entropy joint distribution;
            # they are identical, so the distance is 0.
            distance = 0.0 if h_xy <= 0.0 else (h_xy - info) / h_xy
            out[i, j] = out[j, i] = min(max(distance, 0.0), 1.0)
    return out

marcenko_pastur_edge

marcenko_pastur_edge(q: float, sigma2: float = 1.0) -> float

Upper edge of the Marchenko-Pastur eigenvalue support.

For a correlation matrix built from n_obs observations of n_features independent variables with variance \(\sigma^2\), the eigenvalues concentrate on \([\lambda_-, \lambda_+]\) with

\[\lambda_\pm = \sigma^2 (1 \pm \sqrt{1/q})^2, \quad q = T/N.\]

Any empirical eigenvalue below \(\lambda_+\) is therefore consistent with pure noise. With 205 features over 5000 bars, \(q \approx 24\) and \(\lambda_+ \approx 1.45\) — a sobering threshold, since the average eigenvalue of any correlation matrix is exactly 1.

Parameters:

Name	Type	Description	Default
q	float	Aspect ratio n_obs / n_features, which must exceed 1.	required
sigma2	float	Noise variance.	1.0

Returns:

Type	Description
float	The upper edge \(\lambda_+\).

Raises:

Type	Description
ValueError	If q <= 1 — with fewer observations than features the sample correlation matrix is singular and this whole approach is not applicable.

Source code in polars_ta/selection.py

def marcenko_pastur_edge(q: float, sigma2: float = 1.0) -> float:
    r"""Upper edge of the Marchenko-Pastur eigenvalue support.

    For a correlation matrix built from `n_obs` observations of `n_features`
    *independent* variables with variance $\sigma^2$, the eigenvalues
    concentrate on $[\lambda_-, \lambda_+]$ with

    $$\lambda_\pm = \sigma^2 (1 \pm \sqrt{1/q})^2, \quad q = T/N.$$

    Any empirical eigenvalue below $\lambda_+$ is therefore consistent
    with pure noise. With 205 features over 5000 bars, $q \approx 24$ and
    $\lambda_+ \approx 1.45$ — a sobering threshold, since the average
    eigenvalue of any correlation matrix is exactly 1.

    Args:
        q: Aspect ratio `n_obs / n_features`, which must exceed 1.
        sigma2: Noise variance.

    Returns:
        The upper edge $\lambda_+$.

    Raises:
        ValueError: If `q <= 1` — with fewer observations than features the
            sample correlation matrix is singular and this whole approach is
            not applicable.
    """
    if q <= 1.0:
        raise ValueError(
            f"need more observations than features (q = T/N > 1), got q = {q:.3f}"
        )
    return sigma2 * (1.0 + math.sqrt(1.0 / q)) ** 2

denoise

denoise(corr: ndarray, q: float, detone: int = 0, bandwidth: float = 0.01) -> ndarray

Marchenko-Pastur denoising of a correlation matrix.

Eigenvalues below the MP edge carry no information beyond what a matrix of independent noise would produce, but they are the smallest ones — so anything that inverts or optimizes over the raw matrix is dominated by them. Clipping replaces the whole noise bulk with its average, which kills the spurious structure while preserving the trace (and therefore the unit diagonal, up to the rescaling applied at the end).

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal.	required
q	float	Aspect ratio n_obs / n_features used for the MP edge.	required
detone	int	Number of leading (largest) eigenvalues to remove entirely. detone=1 strips the dominant mode, which for technical indicators on one asset is essentially "trend/level" and swamps everything else — worth doing when you want features that say something beyond direction. Off by default because it is a modelling choice, not a correction.	0
bandwidth	float	Kernel bandwidth for the eigenvalue density fit that estimates the noise variance. The default suits the hundreds of eigenvalues a wide indicator matrix produces; widen it for a narrow matrix, where the empirical density is spiky.	0.01

Returns:

Type	Description
ndarray	The denoised correlation matrix, symmetric with a unit diagonal.

Raises:

Type	Description
ValueError	If corr is not a valid correlation matrix, if q <= 1, or if detone is outside [0, n_features).

Source code in polars_ta/selection.py

def denoise(
    corr: np.ndarray, q: float, detone: int = 0, bandwidth: float = 0.01
) -> np.ndarray:
    """Marchenko-Pastur denoising of a correlation matrix.

    Eigenvalues below the MP edge carry no information beyond what a matrix of
    independent noise would produce, but they are the *smallest* ones — so
    anything that inverts or optimizes over the raw matrix is dominated by
    them. Clipping replaces the whole noise bulk with its average, which kills
    the spurious structure while preserving the trace (and therefore the unit
    diagonal, up to the rescaling applied at the end).

    Args:
        corr: Symmetric correlation matrix with a unit diagonal.
        q: Aspect ratio `n_obs / n_features` used for the MP edge.
        detone: Number of leading (largest) eigenvalues to remove entirely.
            `detone=1` strips the dominant mode, which for technical indicators
            on one asset is essentially "trend/level" and swamps everything
            else — worth doing when you want features that say something
            *beyond* direction. Off by default because it is a modelling
            choice, not a correction.
        bandwidth: Kernel bandwidth for the eigenvalue density fit that
            estimates the noise variance. The default suits the hundreds of
            eigenvalues a wide indicator matrix produces; widen it for a
            narrow matrix, where the empirical density is spiky.

    Returns:
        The denoised correlation matrix, symmetric with a unit diagonal.

    Raises:
        ValueError: If `corr` is not a valid correlation matrix, if `q <= 1`,
            or if `detone` is outside `[0, n_features)`.
    """
    corr = _check_corr(corr)
    n = corr.shape[0]
    if not 0 <= detone < n:
        raise ValueError(f"detone must be in [0, {n}), got {detone}")

    eigvals, eigvecs = np.linalg.eigh(corr)  # ascending
    edge = marcenko_pastur_edge(q, _noise_variance(eigvals, q, bandwidth))

    clipped = eigvals.copy()
    noise = eigvals <= edge
    # An all-signal matrix leaves nothing to clip; the assignment is then a
    # no-op on an empty selection, and `mean` of an empty slice is avoided.
    clipped[noise] = eigvals[noise].mean() if noise.any() else 0.0
    if detone > 0:
        clipped[n - detone :] = 0.0

    out = (eigvecs * clipped) @ eigvecs.T
    # Rescale back to a correlation matrix: clipping preserves the trace but
    # not each diagonal entry, and detoning shrinks them outright.
    scale = np.sqrt(np.clip(np.diag(out), 1e-12, None))
    out = out / np.outer(scale, scale)
    out = np.clip((out + out.T) / 2.0, -1.0, 1.0)
    np.fill_diagonal(out, 1.0)
    return out

signal_rank

signal_rank(corr: ndarray, q: float, bandwidth: float = 0.01) -> int

Number of eigenvalues above the Marchenko-Pastur edge.

A first answer to "how many features do I actually have?". For a matrix of pure independent noise this is 0; for 205 technical indicators on one asset, expect somewhere in the teens.

Read this as a ceiling, not a measurement

The Marchenko-Pastur law assumes i.i.d. rows, and rolling indicators violate that badly — a 14-period ATR has lag-1 autocorrelation around 0.997, so n_obs wildly overstates the effective sample size and the edge lands too low.

Measured on the BTCUSDT 5m fixture with 16 indicators: the observed count is 4, but a :func:block_permute null that keeps each feature's autocorrelation and destroys only the cross-feature alignment produces 5.1 ± 1.0 — so the count itself is not distinguishable from chance (p ≈ 1.0). The same null leaves the top eigenvalue (6.07 vs 1.54) and :func:effective_rank (5.66 vs 15.5) overwhelmingly significant.

In other words: the existence of strong shared structure is real, and the eigenvalues far above the edge are real. The precise count of those hovering at it is not. Use :func:permutation_test to find out which regime you are in rather than trusting this number on its own, and prefer :func:effective_rank as the headline figure.

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal.	required
q	float	Aspect ratio n_obs / n_features.	required
bandwidth	float	Kernel bandwidth for the noise-variance fit; see :func:denoise.	0.01

Returns:

Type	Description
int	The count of eigenvalues strictly above the estimated edge.

Source code in polars_ta/selection.py

def signal_rank(corr: np.ndarray, q: float, bandwidth: float = 0.01) -> int:
    """Number of eigenvalues above the Marchenko-Pastur edge.

    A first answer to "how many features do I actually have?". For a matrix of
    pure independent noise this is 0; for 205 technical indicators on one asset,
    expect somewhere in the teens.

    !!! warning "Read this as a ceiling, not a measurement"
        The Marchenko-Pastur law assumes **i.i.d. rows**, and rolling
        indicators violate that badly — a 14-period ATR has lag-1
        autocorrelation around 0.997, so `n_obs` wildly overstates the
        *effective* sample size and the edge lands too low.

        Measured on the BTCUSDT 5m fixture with 16 indicators: the observed
        count is 4, but a :func:`block_permute` null that keeps each feature's
        autocorrelation and destroys only the cross-feature alignment produces
        **5.1 ± 1.0** — so the count itself is not distinguishable from chance
        (p ≈ 1.0). The same null leaves the top eigenvalue (6.07 vs 1.54) and
        :func:`effective_rank` (5.66 vs 15.5) overwhelmingly significant.

        In other words: the *existence* of strong shared structure is real, and
        the eigenvalues far above the edge are real. The precise count of those
        hovering *at* it is not. Use :func:`permutation_test` to find out which
        regime you are in rather than trusting this number on its own, and
        prefer :func:`effective_rank` as the headline figure.

    Args:
        corr: Symmetric correlation matrix with a unit diagonal.
        q: Aspect ratio `n_obs / n_features`.
        bandwidth: Kernel bandwidth for the noise-variance fit; see
            :func:`denoise`.

    Returns:
        The count of eigenvalues strictly above the estimated edge.
    """
    corr = _check_corr(corr)
    eigvals = np.linalg.eigvalsh(corr)
    edge = marcenko_pastur_edge(q, _noise_variance(eigvals, q, bandwidth))
    return int((eigvals > edge).sum())

effective_rank

effective_rank(corr: ndarray) -> float

Entropy-based effective rank (Roy & Vetterli, 2007).

\(\exp(H)\) where \(H = -\sum_i p_i \ln p_i\) and \(p_i = \lambda_i / \sum_j \lambda_j\). A continuous companion to :func:signal_rank that needs no noise model: it is n_features when the features are perfectly independent and 1 when they are perfectly collinear, so the ratio to n_features reads directly as "fraction of the matrix that is genuinely distinct".

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal.	required

Returns:

Type	Description
float	The effective rank, in [1, n_features].

Source code in polars_ta/selection.py

def effective_rank(corr: np.ndarray) -> float:
    r"""Entropy-based effective rank (Roy & Vetterli, 2007).

    $\exp(H)$ where $H = -\sum_i p_i \ln p_i$ and
    $p_i = \lambda_i / \sum_j \lambda_j$. A continuous companion to
    :func:`signal_rank` that needs no noise model: it is `n_features` when the
    features are perfectly independent and 1 when they are perfectly
    collinear, so the ratio to `n_features` reads directly as "fraction of the
    matrix that is genuinely distinct".

    Args:
        corr: Symmetric correlation matrix with a unit diagonal.

    Returns:
        The effective rank, in ``[1, n_features]``.
    """
    corr = _check_corr(corr)
    eigvals = np.clip(np.linalg.eigvalsh(corr), 0.0, None)
    p = eigvals / eigvals.sum()
    p = p[p > 0.0]
    return float(np.exp(-(p * np.log(p)).sum()))

distance_matrix

distance_matrix(corr: ndarray) -> ndarray

Correlation-to-distance metric \(d_{ij} = \sqrt{(1-\rho_{ij})/2}\).

Unlike \(1-\rho\), this is a true metric (Mantegna, 1999) — it satisfies the triangle inequality, which is what makes hierarchical clustering on it meaningful rather than merely suggestive. It maps \(\rho = 1\) to 0, \(\rho = 0\) to \(1/\sqrt{2}\), and \(\rho = -1\) to 1.

Note that perfectly anti-correlated features land far apart, even though they carry the same information with a flipped sign. If that matters for your feature set, pass np.abs(corr) in.

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal.	required

Returns:

Type	Description
ndarray	A symmetric distance matrix with a zero diagonal.

Source code in polars_ta/selection.py

def distance_matrix(corr: np.ndarray) -> np.ndarray:
    r"""Correlation-to-distance metric $d_{ij} = \sqrt{(1-\rho_{ij})/2}$.

    Unlike $1-\rho$, this is a *true metric* (Mantegna, 1999) — it
    satisfies the triangle inequality, which is what makes hierarchical
    clustering on it meaningful rather than merely suggestive. It maps
    $\rho = 1$ to 0, $\rho = 0$ to $1/\sqrt{2}$, and
    $\rho = -1$ to 1.

    Note that perfectly anti-correlated features land far apart, even though
    they carry the same information with a flipped sign. If that matters for
    your feature set, pass `np.abs(corr)` in.

    Args:
        corr: Symmetric correlation matrix with a unit diagonal.

    Returns:
        A symmetric distance matrix with a zero diagonal.
    """
    corr = _check_corr(corr)
    dist = np.sqrt(np.clip((1.0 - corr) / 2.0, 0.0, None))
    np.fill_diagonal(dist, 0.0)
    return dist

silhouette_score

silhouette_score(dist: ndarray, labels: ndarray) -> float

Mean silhouette of a partition (Rousseeuw, 1987).

For each feature, a is its mean distance to its own cluster and b its mean distance to the nearest other cluster; the silhouette is (b - a) / max(a, b). Near 1 means the clusters are tight and well-separated, near 0 means they overlap, negative means features are closer to a neighbouring cluster than their own.

Members of singleton clusters score 0, the usual convention — a lone feature is neither well nor badly placed.

Parameters:

Name	Type	Description	Default
dist	ndarray	Symmetric distance matrix.	required
labels	ndarray	Cluster id per feature.	required

Returns:

Type	Description
float	The mean silhouette, in [-1, 1]. Returns 0.0 for a degenerate
float	partition (one cluster, or every feature in its own).

Source code in polars_ta/selection.py

def silhouette_score(dist: np.ndarray, labels: np.ndarray) -> float:
    """Mean silhouette of a partition (Rousseeuw, 1987).

    For each feature, `a` is its mean distance to its own cluster and `b` its
    mean distance to the nearest *other* cluster; the silhouette is
    `(b - a) / max(a, b)`. Near 1 means the clusters are tight and
    well-separated, near 0 means they overlap, negative means features are
    closer to a neighbouring cluster than their own.

    Members of singleton clusters score 0, the usual convention — a lone
    feature is neither well nor badly placed.

    Args:
        dist: Symmetric distance matrix.
        labels: Cluster id per feature.

    Returns:
        The mean silhouette, in ``[-1, 1]``. Returns 0.0 for a degenerate
        partition (one cluster, or every feature in its own).
    """
    n = len(labels)
    unique = np.unique(labels)
    if unique.size < 2 or unique.size == n:
        return 0.0

    scores = np.zeros(n)
    masks = {c: labels == c for c in unique}
    for i in range(n):
        own = masks[labels[i]]
        n_own = int(own.sum())
        if n_own == 1:
            continue  # singleton: silhouette 0 by convention
        a = dist[i, own].sum() / (n_own - 1)
        b = min(dist[i, masks[c]].mean() for c in unique if c != labels[i])
        spread = max(a, b)
        # All-identical features give a == b == 0; call that 0, not 0/0.
        scores[i] = 0.0 if spread == 0.0 else (b - a) / spread
    return float(scores.mean())

cluster_features

cluster_features(corr: ndarray, max_k: int | None = None, k: int | None = None) -> tuple[ndarray, float]

Cluster features by correlation distance, picking k by silhouette.

Average-linkage agglomerative clustering on :func:distance_matrix, cut at every k from 2 to max_k and scored by :func:silhouette_score; the best-scoring cut wins. This is the correlation-clustering step of López de Prado's ONC, minus the recursive re-clustering of low-quality clusters.

Parameters:

Name	Type	Description	Default
corr	ndarray	Symmetric correlation matrix with a unit diagonal — denoise it first, otherwise you are clustering on noise.	required
max_k	int | None	Largest number of clusters to consider. Defaults to n_features - 1. Ignored when k is given.	None
k	int | None	Force an exact number of clusters instead of choosing by silhouette. The natural argument is :func:signal_rank, which says how many dimensions the matrix actually has — silhouette answers a different question ("which cut is tidiest?") and on a smooth correlation structure it often settles on a very coarse 2.	None

Returns:

Type	Description
ndarray	(labels, silhouette) — the cluster id per feature and the mean
float	silhouette of the chosen partition.

Raises:

Type	Description
ValueError	If there are fewer than 3 features, if max_k < 2, or if k is outside [2, n_features].

Source code in polars_ta/selection.py

def cluster_features(
    corr: np.ndarray, max_k: int | None = None, k: int | None = None
) -> tuple[np.ndarray, float]:
    """Cluster features by correlation distance, picking `k` by silhouette.

    Average-linkage agglomerative clustering on
    :func:`distance_matrix`, cut at every `k` from 2 to `max_k` and scored by
    :func:`silhouette_score`; the best-scoring cut wins. This is the
    correlation-clustering step of López de Prado's ONC, minus the recursive
    re-clustering of low-quality clusters.

    Args:
        corr: Symmetric correlation matrix with a unit diagonal — denoise it
            first, otherwise you are clustering on noise.
        max_k: Largest number of clusters to consider. Defaults to
            `n_features - 1`. Ignored when `k` is given.
        k: Force an exact number of clusters instead of choosing by silhouette.
            The natural argument is :func:`signal_rank`, which says how many
            dimensions the matrix actually has — silhouette answers a different
            question ("which cut is tidiest?") and on a smooth correlation
            structure it often settles on a very coarse 2.

    Returns:
        `(labels, silhouette)` — the cluster id per feature and the mean
        silhouette of the chosen partition.

    Raises:
        ValueError: If there are fewer than 3 features, if `max_k < 2`, or if
            `k` is outside `[2, n_features]`.
    """
    return cluster_by_distance(distance_matrix(corr), max_k=max_k, k=k)

cluster_by_distance

cluster_by_distance(dist: ndarray, max_k: int | None = None, k: int | None = None) -> tuple[ndarray, float]

Cluster on any feature distance matrix.

The engine behind :func:cluster_features, exposed separately so the distance need not come from a correlation. Pass :func:variation_of_information to cluster on non-linear dependence instead — the right choice when the feature set includes the sparse, discrete candlestick patterns, which a rank correlation barely sees.

Parameters:

Name	Type	Description	Default
dist	ndarray	Symmetric distance matrix with a zero diagonal. It should be a true metric; the correlation distance and normalized variation of information both are.	required
max_k	int | None	Largest number of clusters to consider. Defaults to n_features - 1. Ignored when k is given.	None
k	int | None	Force an exact number of clusters instead of choosing by silhouette.	None

Returns:

Type	Description
ndarray	(labels, silhouette) — the cluster id per feature and the mean
float	silhouette of the chosen partition.

Raises:

Type	Description
ValueError	If dist is not a square, symmetric, zero-diagonal matrix, if there are fewer than 3 features, if max_k < 2, or if k is outside [2, n_features].

Source code in polars_ta/selection.py

def cluster_by_distance(
    dist: np.ndarray, max_k: int | None = None, k: int | None = None
) -> tuple[np.ndarray, float]:
    """Cluster on **any** feature distance matrix.

    The engine behind :func:`cluster_features`, exposed separately so the
    distance need not come from a correlation. Pass
    :func:`variation_of_information` to cluster on non-linear dependence
    instead — the right choice when the feature set includes the sparse,
    discrete candlestick patterns, which a rank correlation barely sees.

    Args:
        dist: Symmetric distance matrix with a zero diagonal. It should be a
            true metric; the correlation distance and normalized variation of
            information both are.
        max_k: Largest number of clusters to consider. Defaults to
            `n_features - 1`. Ignored when `k` is given.
        k: Force an exact number of clusters instead of choosing by silhouette.

    Returns:
        `(labels, silhouette)` — the cluster id per feature and the mean
        silhouette of the chosen partition.

    Raises:
        ValueError: If `dist` is not a square, symmetric, zero-diagonal matrix,
            if there are fewer than 3 features, if `max_k < 2`, or if `k` is
            outside `[2, n_features]`.
    """
    dist = np.asarray(dist, dtype=np.float64)
    if dist.ndim != 2 or dist.shape[0] != dist.shape[1]:
        raise ValueError(f"dist must be square, got shape {dist.shape}")
    if not np.allclose(dist, dist.T):
        raise ValueError("dist must be symmetric")
    if not np.allclose(np.diag(dist), 0.0):
        raise ValueError("dist must have a zero diagonal")

    n = dist.shape[0]
    if n < 3:
        raise ValueError(f"need at least 3 features to cluster, got {n}")

    labelings = _average_linkage(dist)

    if k is not None:
        if not 2 <= k <= n:
            raise ValueError(f"k must be in [2, {n}], got {k}")
        return labelings[k], silhouette_score(dist, labelings[k])

    max_k = n - 1 if max_k is None else min(max_k, n - 1)
    if max_k < 2:
        raise ValueError(f"max_k must be at least 2, got {max_k}")

    best_labels = labelings[2]
    best_score = -np.inf
    for cut in range(2, max_k + 1):
        score = silhouette_score(dist, labelings[cut])
        if score > best_score:
            best_score, best_labels = score, labelings[cut]
    return best_labels, float(best_score)

select_representatives

select_representatives(dist: ndarray, labels: ndarray, names: Sequence[str], scores: Sequence[float] | None = None) -> list[str]

Reduce each cluster to a single feature.

Parameters:

Name	Type	Description	Default
dist	ndarray	Symmetric distance matrix.	required
labels	ndarray	Cluster id per feature.	required
names	Sequence[str]	Feature names, aligned with labels.	required
scores	Sequence[float] | None	Optional per-feature quality score. When given, each cluster elects its highest-|score| member — pass abs information coefficients (see :func:polars_ta.quant.rolling_ic) to keep the most predictive member of each group. When omitted, each cluster elects its medoid: the member with the smallest total distance to the rest, i.e. the most representative one. The medoid is the safe default because it never looks at a target and so cannot leak a label.	None

Returns:

Type	Description
list[str]	One name per cluster, ordered by cluster id.

Source code in polars_ta/selection.py

def select_representatives(
    dist: np.ndarray,
    labels: np.ndarray,
    names: Sequence[str],
    scores: Sequence[float] | None = None,
) -> list[str]:
    """Reduce each cluster to a single feature.

    Args:
        dist: Symmetric distance matrix.
        labels: Cluster id per feature.
        names: Feature names, aligned with `labels`.
        scores: Optional per-feature quality score. When given, each cluster
            elects its **highest-|score|** member — pass `abs` information
            coefficients (see :func:`polars_ta.quant.rolling_ic`) to keep the
            most predictive member of each group. When omitted, each cluster
            elects its **medoid**: the member with the smallest total distance
            to the rest, i.e. the most representative one. The medoid is the
            safe default because it never looks at a target and so cannot leak
            a label.

    Returns:
        One name per cluster, ordered by cluster id.
    """
    picked = []
    for cluster in np.unique(labels):
        members = np.flatnonzero(labels == cluster)
        if scores is None:
            within = dist[np.ix_(members, members)].sum(axis=1)
            choice = members[int(np.argmin(within))]
        else:
            quality = np.abs(np.asarray(scores, dtype=np.float64)[members])
            choice = members[int(np.argmax(quality))]
        picked.append(names[choice])
    return picked

select_features

select_features(df: DataFrame, features: Sequence[str] | None = None, method: str = 'spearman', detone: int = 0, max_k: int | None = None, k: int | None = None, scores: Mapping[str, float] | None = None, bandwidth: float = 0.01) -> SelectionResult

Screen, correlate, denoise, cluster, and elect one feature per cluster.

The end-to-end redundancy pass. No target is involved, so nothing here can leak a label — run it first, then do model-based importance (with purged cross-validation) on the survivors, where the substitution effect can no longer scramble the scores.

import polars as pl
from polars_ta import momentum, trend, volatility, selection

feats = df.with_columns(
    momentum.rsi("close").alias("rsi"),
    momentum.roc("close").alias("roc"),
    trend.macd("close").alias("macd"),
    volatility.average_true_range("high", "low", "close").alias("atr"),
)

result = selection.select_features(feats, ["rsi", "roc", "macd", "atr"])
print(result.selected)      # one feature per cluster
print(result.signal_rank)   # how many dimensions there really were

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame containing the computed indicator columns.	required
features	Sequence[str] | None	Column names to consider. Defaults to every numeric column — check that this excludes raw OHLCV and any target column.	None
method	str	Correlation method, "spearman" (default) or "pearson".	'spearman'
detone	int	Leading eigenvalues to strip; see :func:denoise.	0
max_k	int | None	Largest number of clusters to consider.	None
k	int | None	Force an exact number of clusters; see :func:cluster_features.	None
scores	Mapping[str, float] | None	Optional {feature_name: score} used to elect cluster representatives by |score| instead of by medoid; see :func:select_representatives. Names not present score 0.	None
bandwidth	float	Kernel bandwidth for the noise-variance fit; see :func:denoise.	0.01

Returns:

Name	Type	Description
A	SelectionResult	class:SelectionResult.

Source code in polars_ta/selection.py

def select_features(
    df: pl.DataFrame,
    features: Sequence[str] | None = None,
    method: str = "spearman",
    detone: int = 0,
    max_k: int | None = None,
    k: int | None = None,
    scores: Mapping[str, float] | None = None,
    bandwidth: float = 0.01,
) -> SelectionResult:
    """Screen, correlate, denoise, cluster, and elect one feature per cluster.

    The end-to-end redundancy pass. No target is involved, so nothing here can
    leak a label — run it first, then do model-based importance (with purged
    cross-validation) on the survivors, where the substitution effect can no
    longer scramble the scores.

    ```python
    import polars as pl
    from polars_ta import momentum, trend, volatility, selection

    feats = df.with_columns(
        momentum.rsi("close").alias("rsi"),
        momentum.roc("close").alias("roc"),
        trend.macd("close").alias("macd"),
        volatility.average_true_range("high", "low", "close").alias("atr"),
    )

    result = selection.select_features(feats, ["rsi", "roc", "macd", "atr"])
    print(result.selected)      # one feature per cluster
    print(result.signal_rank)   # how many dimensions there really were
    ```

    Args:
        df: Frame containing the computed indicator columns.
        features: Column names to consider. Defaults to every numeric column —
            check that this excludes raw OHLCV and any target column.
        method: Correlation method, `"spearman"` (default) or `"pearson"`.
        detone: Leading eigenvalues to strip; see :func:`denoise`.
        max_k: Largest number of clusters to consider.
        k: Force an exact number of clusters; see :func:`cluster_features`.
        scores: Optional `{feature_name: score}` used to elect cluster
            representatives by `|score|` instead of by medoid; see
            :func:`select_representatives`. Names not present score 0.
        bandwidth: Kernel bandwidth for the noise-variance fit; see
            :func:`denoise`.

    Returns:
        A :class:`SelectionResult`.
    """
    screened = screen(df, features)
    corr = corr_matrix(screened.values, method=method)
    q = screened.q
    denoised = denoise(corr, q, detone=detone, bandwidth=bandwidth)

    labels, silhouette = cluster_features(denoised, max_k=max_k, k=k)
    dist = distance_matrix(denoised)
    aligned = None if scores is None else [scores.get(n, 0.0) for n in screened.names]
    selected = select_representatives(dist, labels, screened.names, aligned)

    clusters: dict[int, list[str]] = {}
    for name, label in zip(screened.names, labels):
        clusters.setdefault(int(label), []).append(name)

    return SelectionResult(
        selected=selected,
        clusters=clusters,
        labels=labels,
        names=screened.names,
        dropped=screened.dropped,
        corr=denoised,
        silhouette=silhouette,
        # Both diagnostics describe the *raw* matrix — they answer "how many
        # features did I really have?", which denoising would flatter.
        signal_rank=signal_rank(corr, q, bandwidth=bandwidth),
        effective_rank=effective_rank(corr),
        n_obs=screened.values.shape[0],
    )

block_permute

block_permute(values: ndarray, block_size: int, rng: Generator) -> ndarray

Circular block permutation, applied independently to each column.

Each column is rebuilt from randomly-placed contiguous blocks of itself, so its autocorrelation survives roughly intact while its alignment with every other column is destroyed. That is exactly the null needed for a cross-feature statistic: it asks "is the structure between features real?" without pretending the features were ever i.i.d. through time.

Parameters:

Name	Type	Description	Default
values	ndarray	Finite (n_obs, n_features) matrix.	required
block_size	int	Length of each block. Should comfortably exceed the longest indicator window in the matrix, or the permutation breaks the very autocorrelation it is meant to preserve.	required
rng	Generator	NumPy generator, so a test is reproducible.	required

Returns:

Type	Description
ndarray	A permuted matrix of the same shape.

Raises:

Type	Description
ValueError	If block_size is not in [1, n_obs].

Source code in polars_ta/selection.py

def block_permute(
    values: np.ndarray, block_size: int, rng: np.random.Generator
) -> np.ndarray:
    """Circular block permutation, applied independently to each column.

    Each column is rebuilt from randomly-placed contiguous blocks of itself, so
    its autocorrelation survives roughly intact while its *alignment* with
    every other column is destroyed. That is exactly the null needed for a
    cross-feature statistic: it asks "is the structure *between* features real?"
    without pretending the features were ever i.i.d. through time.

    Args:
        values: Finite `(n_obs, n_features)` matrix.
        block_size: Length of each block. Should comfortably exceed the longest
            indicator window in the matrix, or the permutation breaks the very
            autocorrelation it is meant to preserve.
        rng: NumPy generator, so a test is reproducible.

    Returns:
        A permuted matrix of the same shape.

    Raises:
        ValueError: If `block_size` is not in `[1, n_obs]`.
    """
    values = np.asarray(values, dtype=np.float64)
    n_obs, n_features = values.shape
    if not 1 <= block_size <= n_obs:
        raise ValueError(f"block_size must be in [1, {n_obs}], got {block_size}")

    n_blocks = int(math.ceil(n_obs / block_size))
    out = np.empty_like(values)
    for j in range(n_features):
        starts = rng.integers(0, n_obs, size=n_blocks)
        # Wrap with modulo so a block starting near the end stays full length,
        # which keeps every block the same length (a truncated tail block would
        # under-represent the series' end).
        take = (starts[:, None] + np.arange(block_size)[None, :]) % n_obs
        out[:, j] = values[:, j][take.ravel()[:n_obs]]
    return out

permutation_test

permutation_test(values: ndarray, statistic: Callable[[ndarray], float], block_size: int = 100, n_draws: int = 200, alternative: str = 'greater', seed: int | None = 0) -> PermutationTest

Block-permutation test for any statistic of a feature matrix.

test = selection.permutation_test(
    screened.values,
    lambda v: selection.signal_rank(selection.corr_matrix(v), q),
    block_size=250,
)
print(test.observed, test.null_mean, test.p_value)

Note that the statistic is recomputed on every draw, so a 200-draw test costs 200 correlation matrices and eigendecompositions. That is seconds for a few dozen features and minutes for a few hundred.

Parameters:

Name	Type	Description	Default
values	ndarray	Finite (n_obs, n_features) matrix.	required
statistic	Callable[[ndarray], float]	Callable mapping a matrix to a scalar — e.g. a closure over :func:signal_rank or :func:effective_rank.	required
block_size	int	Block length for :func:block_permute.	100
n_draws	int	Number of permuted draws.	200
alternative	str	"greater" when real structure makes the statistic rise (signal_rank), "less" when it makes it fall (effective_rank, which is maximal for independent features).	'greater'
seed	int | None	Seed for the generator; None for a non-reproducible run.	0

Returns:

Name	Type	Description
A	PermutationTest	class:PermutationTest.

Raises:

Type	Description
ValueError	If n_draws < 1 or alternative is not "greater" or "less".

Source code in polars_ta/selection.py

def permutation_test(
    values: np.ndarray,
    statistic: Callable[[np.ndarray], float],
    block_size: int = 100,
    n_draws: int = 200,
    alternative: str = "greater",
    seed: int | None = 0,
) -> PermutationTest:
    """Block-permutation test for any statistic of a feature matrix.

    ```python
    test = selection.permutation_test(
        screened.values,
        lambda v: selection.signal_rank(selection.corr_matrix(v), q),
        block_size=250,
    )
    print(test.observed, test.null_mean, test.p_value)
    ```

    Note that the statistic is recomputed on every draw, so a 200-draw test
    costs 200 correlation matrices and eigendecompositions. That is seconds for
    a few dozen features and minutes for a few hundred.

    Args:
        values: Finite `(n_obs, n_features)` matrix.
        statistic: Callable mapping a matrix to a scalar — e.g. a closure over
            :func:`signal_rank` or :func:`effective_rank`.
        block_size: Block length for :func:`block_permute`.
        n_draws: Number of permuted draws.
        alternative: `"greater"` when real structure makes the statistic *rise*
            (`signal_rank`), `"less"` when it makes it *fall*
            (`effective_rank`, which is maximal for independent features).
        seed: Seed for the generator; `None` for a non-reproducible run.

    Returns:
        A :class:`PermutationTest`.

    Raises:
        ValueError: If `n_draws < 1` or `alternative` is not `"greater"` or
            `"less"`.
    """
    if n_draws < 1:
        raise ValueError(f"n_draws must be at least 1, got {n_draws}")
    if alternative not in ("greater", "less"):
        raise ValueError(
            f"alternative must be 'greater' or 'less', got {alternative!r}"
        )

    rng = np.random.default_rng(seed)
    observed = float(statistic(np.asarray(values, dtype=np.float64)))
    null = np.array(
        [
            float(statistic(block_permute(values, block_size, rng)))
            for _ in range(n_draws)
        ]
    )
    if alternative == "greater":
        at_least_as_extreme = int((null >= observed).sum())
    else:
        at_least_as_extreme = int((null <= observed).sum())
    return PermutationTest(
        observed=observed,
        null=null,
        p_value=(1 + at_least_as_extreme) / (1 + n_draws),
        alternative=alternative,
    )

purged_kfold

purged_kfold(n_samples: int, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0) -> list[tuple[ndarray, ndarray]]

Contiguous K-fold splits with purging and an embargo.

Two corrections to an ordinary K-fold, both from López de Prado (2018, ch. 7), and both mandatory on overlapping financial features:

• Purging removes training samples whose label window overlaps the test fold. With a horizon-h forward return, sample i peeks at bars i+1 … i+h; if any of those fall in the test fold, training on i is training on the answer.
• Embargo additionally drops the embargo samples immediately after the test fold. Purging handles the label overlap; the embargo handles serial correlation in the features themselves, which leaks in the same direction. Set it to at least the longest indicator window.

Folds are contiguous blocks, never shuffled — shuffling a time series destroys the ordering that makes purging meaningful in the first place.

cv = selection.purged_kfold(
    n_samples=len(y), n_splits=5, label_horizon=12, embargo=100,
)

Parameters:

Name	Type	Description	Default
n_samples	int	Number of rows, in time order.	required
n_splits	int	Number of folds.	5
label_horizon	int	How many bars ahead each label looks.	1
embargo	int	Extra samples to drop after each test fold.	0

Returns:

Type	Description
list[tuple[ndarray, ndarray]]	A list of (train_indices, test_indices) pairs.

Raises:

Type	Description
ValueError	If n_splits is outside [2, n_samples], if label_horizon or embargo is negative, or if the purge and embargo leave a fold with no training data at all.

Source code in polars_ta/selection.py

def purged_kfold(
    n_samples: int, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0
) -> list[tuple[np.ndarray, np.ndarray]]:
    """Contiguous K-fold splits with purging and an embargo.

    Two corrections to an ordinary K-fold, both from López de Prado (2018,
    ch. 7), and both mandatory on overlapping financial features:

    - **Purging** removes training samples whose *label* window overlaps the
      test fold. With a horizon-`h` forward return, sample `i` peeks at bars
      `i+1 … i+h`; if any of those fall in the test fold, training on `i` is
      training on the answer.
    - **Embargo** additionally drops the `embargo` samples immediately *after*
      the test fold. Purging handles the label overlap; the embargo handles
      serial correlation in the features themselves, which leaks in the same
      direction. Set it to at least the longest indicator window.

    Folds are contiguous blocks, never shuffled — shuffling a time series
    destroys the ordering that makes purging meaningful in the first place.

    ```python
    cv = selection.purged_kfold(
        n_samples=len(y), n_splits=5, label_horizon=12, embargo=100,
    )
    ```

    Args:
        n_samples: Number of rows, in time order.
        n_splits: Number of folds.
        label_horizon: How many bars ahead each label looks.
        embargo: Extra samples to drop after each test fold.

    Returns:
        A list of `(train_indices, test_indices)` pairs.

    Raises:
        ValueError: If `n_splits` is outside `[2, n_samples]`, if
            `label_horizon` or `embargo` is negative, or if the purge and
            embargo leave a fold with no training data at all.
    """
    if not 2 <= n_splits <= n_samples:
        raise ValueError(f"n_splits must be in [2, {n_samples}], got {n_splits}")
    if label_horizon < 0:
        raise ValueError(f"label_horizon must be >= 0, got {label_horizon}")
    if embargo < 0:
        raise ValueError(f"embargo must be >= 0, got {embargo}")

    indices = np.arange(n_samples)
    bounds = np.linspace(0, n_samples, n_splits + 1).astype(int)

    splits = []
    for fold in range(n_splits):
        start, stop = bounds[fold], bounds[fold + 1]
        test = indices[start:stop]
        # Everything from `label_horizon` before the fold to `embargo` after it
        # is contaminated; the rest is safe to train on.
        blocked = (indices >= start - label_horizon) & (indices < stop + embargo)
        train = indices[~blocked]
        if train.size == 0:
            raise ValueError(
                f"fold {fold} has no training data left after purging "
                f"({label_horizon}) and embargo ({embargo}); reduce them or "
                f"use fewer splits"
            )
        splits.append((train, test))
    return splits

clustered_mda

clustered_mda(df: DataFrame, features: Sequence[str], target: str, labels: ndarray | None = None, model_factory: Callable[[], object] | None = None, scorer: Callable[[ndarray, ndarray], float] | None = None, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0, seed: int | None = 0) -> ImportanceResult

Clustered Mean-Decrease-Accuracy importance under purged cross-validation.

Ordinary MDA shuffles one feature at a time and measures the drop in out-of-sample score. On correlated features that is close to meaningless: shuffle rsi_14 and the model simply leans on rsi_21, so both score zero. Clustered MDA shuffles every member of a cluster jointly, with the same row permutation, so the within-cluster structure survives while the cluster's relationship to the target is destroyed. What it measures is then attributable to the cluster as a whole (López de Prado, 2020, ch. 6).

The cross-validation is :func:purged_kfold, not a plain K-fold — see that function for why the difference is not optional.

scored = feats.with_columns(
    fwd=pl.col("close").pct_change().shift(-12),
)
sel = selection.select_features(scored, FEATURES, k=4)
imp = selection.clustered_mda(
    scored, sel.names, "fwd", labels=sel.labels,
    label_horizon=12, embargo=100,
)
for cluster in imp.ranked:
    print(imp.clusters[cluster], round(imp.importance[cluster], 4))

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column — typically a forward return, which is look-ahead by construction and therefore research-only.	required
labels	ndarray | None	Cluster id per feature, e.g. SelectionResult.labels. Defaults to one cluster per feature, which reintroduces the substitution effect and is only sensible on an already-decorrelated set.	None
model_factory	Callable[[], object] | None	Zero-argument callable returning a fresh object with fit(X, y) and predict(X). Defaults to :class:RidgeRegressor.	None
scorer	Callable[[ndarray, ndarray], float] | None	(y_true, y_pred) -> float, higher is better. Defaults to the out-of-sample information coefficient.	None
n_splits	int	Cross-validation folds.	5
label_horizon	int	Bars the target looks ahead — must match how target was built, or purging removes the wrong rows.	1
embargo	int	Extra samples dropped after each test fold; set it to at least the longest indicator window in features.	0
seed	int | None	Seed for the shuffles.	0

Returns:

Name	Type	Description
An	ImportanceResult	class:ImportanceResult.

Raises:

Type	Description
ValueError	If fewer than two features are given, if labels does not align with features, or if too few rows survive alignment.

Source code in polars_ta/selection.py

def clustered_mda(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    labels: np.ndarray | None = None,
    model_factory: Callable[[], object] | None = None,
    scorer: Callable[[np.ndarray, np.ndarray], float] | None = None,
    n_splits: int = 5,
    label_horizon: int = 1,
    embargo: int = 0,
    seed: int | None = 0,
) -> ImportanceResult:
    """Clustered Mean-Decrease-Accuracy importance under purged cross-validation.

    Ordinary MDA shuffles one feature at a time and measures the drop in
    out-of-sample score. On correlated features that is close to meaningless:
    shuffle `rsi_14` and the model simply leans on `rsi_21`, so both score
    zero. **Clustered** MDA shuffles every member of a cluster *jointly*, with
    the same row permutation, so the within-cluster structure survives while
    the cluster's relationship to the target is destroyed. What it measures is
    then attributable to the cluster as a whole (López de Prado, 2020, ch. 6).

    The cross-validation is :func:`purged_kfold`, not a plain K-fold — see that
    function for why the difference is not optional.

    ```python
    scored = feats.with_columns(
        fwd=pl.col("close").pct_change().shift(-12),
    )
    sel = selection.select_features(scored, FEATURES, k=4)
    imp = selection.clustered_mda(
        scored, sel.names, "fwd", labels=sel.labels,
        label_horizon=12, embargo=100,
    )
    for cluster in imp.ranked:
        print(imp.clusters[cluster], round(imp.importance[cluster], 4))
    ```

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column — typically a *forward* return, which is
            look-ahead by construction and therefore research-only.
        labels: Cluster id per feature, e.g. `SelectionResult.labels`. Defaults
            to one cluster per feature, which reintroduces the substitution
            effect and is only sensible on an already-decorrelated set.
        model_factory: Zero-argument callable returning a fresh object with
            `fit(X, y)` and `predict(X)`. Defaults to :class:`RidgeRegressor`.
        scorer: `(y_true, y_pred) -> float`, higher is better. Defaults to the
            out-of-sample information coefficient.
        n_splits: Cross-validation folds.
        label_horizon: Bars the target looks ahead — **must** match how
            `target` was built, or purging removes the wrong rows.
        embargo: Extra samples dropped after each test fold; set it to at least
            the longest indicator window in `features`.
        seed: Seed for the shuffles.

    Returns:
        An :class:`ImportanceResult`.

    Raises:
        ValueError: If fewer than two features are given, if `labels` does not
            align with `features`, or if too few rows survive alignment.
    """
    features, labels = _check_features_and_labels(features, labels)
    x, y = _aligned_xy(df, features, target, n_splits)

    model_factory = RidgeRegressor if model_factory is None else model_factory
    scorer = _information_coefficient if scorer is None else scorer
    rng = np.random.default_rng(seed)
    splits = purged_kfold(x.shape[0], n_splits, label_horizon, embargo)

    unique = np.unique(labels)
    drops: dict[int, list[float]] = {int(c): [] for c in unique}
    baselines = []
    for train, test in splits:
        model = model_factory()
        model.fit(x[train], y[train])
        baseline = scorer(y[test], model.predict(x[test]))
        baselines.append(baseline)
        for cluster in unique:
            columns = np.flatnonzero(labels == cluster)
            shuffled = x[test].copy()
            # One permutation for the whole cluster: this is what separates
            # clustered MDA from shuffling each column independently.
            order = rng.permutation(shuffled.shape[0])
            shuffled[:, columns] = shuffled[np.ix_(order, columns)]
            drops[int(cluster)].append(
                baseline - scorer(y[test], model.predict(shuffled))
            )

    return ImportanceResult(
        clusters=_clusters_of(features, labels),
        importance={c: float(np.mean(v)) for c, v in drops.items()},
        std={c: float(np.std(v)) for c, v in drops.items()},
        baseline=float(np.mean(baselines)),
        n_splits=len(splits),
    )

target_screen

target_screen(df: DataFrame, features: Sequence[str], target: str, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0, bins: int | None = None) -> DataFrame

Model-free univariate screen of every feature against the target.

The cheapest question worth asking, and the one to ask first: before fitting anything, does this feature co-move with the label at all? Two complementary statistics, both computed out of sample on the test folds of a :func:purged_kfold, so a feature cannot score well by being memorized:

• Information coefficient — the Spearman rank correlation between the feature and the target, the standard unit of predictive skill on a research desk (Grinold & Kahn, 1999). Rank rather than Pearson because a single outlier bar otherwise sets the number, and because most indicators are monotone transforms of each other anyway. On real forward returns a mean |IC| of 0.02-0.05 is a genuinely useful feature; anything above 0.1 should be treated as a bug report against your label until proven otherwise.
• Mutual information — normalized, from the same histogram estimator as :func:mutual_information. It catches what the IC structurally cannot: a U-shaped relationship (volatility features often have one — both tails predict, the middle does not) scores an IC of zero and an MI well above it. A feature with high MI and near-zero IC is not useless, it is non-monotone, and needs a model that can represent that.

The per-fold standard deviation matters as much as the mean. A feature with IC 0.04 ± 0.01 is a feature; IC 0.04 ± 0.09 is one lucky fold, and the ic_t_stat column is there to make that distinction impossible to miss.

ic_t_stat measures consistency, not size

Read it next to abs_ic, never instead of it. A feature whose IC is negligible but reliably signed gets a large t-stat — in this library's own tests a pure-noise feature scores |IC| 0.008 with a t-stat near 5, because 0.008 is tiny against the label but not against its own fold-to-fold spread. That is a stable nothing. Sort by abs_ic, then use the t-stat to discard whatever near the top is not reproducible.

A screen, not a selector

Univariate scores are blind to interaction — a feature that only predicts conditional on the regime scores zero here — and they are blind to redundancy, so the top 10 will typically be 10 readings of one quantity. Use this to cut the obvious dead weight, then let :func:clustered_mda do the actual ranking.

screen = selection.target_screen(
    scored, FEATURES, "fwd", label_horizon=12, embargo=100,
)
screen.head(10)                              # the most promising features
screen.filter(pl.col("ic_t_stat").abs() > 2) # ...that are also stable

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column, typically a forward return.	required
n_splits	int	Cross-validation folds. Scores are averaged over the folds' test sets.	5
label_horizon	int	Bars the target looks ahead; see :func:purged_kfold.	1
embargo	int	Extra samples dropped after each test fold.	0
bins	int | None	Histogram bins for the mutual information; defaults to the Hacine-Gharbi optimum for the fold size.	None

Returns:

Type	Description
DataFrame	One row per feature, sorted by descending abs_ic, with columns
DataFrame	feature, ic (mean rank correlation across folds), ic_std,
DataFrame	abs_ic, ic_t_stat (ic / (ic_std / sqrt(n_splits))) and
DataFrame	mutual_info.

Raises:

Type	Description
ValueError	If fewer than two features are given or too few rows survive alignment.

Source code in polars_ta/selection.py

def target_screen(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    n_splits: int = 5,
    label_horizon: int = 1,
    embargo: int = 0,
    bins: int | None = None,
) -> pl.DataFrame:
    """Model-free univariate screen of every feature against the target.

    The cheapest question worth asking, and the one to ask *first*: before
    fitting anything, does this feature co-move with the label at all? Two
    complementary statistics, both computed **out of sample** on the test folds
    of a :func:`purged_kfold`, so a feature cannot score well by being
    memorized:

    - **Information coefficient** — the Spearman rank correlation between the
      feature and the target, the standard unit of predictive skill on a
      research desk (Grinold & Kahn, 1999). Rank rather than Pearson because a
      single outlier bar otherwise sets the number, and because most indicators
      are monotone transforms of each other anyway. On real forward returns a
      *mean* |IC| of 0.02-0.05 is a genuinely useful feature; anything above 0.1
      should be treated as a bug report against your label until proven
      otherwise.
    - **Mutual information** — normalized, from the same histogram estimator as
      :func:`mutual_information`. It catches what the IC structurally cannot: a
      U-shaped relationship (volatility features often have one — both tails
      predict, the middle does not) scores an IC of zero and an MI well above
      it. A feature with high MI and near-zero IC is not useless, it is
      *non-monotone*, and needs a model that can represent that.

    The per-fold standard deviation matters as much as the mean. A feature with
    IC 0.04 ± 0.01 is a feature; IC 0.04 ± 0.09 is one lucky fold, and the
    `ic_t_stat` column is there to make that distinction impossible to miss.

    !!! warning "`ic_t_stat` measures consistency, not size"
        Read it *next to* `abs_ic`, never instead of it. A feature whose IC is
        negligible but reliably signed gets a large t-stat — in this library's
        own tests a pure-noise feature scores |IC| 0.008 with a t-stat near 5,
        because 0.008 is tiny against the label but not against its own
        fold-to-fold spread. That is a stable nothing. Sort by `abs_ic`, then
        use the t-stat to discard whatever near the top is not reproducible.

    !!! warning "A screen, not a selector"
        Univariate scores are blind to interaction — a feature that only
        predicts *conditional on the regime* scores zero here — and they are
        blind to redundancy, so the top 10 will typically be 10 readings of one
        quantity. Use this to cut the obvious dead weight, then let
        :func:`clustered_mda` do the actual ranking.

    ```python
    screen = selection.target_screen(
        scored, FEATURES, "fwd", label_horizon=12, embargo=100,
    )
    screen.head(10)                              # the most promising features
    screen.filter(pl.col("ic_t_stat").abs() > 2) # ...that are also stable
    ```

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column, typically a forward return.
        n_splits: Cross-validation folds. Scores are averaged over the folds'
            test sets.
        label_horizon: Bars the target looks ahead; see :func:`purged_kfold`.
        embargo: Extra samples dropped after each test fold.
        bins: Histogram bins for the mutual information; defaults to the
            Hacine-Gharbi optimum for the fold size.

    Returns:
        One row per feature, sorted by descending `abs_ic`, with columns
        `feature`, `ic` (mean rank correlation across folds), `ic_std`,
        `abs_ic`, `ic_t_stat` (`ic / (ic_std / sqrt(n_splits))`) and
        `mutual_info`.

    Raises:
        ValueError: If fewer than two features are given or too few rows survive
            alignment.
    """
    features, _ = _check_features_and_labels(features, None)
    x, y = _aligned_xy(df, features, target, n_splits)
    splits = purged_kfold(x.shape[0], n_splits, label_horizon, embargo)

    rows = []
    for j, name in enumerate(features):
        ics, infos = [], []
        for _, test in splits:
            feature, label = x[test, j], y[test]
            if feature.std() == 0.0 or label.std() == 0.0:
                ics.append(0.0)
                infos.append(0.0)
                continue
            ics.append(
                float(np.corrcoef(_average_ranks(feature), _average_ranks(label))[0, 1])
            )
            n_bins = _bin_count(test.size) if bins is None else bins
            h_x, h_y, h_xy = _entropies(feature, label, n_bins)
            floor = min(h_x, h_y)
            infos.append(0.0 if floor <= 0.0 else (h_x + h_y - h_xy) / floor)

        mean_ic = float(np.mean(ics))
        std_ic = float(np.std(ics))
        rows.append(
            {
                "feature": name,
                "ic": mean_ic,
                "ic_std": std_ic,
                "abs_ic": abs(mean_ic),
                "ic_t_stat": (
                    0.0
                    if std_ic == 0.0
                    else mean_ic / (std_ic / math.sqrt(len(splits)))
                ),
                "mutual_info": float(np.mean(infos)),
            }
        )
    return pl.DataFrame(rows).sort("abs_ic", descending=True)

clustered_mdi

clustered_mdi(df: DataFrame, features: Sequence[str], target: str, model_factory: Callable[[], object], labels: ndarray | None = None, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0) -> ImportanceResult

Clustered Mean-Decrease-Impurity importance (López de Prado, 2020, ch. 6).

MDI is what a tree ensemble reports for free after fitting: for each feature, the total impurity reduction across every node that split on it, normalized to sum to 1. Clustered MDI sums those shares within each cluster, which repairs MDI's worst failure mode — with rsi_14 and rsi_21 both available, each tree picks one arbitrarily and the pair's shared importance is split into two unimpressive halves.

Read it as what the model used, and be clear about what that is not:

• It is in-sample. MDI is computed on the training data, so a feature the model overfit to scores high. This is the essential difference from :func:clustered_mda, which measures out-of-sample degradation, and it is why the two disagreeing is diagnostic: high MDI with zero MDA means the model memorized that feature.
• It is biased toward high-cardinality features. A continuous feature offers a tree far more places to split than a 0/±100 candlestick pattern does, and MDI rewards it for that regardless of predictive content. On a mixed feature set — which this library's is — do not compare a continuous feature's MDI against a discrete one's.
• Every feature gets a non-zero score. MDI has no null: a pure-noise feature still gets used somewhere and still scores above zero. Only MDA and :func:null_importance can say "this is worth nothing".

Measured on the full library, MDI ranked it backwards

Over 183 indicators in 8 clusters against a 12-bar forward return, a RandomForestRegressor spent 47% of its impurity budget on the cluster whose out-of-sample MDA was 0.0001 and whose Shapley value was −0.021, and 0.001 on the three candlestick clusters — one of which had the best MDA and the lowest p-value in the study. Both failure modes above, in one table. Use this to detect memorization by disagreeing with :func:clustered_mda, never to pick features on its own.

In exchange it is essentially free (no refitting, no shuffling) and it is the only one of the three that sees interactions directly, since a split deep in a tree is conditional on every split above it.

Requires a model exposing feature_importances_ after fit — any scikit-learn tree ensemble, LightGBM or XGBoost. There is no default: :class:RidgeRegressor has no impurity to decrease, and silently substituting coefficients would be a different measure wearing this one's name.

from sklearn.ensemble import RandomForestRegressor

mdi = selection.clustered_mdi(
    scored, sel.names, "fwd",
    model_factory=lambda: RandomForestRegressor(n_estimators=200, max_depth=4),
    labels=sel.labels, label_horizon=12, embargo=100,
)

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column.	required
model_factory	Callable[[], object]	Zero-argument callable returning a fresh model whose fitted form exposes feature_importances_. Required.	required
labels	ndarray | None	Cluster id per feature, e.g. SelectionResult.labels. Defaults to one cluster per feature, i.e. plain per-feature MDI.	None
n_splits	int	Folds. Each fold fits on its (purged) training set and contributes one set of importances, so std measures how stable the model's choice of feature is across the sample.	5
label_horizon	int	Bars the target looks ahead; see :func:purged_kfold.	1
embargo	int	Extra samples dropped after each test fold.	0

Returns:

Name	Type	Description
An	ImportanceResult	class:ImportanceResult whose importance values are impurity
	ImportanceResult	shares summing to 1 across clusters, and whose baseline is the mean
	ImportanceResult	in-sample score of the fitted models — a number worth glancing at, since
	ImportanceResult	a baseline far above what :func:clustered_mda reports out of sample is
	ImportanceResult	overfitting stated in one figure.

Raises:

Type	Description
ValueError	If fewer than two features are given, if labels does not align with features, if too few rows survive alignment, or if the fitted model exposes no feature_importances_.

Source code in polars_ta/selection.py

def clustered_mdi(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    model_factory: Callable[[], object],
    labels: np.ndarray | None = None,
    n_splits: int = 5,
    label_horizon: int = 1,
    embargo: int = 0,
) -> ImportanceResult:
    """Clustered Mean-Decrease-Impurity importance (López de Prado, 2020, ch. 6).

    MDI is what a tree ensemble reports for free after fitting: for each
    feature, the total impurity reduction across every node that split on it,
    normalized to sum to 1. **Clustered** MDI sums those shares within each
    cluster, which repairs MDI's worst failure mode — with `rsi_14` and
    `rsi_21` both available, each tree picks one arbitrarily and the pair's
    shared importance is split into two unimpressive halves.

    Read it as *what the model used*, and be clear about what that is not:

    - **It is in-sample.** MDI is computed on the training data, so a feature
      the model overfit to scores high. This is the essential difference from
      :func:`clustered_mda`, which measures out-of-sample degradation, and it is
      why the two disagreeing is diagnostic: high MDI with zero MDA means the
      model memorized that feature.
    - **It is biased toward high-cardinality features.** A continuous feature
      offers a tree far more places to split than a 0/±100 candlestick pattern
      does, and MDI rewards it for that regardless of predictive content. On a
      mixed feature set — which this library's is — do not compare a continuous
      feature's MDI against a discrete one's.
    - **Every feature gets a non-zero score.** MDI has no null: a pure-noise
      feature still gets used somewhere and still scores above zero. Only MDA
      and :func:`null_importance` can say "this is worth nothing".

    !!! danger "Measured on the full library, MDI ranked it backwards"
        Over 183 indicators in 8 clusters against a 12-bar forward return, a
        `RandomForestRegressor` spent **47%** of its impurity budget on the
        cluster whose out-of-sample MDA was 0.0001 and whose Shapley value was
        −0.021, and **0.001** on the three candlestick clusters — one of which
        had the best MDA and the lowest p-value in the study. Both failure modes
        above, in one table. Use this to *detect* memorization by disagreeing
        with :func:`clustered_mda`, never to pick features on its own.

    In exchange it is essentially free (no refitting, no shuffling) and it is
    the only one of the three that sees *interactions* directly, since a split
    deep in a tree is conditional on every split above it.

    Requires a model exposing `feature_importances_` after `fit` — any
    scikit-learn tree ensemble, LightGBM or XGBoost. There is no default:
    :class:`RidgeRegressor` has no impurity to decrease, and silently
    substituting coefficients would be a different measure wearing this one's
    name.

    ```python
    from sklearn.ensemble import RandomForestRegressor

    mdi = selection.clustered_mdi(
        scored, sel.names, "fwd",
        model_factory=lambda: RandomForestRegressor(n_estimators=200, max_depth=4),
        labels=sel.labels, label_horizon=12, embargo=100,
    )
    ```

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column.
        model_factory: Zero-argument callable returning a fresh model whose
            fitted form exposes `feature_importances_`. Required.
        labels: Cluster id per feature, e.g. `SelectionResult.labels`. Defaults
            to one cluster per feature, i.e. plain per-feature MDI.
        n_splits: Folds. Each fold fits on its (purged) training set and
            contributes one set of importances, so `std` measures how stable the
            model's choice of feature is across the sample.
        label_horizon: Bars the target looks ahead; see :func:`purged_kfold`.
        embargo: Extra samples dropped after each test fold.

    Returns:
        An :class:`ImportanceResult` whose `importance` values are impurity
        shares summing to 1 across clusters, and whose `baseline` is the mean
        in-sample score of the fitted models — a number worth glancing at, since
        a baseline far above what :func:`clustered_mda` reports out of sample is
        overfitting stated in one figure.

    Raises:
        ValueError: If fewer than two features are given, if `labels` does not
            align with `features`, if too few rows survive alignment, or if the
            fitted model exposes no `feature_importances_`.
    """
    features, labels = _check_features_and_labels(features, labels)
    x, y = _aligned_xy(df, features, target, n_splits)
    splits = purged_kfold(x.shape[0], n_splits, label_horizon, embargo)

    unique = np.unique(labels)
    shares: dict[int, list[float]] = {int(c): [] for c in unique}
    baselines = []
    for train, _ in splits:
        model = model_factory()
        model.fit(x[train], y[train])
        raw = getattr(model, "feature_importances_", None)
        if raw is None:
            raise ValueError(
                f"{type(model).__name__} exposes no feature_importances_ after "
                "fit — MDI needs a tree-based model (RandomForest, "
                "GradientBoosting, LightGBM, XGBoost). Use clustered_mda for "
                "models that do not report impurity."
            )
        raw = np.asarray(raw, dtype=np.float64)
        if raw.shape != (len(features),):
            raise ValueError(
                f"feature_importances_ has shape {raw.shape}, expected "
                f"({len(features)},)"
            )
        # Normalize per fold, so a fold whose model reports unnormalized gains
        # cannot dominate the average purely by scale.
        total = raw.sum()
        raw = raw / total if total > 0.0 else np.full_like(raw, 1.0 / raw.size)
        for cluster in unique:
            shares[int(cluster)].append(float(raw[labels == cluster].sum()))
        baselines.append(_information_coefficient(y[train], model.predict(x[train])))

    return ImportanceResult(
        clusters=_clusters_of(features, labels),
        importance={c: float(np.mean(v)) for c, v in shares.items()},
        std={c: float(np.std(v)) for c, v in shares.items()},
        baseline=float(np.mean(baselines)),
        n_splits=len(splits),
    )

single_feature_importance

single_feature_importance(df: DataFrame, features: Sequence[str], target: str, labels: ndarray | None = None, model_factory: Callable[[], object] | None = None, scorer: Callable[[ndarray, ndarray], float] | None = None, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0) -> ImportanceResult

Single-feature importance: fit each cluster alone (SFI).

The third leg of López de Prado's importance set (2018, ch. 8). MDI and MDA both score a feature in the presence of all the others, which is what makes them vulnerable to substitution — and what :func:clustered_mda's clustering only partly repairs, since features can substitute across cluster boundaries too. SFI sidesteps the problem structurally: fit a model on one cluster's columns and nothing else, and score it out of sample. There is nothing left to substitute with, so the number is unambiguous.

The price is exactly symmetric: SFI is blind to everything joint. A feature that is worthless alone but unlocks another one — a volatility measure that tells you when to trust a momentum signal is the canonical case — scores zero here and high under MDA. That disagreement is the point of computing both:

SFI	MDA	Reading
high	high	Genuinely predictive. Keep.
high	low	Redundant — something else already carries it.
low	high	Only works in combination. Keep, but never alone.
low	low	Nothing there.

Because each fit sees one cluster, SFI is also the only one of the three that is immune to the curse of dimensionality in the fit itself: a model on 200 collinear features is badly conditioned no matter how it is scored, and a model on 3 is not.

sfi = selection.single_feature_importance(
    scored, sel.names, "fwd",
    labels=sel.labels, label_horizon=12, embargo=100,
)
mda = selection.clustered_mda(scored, sel.names, "fwd", labels=sel.labels)
for c in sfi.ranked:
    print(sfi.clusters[c], round(sfi.importance[c], 4), round(mda.importance[c], 4))

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column.	required
labels	ndarray | None	Cluster id per feature. Defaults to one cluster per feature, which is classic per-feature SFI; pass SelectionResult.labels to score whole clusters together instead.	None
model_factory	Callable[[], object] | None	Zero-argument callable returning a fresh object with fit(X, y) and predict(X). Defaults to :class:RidgeRegressor.	None
scorer	Callable[[ndarray, ndarray], float] | None	(y_true, y_pred) -> float, higher is better. Defaults to the out-of-sample information coefficient.	None
n_splits	int	Cross-validation folds.	5
label_horizon	int	Bars the target looks ahead; see :func:purged_kfold.	1
embargo	int	Extra samples dropped after each test fold.	0

Returns:

Name	Type	Description
An	ImportanceResult	class:ImportanceResult whose importance is each cluster's mean
	ImportanceResult	out-of-sample score on its own — directly comparable to baseline,
	ImportanceResult	which is the score of the model fitted on all features. A cluster
	ImportanceResult	scoring near the baseline by itself is doing all the work; every cluster
	ImportanceResult	scoring far below it means the edge is genuinely joint.

Raises:

Type	Description
ValueError	If fewer than two features are given, if labels does not align with features, or if too few rows survive alignment.

Source code in polars_ta/selection.py

def single_feature_importance(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    labels: np.ndarray | None = None,
    model_factory: Callable[[], object] | None = None,
    scorer: Callable[[np.ndarray, np.ndarray], float] | None = None,
    n_splits: int = 5,
    label_horizon: int = 1,
    embargo: int = 0,
) -> ImportanceResult:
    """Single-feature importance: fit each cluster **alone** (SFI).

    The third leg of López de Prado's importance set (2018, ch. 8). MDI and MDA
    both score a feature *in the presence of all the others*, which is what
    makes them vulnerable to substitution — and what
    :func:`clustered_mda`'s clustering only partly repairs, since features can
    substitute across cluster boundaries too. SFI sidesteps the problem
    structurally: fit a model on one cluster's columns and nothing else, and
    score it out of sample. There is nothing left to substitute with, so the
    number is unambiguous.

    The price is exactly symmetric: SFI is blind to everything joint. A feature
    that is worthless alone but unlocks another one — a volatility measure that
    tells you *when* to trust a momentum signal is the canonical case — scores
    zero here and high under MDA. That disagreement is the point of computing
    both:

    | SFI | MDA | Reading |
    |---|---|---|
    | high | high | Genuinely predictive. Keep. |
    | high | low | Redundant — something else already carries it. |
    | low | high | Only works in combination. Keep, but never alone. |
    | low | low | Nothing there. |

    Because each fit sees one cluster, SFI is also the only one of the three
    that is immune to the *curse of dimensionality* in the fit itself: a model
    on 200 collinear features is badly conditioned no matter how it is scored,
    and a model on 3 is not.

    ```python
    sfi = selection.single_feature_importance(
        scored, sel.names, "fwd",
        labels=sel.labels, label_horizon=12, embargo=100,
    )
    mda = selection.clustered_mda(scored, sel.names, "fwd", labels=sel.labels)
    for c in sfi.ranked:
        print(sfi.clusters[c], round(sfi.importance[c], 4), round(mda.importance[c], 4))
    ```

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column.
        labels: Cluster id per feature. Defaults to one cluster per feature,
            which is *classic* per-feature SFI; pass `SelectionResult.labels` to
            score whole clusters together instead.
        model_factory: Zero-argument callable returning a fresh object with
            `fit(X, y)` and `predict(X)`. Defaults to :class:`RidgeRegressor`.
        scorer: `(y_true, y_pred) -> float`, higher is better. Defaults to the
            out-of-sample information coefficient.
        n_splits: Cross-validation folds.
        label_horizon: Bars the target looks ahead; see :func:`purged_kfold`.
        embargo: Extra samples dropped after each test fold.

    Returns:
        An :class:`ImportanceResult` whose `importance` is each cluster's mean
        out-of-sample score *on its own* — directly comparable to `baseline`,
        which is the score of the model fitted on **all** features. A cluster
        scoring near the baseline by itself is doing all the work; every cluster
        scoring far below it means the edge is genuinely joint.

    Raises:
        ValueError: If fewer than two features are given, if `labels` does not
            align with `features`, or if too few rows survive alignment.
    """
    features, labels = _check_features_and_labels(features, labels)
    x, y = _aligned_xy(df, features, target, n_splits)

    model_factory = RidgeRegressor if model_factory is None else model_factory
    scorer = _information_coefficient if scorer is None else scorer
    splits = purged_kfold(x.shape[0], n_splits, label_horizon, embargo)

    unique = np.unique(labels)
    scores: dict[int, list[float]] = {int(c): [] for c in unique}
    baselines = []
    for train, test in splits:
        full = model_factory()
        full.fit(x[train], y[train])
        baselines.append(scorer(y[test], full.predict(x[test])))
        for cluster in unique:
            columns = np.flatnonzero(labels == cluster)
            model = model_factory()
            model.fit(x[np.ix_(train, columns)], y[train])
            scores[int(cluster)].append(
                scorer(y[test], model.predict(x[np.ix_(test, columns)]))
            )

    return ImportanceResult(
        clusters=_clusters_of(features, labels),
        importance={c: float(np.mean(v)) for c, v in scores.items()},
        std={c: float(np.std(v)) for c, v in scores.items()},
        baseline=float(np.mean(baselines)),
        n_splits=len(splits),
    )

clustered_shapley

clustered_shapley(df: DataFrame, features: Sequence[str], target: str, labels: ndarray | None = None, model_factory: Callable[[], object] | None = None, scorer: Callable[[ndarray, ndarray], float] | None = None, n_splits: int = 5, label_horizon: int = 1, embargo: int = 0, max_clusters: int = 12) -> ImportanceResult

Exact Shapley attribution of out-of-sample score across clusters.

MDA and SFI each answer half the question. MDA asks "what breaks if I remove this last?", SFI asks "what does this earn first?", and on interacting features those give different — sometimes wildly different — answers, with no principled way to reconcile them. The Shapley value (Shapley, 1953) is that reconciliation: it averages a cluster's marginal contribution over every possible order in which the clusters could have been added,

\[\phi_i = \sum_{S \subseteq N \setminus \{i\}} \frac{|S|!\,(|N|-|S|-1)!}{|N|!}\,\bigl[v(S \cup \{i\}) - v(S)\bigr]\]

where \(v(S)\) is the out-of-sample score of a model fitted on the clusters in \(S\). SFI is the \(S = \emptyset\) term and MDA is close to the \(S = N \setminus \{i\}\) term; Shapley is the properly weighted average of those and everything in between.

What this buys you is the property neither of the others has — efficiency: the values sum exactly to \(v(N) - v(\emptyset)\), the full model's score over an empty one. So "this cluster is 40% of the edge" is a statement you can defend, and the values can be compared, added and budgeted. It is also the only one of the four that can report a negative importance honestly: a cluster that consistently makes the model worse gets a negative \(\phi\), not a floor of zero.

Cost is exponential in the number of clusters

Every subset must be fitted: \(2^k\) models per fold. That is 32 for k=5 and 4096 for k=12, times n_splits. max_clusters guards against an accidental 200-cluster call that would never return; raise it only if you know what you are asking for. This is why the function takes clusters rather than features — running it per-feature on a wide matrix is not merely slow, it is infeasible, and cluster-level attribution is the more meaningful question anyway.

shap = selection.clustered_shapley(
    scored, sel.names, "fwd",
    labels=sel.labels, label_horizon=12, embargo=100,
)
total = sum(shap.importance.values())
for c in shap.ranked:
    print(shap.clusters[c], f"{shap.importance[c] / total:.0%} of the edge")

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column.	required
labels	ndarray | None	Cluster id per feature, e.g. SelectionResult.labels. Defaults to one cluster per feature — fine for a handful of features, and rejected by max_clusters beyond that.	None
model_factory	Callable[[], object] | None	Zero-argument callable returning a fresh object with fit(X, y) and predict(X). Defaults to :class:RidgeRegressor.	None
scorer	Callable[[ndarray, ndarray], float] | None	(y_true, y_pred) -> float, higher is better. Defaults to the out-of-sample information coefficient.	None
n_splits	int	Cross-validation folds.	5
label_horizon	int	Bars the target looks ahead; see :func:purged_kfold.	1
embargo	int	Extra samples dropped after each test fold.	0
max_clusters	int	Refuse to run beyond this many clusters, since the cost is \(2^k\) fits per fold.	12

Returns:

Name	Type	Description
An	ImportanceResult	class:ImportanceResult whose importance values sum to
	ImportanceResult	baseline minus the empty-model score, and whose std is the spread of
	ImportanceResult	each value across folds.

Raises:

Type	Description
ValueError	If fewer than two features are given, if labels does not align with features, if too few rows survive alignment, or if the number of clusters exceeds max_clusters.

Source code in polars_ta/selection.py

def clustered_shapley(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    labels: np.ndarray | None = None,
    model_factory: Callable[[], object] | None = None,
    scorer: Callable[[np.ndarray, np.ndarray], float] | None = None,
    n_splits: int = 5,
    label_horizon: int = 1,
    embargo: int = 0,
    max_clusters: int = 12,
) -> ImportanceResult:
    r"""Exact Shapley attribution of out-of-sample score across clusters.

    MDA and SFI each answer half the question. MDA asks "what breaks if I remove
    this *last*?", SFI asks "what does this earn *first*?", and on interacting
    features those give different — sometimes wildly different — answers, with
    no principled way to reconcile them. The Shapley value (Shapley, 1953) is
    that reconciliation: it averages a cluster's marginal contribution over
    **every possible order** in which the clusters could have been added,

    $$\phi_i = \sum_{S \subseteq N \setminus \{i\}}
      \frac{|S|!\,(|N|-|S|-1)!}{|N|!}\,\bigl[v(S \cup \{i\}) - v(S)\bigr]$$

    where $v(S)$ is the out-of-sample score of a model fitted on the clusters in
    $S$. SFI is the $S = \emptyset$ term and MDA is close to the $S = N
    \setminus \{i\}$ term; Shapley is the properly weighted average of those and
    everything in between.

    What this buys you is the property neither of the others has —
    **efficiency**: the values sum exactly to $v(N) - v(\emptyset)$, the full
    model's score over an empty one. So "this cluster is 40% of the edge" is a
    statement you can defend, and the values can be compared, added and
    budgeted. It is also the only one of the four that can report a *negative*
    importance honestly: a cluster that consistently makes the model worse gets
    a negative $\phi$, not a floor of zero.

    !!! warning "Cost is exponential in the number of clusters"
        Every subset must be fitted: $2^k$ models per fold. That is 32 for
        `k=5` and 4096 for `k=12`, times `n_splits`. `max_clusters` guards
        against an accidental 200-cluster call that would never return; raise it
        only if you know what you are asking for. This is why the function takes
        **clusters** rather than features — running it per-feature on a wide
        matrix is not merely slow, it is infeasible, and cluster-level
        attribution is the more meaningful question anyway.

    ```python
    shap = selection.clustered_shapley(
        scored, sel.names, "fwd",
        labels=sel.labels, label_horizon=12, embargo=100,
    )
    total = sum(shap.importance.values())
    for c in shap.ranked:
        print(shap.clusters[c], f"{shap.importance[c] / total:.0%} of the edge")
    ```

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column.
        labels: Cluster id per feature, e.g. `SelectionResult.labels`. Defaults
            to one cluster per feature — fine for a handful of features, and
            rejected by `max_clusters` beyond that.
        model_factory: Zero-argument callable returning a fresh object with
            `fit(X, y)` and `predict(X)`. Defaults to :class:`RidgeRegressor`.
        scorer: `(y_true, y_pred) -> float`, higher is better. Defaults to the
            out-of-sample information coefficient.
        n_splits: Cross-validation folds.
        label_horizon: Bars the target looks ahead; see :func:`purged_kfold`.
        embargo: Extra samples dropped after each test fold.
        max_clusters: Refuse to run beyond this many clusters, since the cost is
            $2^k$ fits per fold.

    Returns:
        An :class:`ImportanceResult` whose `importance` values sum to
        `baseline` minus the empty-model score, and whose `std` is the spread of
        each value across folds.

    Raises:
        ValueError: If fewer than two features are given, if `labels` does not
            align with `features`, if too few rows survive alignment, or if the
            number of clusters exceeds `max_clusters`.
    """
    features, labels = _check_features_and_labels(features, labels)
    unique = np.unique(labels)
    k = unique.size
    if k > max_clusters:
        raise ValueError(
            f"{k} clusters would need {2**k} model fits per fold; raise "
            f"max_clusters (currently {max_clusters}) only if you mean it, or "
            "pass coarser `labels`"
        )

    x, y = _aligned_xy(df, features, target, n_splits)
    model_factory = RidgeRegressor if model_factory is None else model_factory
    scorer = _information_coefficient if scorer is None else scorer
    splits = purged_kfold(x.shape[0], n_splits, label_horizon, embargo)

    columns_of = {i: np.flatnonzero(labels == c) for i, c in enumerate(unique)}
    weight = {
        size: math.factorial(size) * math.factorial(k - size - 1) / math.factorial(k)
        for size in range(k)
    }

    values: dict[int, list[float]] = {int(c): [] for c in unique}
    baselines = []
    for train, test in splits:
        # Score every subset once, keyed by its bitmask, then read each marginal
        # contribution off the table — fitting per (subset, cluster) pair would
        # refit the same models k times over.
        scored_subset: dict[int, float] = {}
        for mask in range(1 << k):
            columns = np.concatenate(
                [columns_of[i] for i in range(k) if mask >> i & 1]
                or [np.empty(0, dtype=int)]
            ).astype(int)
            if columns.size == 0:
                # The empty coalition: no features, so the best any model can do
                # is predict a constant, which correlates with nothing.
                scored_subset[mask] = 0.0
                continue
            model = model_factory()
            model.fit(x[np.ix_(train, columns)], y[train])
            scored_subset[mask] = scorer(
                y[test], model.predict(x[np.ix_(test, columns)])
            )
        baselines.append(scored_subset[(1 << k) - 1])

        for i, cluster in enumerate(unique):
            phi = 0.0
            for mask in range(1 << k):
                if mask >> i & 1:
                    continue  # iterate over coalitions *without* i
                size = bin(mask).count("1")
                phi += weight[size] * (
                    scored_subset[mask | (1 << i)] - scored_subset[mask]
                )
            values[int(cluster)].append(phi)

    return ImportanceResult(
        clusters=_clusters_of(features, labels),
        importance={c: float(np.mean(v)) for c, v in values.items()},
        std={c: float(np.std(v)) for c, v in values.items()},
        baseline=float(np.mean(baselines)),
        n_splits=len(splits),
    )

null_importance

null_importance(df: DataFrame, features: Sequence[str], target: str, labels: ndarray | None = None, importance_fn: Callable[..., ImportanceResult] | None = None, n_draws: int = 50, block_size: int | None = None, seed: int | None = 0, **kwargs: object) -> NullImportanceResult

Turn any importance score into a p-value, by shuffling the label.

An importance of 0.03 is not a finding until you know what 0.03 looks like when nothing is there. It usually is not zero: with 200 candidates, a flexible model and a weak label, the best cluster will always show a positive score, and picking it is how backtests get overfitted before a single trade is placed.

The fix (Altmann et al., 2010) is a null built by breaking the one link that matters — feature to label — while leaving everything else intact. The features keep their correlation structure, the model keeps its capacity, the cross-validation keeps its geometry; only the target is scrambled. Whatever importance survives that is what the method manufactures out of nothing, and a real score has to clear it.

The label is scrambled with :func:block_permute, not a plain shuffle. A forward return is autocorrelated (volatility clusters, so nearby labels have similar magnitude), and destroying that alongside the feature-label link makes the null far too easy to beat — the same trap :func:permutation_test avoids on the redundancy side.

null = selection.null_importance(
    scored, sel.names, "fwd", labels=sel.labels,
    n_draws=100, label_horizon=12, embargo=100,
)
null.report()                    # observed vs chance, per cluster
null.significant(alpha=0.05)     # the ids that survive Bonferroni

Budget the cost before you call it

This is n_draws + 1 full runs of importance_fn. With the default :func:clustered_mda and a ridge model that is seconds; with a gradient-boosted model and 100 draws it is a coffee break, and with :func:clustered_shapley it is n_draws * 2^k * n_splits fits — check that number before starting.

Parameters:

Name	Type	Description	Default
df	DataFrame	Frame holding the features and the target.	required
features	Sequence[str]	Feature column names.	required
target	str	Target column.	required
labels	ndarray | None	Cluster id per feature, e.g. SelectionResult.labels.	None
importance_fn	Callable[..., ImportanceResult] | None	The importance routine to calibrate. Any callable with :func:clustered_mda's signature — the other three in this module all qualify. Defaults to :func:clustered_mda.	None
n_draws	int	Number of shuffled-label refits. The smallest p-value obtainable is 1/(n_draws+1), so 50 draws cannot resolve anything below 0.02; use 100+ if you intend to apply a Bonferroni correction across many clusters.	50
block_size	int | None	Block length for the label permutation. Defaults to max(kwargs["label_horizon"], n_obs // 20), which keeps each block comfortably longer than the label's own horizon.	None
seed	int | None	Seed for the shuffles.	0
**kwargs	object	Forwarded to importance_fn — model_factory, scorer, n_splits, label_horizon, embargo, and so on.	{}

Returns:

Name	Type	Description
A	NullImportanceResult	class:NullImportanceResult.

Raises:

Type	Description
ValueError	If n_draws < 1, or from the underlying importance_fn.

Source code in polars_ta/selection.py

def null_importance(
    df: pl.DataFrame,
    features: Sequence[str],
    target: str,
    labels: np.ndarray | None = None,
    importance_fn: Callable[..., ImportanceResult] | None = None,
    n_draws: int = 50,
    block_size: int | None = None,
    seed: int | None = 0,
    **kwargs: object,
) -> NullImportanceResult:
    """Turn any importance score into a p-value, by shuffling the **label**.

    An importance of 0.03 is not a finding until you know what 0.03 looks like
    when nothing is there. It usually is not zero: with 200 candidates, a
    flexible model and a weak label, the *best* cluster will always show a
    positive score, and picking it is how backtests get overfitted before a
    single trade is placed.

    The fix (Altmann et al., 2010) is a null built by breaking the one link that
    matters — feature to label — while leaving everything else intact. The
    features keep their correlation structure, the model keeps its capacity, the
    cross-validation keeps its geometry; only the target is scrambled. Whatever
    importance survives that is what the *method* manufactures out of nothing,
    and a real score has to clear it.

    The label is scrambled with :func:`block_permute`, not a plain shuffle. A
    forward return is autocorrelated (volatility clusters, so nearby labels have
    similar magnitude), and destroying that alongside the feature-label link
    makes the null far too easy to beat — the same trap
    :func:`permutation_test` avoids on the redundancy side.

    ```python
    null = selection.null_importance(
        scored, sel.names, "fwd", labels=sel.labels,
        n_draws=100, label_horizon=12, embargo=100,
    )
    null.report()                    # observed vs chance, per cluster
    null.significant(alpha=0.05)     # the ids that survive Bonferroni
    ```

    !!! note "Budget the cost before you call it"
        This is `n_draws + 1` full runs of `importance_fn`. With the default
        :func:`clustered_mda` and a ridge model that is seconds; with a
        gradient-boosted model and 100 draws it is a coffee break, and with
        :func:`clustered_shapley` it is `n_draws * 2^k * n_splits` fits — check
        that number before starting.

    Args:
        df: Frame holding the features and the target.
        features: Feature column names.
        target: Target column.
        labels: Cluster id per feature, e.g. `SelectionResult.labels`.
        importance_fn: The importance routine to calibrate. Any callable with
            :func:`clustered_mda`'s signature — the other three in this module
            all qualify. Defaults to :func:`clustered_mda`.
        n_draws: Number of shuffled-label refits. The smallest p-value
            obtainable is `1/(n_draws+1)`, so 50 draws cannot resolve anything
            below 0.02; use 100+ if you intend to apply a Bonferroni correction
            across many clusters.
        block_size: Block length for the label permutation. Defaults to
            `max(kwargs["label_horizon"], n_obs // 20)`, which keeps each block
            comfortably longer than the label's own horizon.
        seed: Seed for the shuffles.
        **kwargs: Forwarded to `importance_fn` — `model_factory`, `scorer`,
            `n_splits`, `label_horizon`, `embargo`, and so on.

    Returns:
        A :class:`NullImportanceResult`.

    Raises:
        ValueError: If `n_draws < 1`, or from the underlying `importance_fn`.
    """
    if n_draws < 1:
        raise ValueError(f"n_draws must be at least 1, got {n_draws}")

    features, labels = _check_features_and_labels(features, labels)
    importance_fn = clustered_mda if importance_fn is None else importance_fn

    real = importance_fn(df, features, target, labels=labels, **kwargs)

    column = df[target].cast(pl.Float64).to_numpy()
    n_obs = column.size
    if block_size is None:
        block_size = max(int(kwargs.get("label_horizon", 1)), max(n_obs // 20, 1))
    block_size = min(block_size, n_obs)

    rng = np.random.default_rng(seed)
    draws: dict[int, list[float]] = {c: [] for c in real.importance}
    for _ in range(n_draws):
        shuffled = block_permute(column.reshape(-1, 1), block_size, rng).ravel()
        permuted = df.with_columns(pl.Series(target, shuffled))
        result = importance_fn(permuted, features, target, labels=labels, **kwargs)
        for cluster, value in result.importance.items():
            draws[cluster].append(value)

    null = {c: np.array(v) for c, v in draws.items()}
    return NullImportanceResult(
        clusters=real.clusters,
        observed=dict(real.importance),
        null=null,
        p_value={
            c: (1 + int((null[c] >= real.importance[c]).sum())) / (1 + n_draws)
            for c in null
        },
    )