diff --git a/ffn/core.py b/ffn/core.py index 0108922..568b0fa 100644 --- a/ffn/core.py +++ b/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) + `_. + + 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 diff --git a/tests/test_core.py b/tests/test_core.py index cd9e680..2062c82 100644 --- a/tests/test_core.py +++ b/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( {