Add calc_fdr_hurdle, the Harvey and Liu (2020) false discovery hurdle

ffn/core.py

@@ -2608,6 +2608,121 @@ def calc_deflated_sharpe_ratio(returns, trial_sharpe_ratios, rf=0.0, nperiods=No
     return scipy.stats.norm.cdf((sr - sr0) * np.sqrt(n - 1) / np.sqrt(variance_adj))
 
 
+def calc_fdr_hurdle(returns, target_fdr=0.05, n_boot=1000, seed=0, grid_step=0.05):
+    """
+    Calculates the t-statistic hurdle that a strategy search implies at a chosen
+    false discovery rate, following `Harvey and Liu (2020)
+    <https://doi.org/10.1111/jofi.12951>`_.
+
+    There is no universal cutoff. The bar a candidate must clear depends on how
+    many candidates were tried and on how correlated they were, so it has to be
+    derived from the search itself rather than read off a table. Pass the return
+    series of **every** trial that was evaluated, one per column: the count and
+    the joint behaviour of those columns set the hurdle, and a panel of only the
+    survivors will understate it.
+
+    The null is built by demeaning each column and resampling the time index,
+    drawing the **same** periods for every column so that cross-trial
+    correlation survives into the null. Every trial is treated as null when
+    counting expected false discoveries, which is conservative in the same way
+    Benjamini-Hochberg is with ``m0 = m``. The hurdle returned is the smallest
+    one for which the target also holds at every stricter hurdle, so it cannot
+    land on a dip in a noisy curve.
+
+    Because the alternative is never represented, this says nothing about
+    discoveries that were *missed*; it bounds false ones only.
+
+    Columns with no dispersion are dropped before the hurdle is computed. Such a
+    column has no t-statistic, yet floating-point residue gives it an enormous
+    finite one, and it would otherwise be counted as the strongest trial in the
+    search. Dropping changes the trial count, which is what sets the hurdle, so
+    the result is the hurdle for the trials that carry information.
+
+    Source: Harvey, C. R. and Liu, Y. (2020), "False (and Missed) Discoveries in
+    Financial Economics", Journal of Finance, 75(5), 2503-2553.
+
+    Args:
+        * returns (DataFrame): Return series of all evaluated trials, one
+            column per trial, sharing a time index.
+        * target_fdr (float): Share of the discoveries that may be false.
+            Harvey and Liu use 0.05.
+        * n_boot (int): Number of bootstrap draws used to build the null.
+        * seed (int): Seed for the bootstrap, so the hurdle is reproducible.
+        * grid_step (float): Resolution of the hurdle grid in t units.
+
+    Returns:
+        * float -- the |t| a trial must reach. Compare it against each column's
+          own t-statistic; columns at or above it are the discoveries. When no
+          trial clears the target the value sits just above the largest observed
+          |t|, marking that the panel holds no discovery rather than describing
+          the multiplicity.
+
+    """
+    if not 0.0 < target_fdr < 1.0:
+        raise ValueError("target_fdr must be between 0 and 1")
+    if n_boot < 100:
+        raise ValueError("n_boot below 100 makes the null too coarse to trust")
+    if grid_step <= 0:
+        raise ValueError("grid_step must be positive")
+
+    x = pd.DataFrame(returns).dropna(axis=0, how="any")
+    n, m = x.shape
+    if m < 2:
+        raise ValueError("need at least 2 trials; the hurdle comes from the multiplicity")
+    if n < 3:
+        return np.nan
+
+    values = x.to_numpy(dtype=float)
+
+    # A column with no dispersion has no t-statistic, but it does not arrive here as a
+    # nan: the standard deviation of a constant series is floating-point residue rather
+    # than an exact zero, so a flat column divides out to something enormous but finite
+    # and would be counted as the strongest trial in the search. Compare against the
+    # resolution of a float at the scale of that column, so a genuinely low-volatility
+    # trial still gets a number, and drop only the columns carrying no information.
+    scale = np.max(np.abs(values), axis=0)
+    keep = values.std(axis=0, ddof=1) > n * np.finfo(float).eps * np.maximum(scale, 1.0)
+    if not keep.any():
+        return np.nan
+    dropped = m - int(keep.sum())
+    values = values[:, keep]
+    m = values.shape[1]
+    if m < 2:
+        raise ValueError(f"only {m} of {m + dropped} trials carry dispersion; need 2")
+
+    root_n = np.sqrt(n)
+    obs = np.sort(np.abs(values.mean(axis=0) / (values.std(axis=0, ddof=1) / root_n)))
+
+    demeaned = values - values.mean(axis=0)
+    rng = np.random.default_rng(seed)
+    null = np.empty(n_boot * m, dtype=float)
+    for b in range(n_boot):
+        idx = rng.integers(0, n, size=n)  # one draw shared by every column
+        sample = demeaned[idx]
+        sd = sample.std(axis=0, ddof=1)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            tb = np.abs(sample.mean(axis=0) / (sd / root_n))
+        null[b * m : (b + 1) * m] = np.where(np.isfinite(tb), tb, 0.0)
+    null.sort()
+
+    grid = np.arange(0.0, max(4.0, obs[-1] + 0.5) + grid_step, grid_step)
+    # Expected false discoveries at each hurdle, and the discoveries actually seen.
+    expected_false = (len(null) - np.searchsorted(null, grid, side="left")) / n_boot
+    discoveries = m - np.searchsorted(obs, grid, side="left")
+    with np.errstate(divide="ignore", invalid="ignore"):
+        fdr = np.where(discoveries > 0, np.minimum(1.0, expected_false / discoveries), 0.0)
+
+    # Walk down from the strictest hurdle and keep the smallest one whose whole tail
+    # still meets the target: a single dip in the curve must not be enough to pass.
+    hurdle = grid[-1]
+    tail_max = 0.0
+    for h, f in zip(grid[::-1], fdr[::-1]):
+        tail_max = max(tail_max, f)
+        if tail_max <= target_fdr:
+            hurdle = h
+    return float(hurdle)
+
+
 def resample_returns(returns, func, seed=0, num_trials=100):
     """
     Resample the returns and calculate any statistic on every new sample.
@@ -2686,6 +2801,7 @@ def extend_pandas():
     PandasObject.calc_sharpe = calc_sharpe
     PandasObject.calc_sharpe_ratio = calc_sharpe
     PandasObject.calc_deflated_sharpe_ratio = calc_deflated_sharpe_ratio
+    PandasObject.calc_fdr_hurdle = calc_fdr_hurdle
     PandasObject.to_excess_returns = to_excess_returns
     PandasObject.to_ulcer_index = to_ulcer_index
     PandasObject.to_ulcer_performance_index = to_ulcer_performance_index

tests/test_core.py

@@ -1,6 +1,7 @@
 import ffn
 import pandas as pd
 import numpy as np
+import pytest
 from pytest import fixture
 from numpy.testing import assert_almost_equal as aae
 from packaging.version import Version
@@ -889,6 +890,107 @@ def test_calc_deflated_sharpe_ratio():
     aae(winner.calc_deflated_sharpe_ratio(sharpes), dsr)
 
 
+def _trial_t_stats(panel):
+    return (panel.mean() / (panel.std(ddof=1) / np.sqrt(len(panel)))).abs()
+
+
+def test_calc_fdr_hurdle():
+    np.random.seed(0)
+    n_periods = 500
+    index = pd.date_range(start="2015-01-01", periods=n_periods, freq="D")
+    # A search over trials that have no skill whatsoever
+    trials = pd.DataFrame(np.random.normal(0, 0.01, (n_periods, 60)), index=index)
+
+    hurdle = ffn.calc_fdr_hurdle(trials, n_boot=400, seed=0)
+    assert hurdle > 0
+    # Nothing in a skill-less search may clear the bar that search itself implies
+    assert (_trial_t_stats(trials) >= hurdle).sum() == 0
+    # ... even though its best trial looks significant on its own
+    assert _trial_t_stats(trials).max() > 2.0
+
+    # Seeded, and attached to pandas objects like the other metrics
+    assert ffn.calc_fdr_hurdle(trials, n_boot=400, seed=0) == hurdle
+    assert trials.calc_fdr_hurdle(n_boot=400, seed=0) == hurdle
+
+
+def test_calc_fdr_hurdle_rises_with_the_size_of_the_search():
+    # The bar has to come from the multiplicity, so it must climb as more
+    # candidates are tried against the same real signal. Shown on a panel that
+    # holds discoveries: where nothing clears the target the hurdle collapses to
+    # just above the best observed t and stops describing the search at all.
+    np.random.seed(11)
+    n_periods = 500
+    panel = np.random.normal(0, 0.01, (n_periods, 60))
+    for j in range(3):
+        panel[:, j] = panel[:, j] + 0.01 * 6.0 / np.sqrt(n_periods)
+    panel = pd.DataFrame(panel)
+
+    hurdles = [
+        ffn.calc_fdr_hurdle(panel.iloc[:, :k], n_boot=400, seed=0)
+        for k in (5, 20, 60)
+    ]
+    assert hurdles[0] < hurdles[1] < hurdles[2]
+    # The planted trials are strong enough to survive the widest search
+    survivors = _trial_t_stats(panel) >= hurdles[-1]
+    assert survivors[:3].all()
+
+
+def test_calc_fdr_hurdle_finds_planted_signal():
+    np.random.seed(1)
+    n_periods, n_trials = 500, 20
+    panel = pd.DataFrame(np.random.normal(0, 0.01, (n_periods, n_trials)))
+    # One trial with an edge far beyond what this search can produce by chance
+    panel[0] = panel[0] + 0.01 * 6.0 / np.sqrt(n_periods)
+
+    hurdle = ffn.calc_fdr_hurdle(panel, n_boot=400, seed=0)
+    discoveries = _trial_t_stats(panel) >= hurdle
+    assert discoveries[0]
+    assert discoveries.sum() == 1
+
+
+def test_calc_fdr_hurdle_zero_dispersion():
+    # A constant column has no t-statistic, but floating-point residue gives it an
+    # enormous finite one; left in, it would be the strongest trial in the search.
+    flat = pd.DataFrame({"a": [0.001] * 250, "b": [0.001] * 250})
+    assert np.isnan(ffn.calc_fdr_hurdle(flat, n_boot=200))
+
+    np.random.seed(2)
+    real = np.random.normal(0, 0.01, (250, 2))
+    mixed = pd.DataFrame({"flat": [0.001] * 250, "a": real[:, 0], "b": real[:, 1]})
+    assert np.isfinite(ffn.calc_fdr_hurdle(mixed, n_boot=200, seed=0))
+
+    # Dropping leaves too few trials to speak of multiplicity at all
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(
+            pd.DataFrame({"flat": [0.001] * 250, "a": real[:, 0]}), n_boot=200
+        )
+
+    # A genuinely low-volatility panel is not degenerate and must still get a number,
+    # across scales and lengths rather than at the single point the floor was set on.
+    for scale in (1e-12, 1e-6, 1.0, 1e3):
+        for n in (50, 250, 1000):
+            np.random.seed(3)
+            tiny = pd.DataFrame(np.random.normal(0, scale, (n, 5)))
+            assert np.isfinite(ffn.calc_fdr_hurdle(tiny, n_boot=200, seed=0))
+
+
+def test_calc_fdr_hurdle_guards():
+    np.random.seed(4)
+    panel = pd.DataFrame(np.random.normal(0, 0.01, (100, 5)))
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(panel, target_fdr=0.0)
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(panel, target_fdr=1.0)
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(panel, n_boot=50)
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(panel, grid_step=0)
+    with pytest.raises(ValueError):
+        ffn.calc_fdr_hurdle(panel.iloc[:, :1])
+    # Too few observations to form a t-statistic
+    assert np.isnan(ffn.calc_fdr_hurdle(panel.iloc[:2], n_boot=200))
+
+
 def test_calc_information_ratio_dataframe():
     returns = pd.DataFrame(
         {