The multiple testing problem: Bonferroni, FDR and the “factor zoo”.

A thousand strategies, five good ones

A p-value answers a specific question: if this strategy had no edge, what would be the probability of seeing a result as good as the one observed? With a threshold of 0.05, a strategy with no edge passes the test one time in twenty. If you test one, it is a reasonable risk. If you test a thousand and present the best, the p-value no longer answers the question that matters, which is a different one: what is the probability that the best of a thousand looks this good by chance? Almost certain.

This is the heart of honest quantitative research, and the reason for lessons 2.2.4 (logging every experiment), 2.5.4 (Deflated Sharpe and PBO) and 2.5.2 (out-of-sample validation). Here it is measured with ruido2_estrategias.parquet: the daily returns of 1,000 synthetic “strategies” from 2022 to 2025, in 50 families of 20 variants correlated with each other, like the ones that come out of sweeping the parameters of a single idea. All of them have zero edge except five, each in a different family, and ruido2_verdad.csv says which: the known answer. They are made-up data; no figure in this lesson says anything about a real strategy.

flowchart TB
  I[One idea] --> V[Twenty variants]
  V --> F[Fifty families]
  F --> T[A thousand tests]
  T --> M[The best is chosen]
  M --> P{"p < 0.05?"}
  P -- almost always --> E[Looks like a finding]

import numpy as np
import pandas as pd
from scipy import stats

ruido = pd.read_parquet("ruido2_estrategias.parquet")
verdad = pd.read_csv("ruido2_verdad.csv", index_col="estrategia")
R = ruido.pivot(index="fecha", columns="estrategia",
                values="rendimiento")
n = len(R)
sharpe = R.mean() / R.std() * np.sqrt(252)
t, p = stats.ttest_1samp(R, 0.0)
t = pd.Series(t, index=R.columns)
p = pd.Series(p, index=R.columns)
print(n, "sessions,", R.shape[1], "strategies")
print("with p < 0.05:", (p < 0.05).sum(),
      "| expected if none had an edge:", 0.05 * R.shape[1])
print("of those, with a positive mean return:",
      ((p < 0.05) & (t > 0)).sum())

1003 sessions, 1000 strategies
with p < 0.05: 46 | expected if none had an edge: 50.0
of those, with a positive mean return: 25

46 strategies clear the 5% threshold in a two-sided test, close to the 50 that pure chance would give, and almost half (21) clear it by being significantly bad. The test is the one for the mean of the daily returns (a Student's t), which is equivalent to testing the Sharpe: with $n$ sessions, $t = \overline{r}/s \cdot \sqrt{n} = \text{SR}_{\text{annual}} \cdot \sqrt{n / 252}$, that is, $t \approx$ annual Sharpe times the square root of the years. Where are the five good ones?

buenas = verdad.index[verdad["sharpe_real"] > 0]
puesto = sharpe.rank(ascending=False).astype(int)
pd.DataFrame({"sharpe_real": verdad.loc[buenas, "sharpe_real"],
              "sharpe_medido": sharpe[buenas], "t": t[buenas],
              "p": p[buenas], "puesto": puesto[buenas]}).round(4)

            sharpe_real  sharpe_medido       t       p  puesto
estrategia
e0080               0.5         0.9508  1.8969  0.0581      31
e0133               0.8         1.4586  2.9099  0.0037       3
e0258               0.7         0.6216  1.2402  0.2152     121
e0518               1.5         1.5791  3.1504  0.0017       1
e0751               0.6         0.8320  1.6599  0.0973      56

Two of the five have p < 0.05; the other three, with true Sharpes of 0.5 to 0.7, cannot be told apart from noise in four years. And at the top of the ranking the truth mixes with chance:

top = sharpe.nlargest(8)
pd.DataFrame({"sharpe": top, "p": p[top.index],
              "familia": verdad.loc[top.index, "familia"],
              "ventaja_real": verdad.loc[top.index, "sharpe_real"] > 0}
             ).round(4)

            sharpe       p familia  ventaja_real
estrategia
e0518       1.5791  0.0017     f03          True
e0598       1.5030  0.0028     f22         False
e0133       1.4586  0.0037     f37          True
e0316       1.4227  0.0046     f22         False
e0124       1.4105  0.0050     f35         False
e0462       1.4089  0.0050     f35         False
e0542       1.3857  0.0058     f22         False
e0401       1.3549  0.0070     f24         False

The figure at the top ranks the best sixty. The second of the 1,000, e0598, has a measured Sharpe of 1.50 over four years and no edge at all: it is the best of the noise. Six of the top eight are noise. An analyst who chose “the best” would have almost the same chance of ending up with an empty strategy as with the good one, and would have no way of knowing it by looking at the numbers.

Correcting: family-wise error and false discoveries

The classic corrections answer two different questions.

Family-wise error rate (FWER): the probability of making at least one false positive among the $m$ tests. Bonferroni controls it by requiring $p_i < \alpha / m$ of each one: with $m = 1{,}000$ and $\alpha = 0.05$, $p < 0.00005$. Holm does the same in a step-down fashion and with somewhat more power: it sorts the p-values from smallest to largest and compares the $k$-th with $\alpha / (m - k + 1)$, stopping at the first one that fails.

False discovery rate (FDR): the expected proportion of false positives among the discoveries. It is a milder requirement, designed for screens in which it is accepted that some of the chosen ones are false. Benjamini-Hochberg (1995) sorts the p-values, $p_{(1)} \le \dots \le p_{(m)}$, and declares the first $k$ to be discoveries, with

It guarantees FDR ≤ q with independent tests or with positive dependence; for arbitrary dependence, Benjamini-Yekutieli also divides q by $\sum_{i=1}^{m} 1/i$ (about 7.5 with $m = 1{,}000$).

flowchart TB
  Q{"Which error hurts?"} -- a single one --> W[FWER]
  Q -- a proportion --> D[FDR]
  W --> B[Bonferroni, Holm]
  D --> H{"Dependence?"}
  H -- positive --> BH[Benjamini-Hochberg]
  H -- any --> BY[Benjamini-Yekutieli]

statsmodels has all four in multipletests:

from statsmodels.stats.multitest import multipletests

es_buena = R.columns.isin(buenas)
pb = p[es_buena]
filas = {"uncorrected": {"descubrimientos": int((p < 0.05).sum()),
                         "con_ventaja": int((pb < 0.05).sum()),
                         "p_min_ajustado": p.min()}}
for metodo, nombre in [("bonferroni", "Bonferroni"), ("holm", "Holm"),
                       ("fdr_bh", "Benjamini-Hochberg"),
                       ("fdr_by", "Benjamini-Yekutieli")]:
    rechaza, p_aj, _, _ = multipletests(p, alpha=0.05, method=metodo)
    filas[nombre] = {"descubrimientos": int(rechaza.sum()),
                     "con_ventaja": int(rechaza[es_buena].sum()),
                     "p_min_ajustado": p_aj.min()}
pd.DataFrame(filas).T.astype({"descubrimientos": int,
                               "con_ventaja": int})

                     descubrimientos  con_ventaja  p_min_ajustado
uncorrected                       46            2        0.001679
Bonferroni                         0            0        1.000000
Holm                               0            0        1.000000
Benjamini-Hochberg                 0            0        0.829541
Benjamini-Yekutieli                0            0        1.000000

None survives any correction, not even the five good ones. The smallest p-value of the 1,000 (0.0017, that of e0518, which does have an edge) becomes 1 after Bonferroni and a value far above 0.05 after Benjamini-Hochberg. Looking at the Benjamini-Hochberg procedure in the top positions helps to understand why:

orden = p.sort_values()
k = np.arange(1, 9)
paso = pd.DataFrame(
    {"p": orden.values[:8],
     "umbral_bh": 0.05 * k / len(p),
     "ventaja_real": verdad.loc[orden.index[:8], "sharpe_real"]
     .to_numpy() > 0},
    index=pd.Index(orden.index[:8], name="estrategia"))
t_bonf = stats.t.isf(0.05 / len(p) / 2, n - 1)
print(f"Bonferroni requires |t| > {t_bonf:.2f}, an annual Sharpe of "
      f"{t_bonf * np.sqrt(252 / n):.2f} over {n} sessions")
print(f"uncorrected, |t| > {stats.t.isf(0.025, n - 1):.2f} is enough: "
      f"Sharpe {stats.t.isf(0.025, n - 1) * np.sqrt(252 / n):.2f}")
paso

Bonferroni requires |t| > 4.07, an annual Sharpe of 2.04 over 1003 sessions
uncorrected, |t| > 1.96 is enough: Sharpe 0.98
                   p  umbral_bh  ventaja_real
estrategia
e0518       0.001679    0.00005          True
e0598       0.002780    0.00010         False
e0133       0.003695    0.00015          True
e0316       0.004628    0.00020         False
e0124       0.004988    0.00025         False
e0462       0.005038    0.00030         False
e0542       0.005807    0.00035         False
e0401       0.006986    0.00040         False

The Benjamini-Hochberg threshold for the first p-value is $0.05/1{,}000 = 0.00005$, the same as Bonferroni's, and it grows slowly: for the eighth, 0.0004. The p-values in the top positions are between 17 and 34 times above their threshold. Translated into Sharpe: uncorrected, a Sharpe of 0.98 is enough over four years; with Bonferroni over 1,000 tests it takes 2.04, and none of the 1,000 gets there, not even the one with a true Sharpe of 1.5.

The price: power and years of data

The corrections protect against false positives and pay for it in power: the probability of detecting an edge that does exist. With $t \approx$ Sharpe · $\sqrt{\text{years}}$, the years needed for a real edge to clear a threshold are $(t_{\text{threshold}} / \text{SR})^2$:

umbrales = {"t 1.96 (one test)": 1.96, "t 3.0 (HLZ)": 3.0,
            f"t {t_bonf:.2f} (Bonferroni, 1,000)": t_bonf}
anos = pd.DataFrame({nombre: {sr: (u / sr) ** 2
                              for sr in [0.5, 0.8, 1.0, 1.5]}
                     for nombre, u in umbrales.items()})
anos.rename_axis("sharpe_verdadero").round(1)

                  t 1.96 (one test)  t 3.0 (HLZ)  t 4.07 (Bonferroni, 1,000)
sharpe_verdadero
0.5                            15.4         36.0                        66.4
0.8                             6.0         14.1                        25.9
1.0                             3.8          9.0                        16.6
1.5                             1.7          4.0                         7.4

A true Sharpe of 0.5, which would be a good strategy in practice, needs 15 years of data to pass a single test at 5% and 66 to pass Bonferroni with 1,000 tests. Even a Sharpe of 1.5 needs more than seven years. These are years of daily data assuming the edge does not change over all that time, something that almost never holds. The conclusion is not that the corrections should be given up, but that testing a thousand things should: the number of tests is a research design decision (lesson 2.2.4) and every extra test raises the bar for all of them.

What the distribution of p-values says

If no strategy had an edge, the p-values would be spread uniformly between 0 and 1: about 50 in each bin of 0.05. The strategies with an edge should pile up near 0. The histogram gives an idea of how many null hypotheses are true, without knowing which:

conteo, _ = np.histogram(p, bins=20, range=(0, 1))
print("p-values per 0.05 bin:", conteo.tolist())
pi0 = (p > 0.5).mean() / 0.5
print(f"estimated proportion of true nulls (Storey, "
      f"lambda = 0.5): {pi0:.2f}")

p-values per 0.05 bin: [46, 45, 49, 48, 54, 53, 42, 57, 58, 49, 48, 48, 56, 47, 39, 57, 47, 55, 50, 52]
estimated proportion of true nulls (Storey, lambda = 0.5): 1.00

The histogram is flat: the first bin has 46 p-values, fewer than the 50 that chance alone would give. Storey's estimate (twice the proportion of p-values above 0.5) comes out close to 1: the data are compatible with practically all the null hypotheses being true, and they are, except five out of a thousand. When there are many real edges, the first bin clearly stands out from the rest; here nothing stands out.

One last caveat: the 20 variants of each family are correlated, so 1,000 tests are not 1,000 independent tests.

familia = verdad.loc[R.columns, "familia"].to_numpy()
C = R.corr().to_numpy()
diag = np.eye(len(familia), dtype=bool)
misma = (familia[:, None] == familia[None, :]) & ~diag
print(f"mean correlation within the family: {C[misma].mean():.2f};"
      f" between families: {C[~misma & ~diag].mean():.2f}")
mejor_de_familia = sharpe.groupby(familia).max()
print("families whose best variant exceeds a Sharpe of 0.98:",
      (mejor_de_familia > 0.98).sum(), "of", len(mejor_de_familia))

mean correlation within the family: 0.50; between families: 0.00
families whose best variant exceeds a Sharpe of 0.98: 13 of 50

With a correlation of 0.5 within each family, the “effective” number of tests is between 50 and 1,000: Bonferroni is even more conservative than necessary. But the underlying problem does not change: in 13 of the 50 families the best variant clears the uncorrected threshold (a Sharpe of 0.98), and only three of those thirteen contain a strategy with an edge; in one of the three, the variant that clears it is not even the good one. Whoever only looks at the best variant of their favorite family does not know that the other nineteen existed. The techniques that take dependence into account (the Deflated Sharpe with an effective number of tests, the PBO, the bootstrap reality check) are the subject of lessons 2.5.4 and 2.5.5.

The “factor zoo”

In academic research on the cross-section of returns, the problem has a name. Harvey, Liu and Zhu (2016) compile hundreds of factors published to explain the expected returns of stocks, argue that with that much data mining the usual criterion of a t above 2 makes no sense and estimate that a new factor should clear a t of 311Note 1 · Abstract of NBER working paper w20592 (October 2014), published in the Review of Financial Studies 29(1), 2016, pp. 5-68: “a newly discovered factor needs to clear a much higher hurdle, with a t-ratio greater than 3.0”.. It is the “factor zoo”.

The level 2 synthetic stock universe makes it possible to do the exercise with the known answer (universo2_caracteristicas.parquet, one row per stock and month-end, point-in-time). The generator plants three premiums: 12-1 momentum, value (book equity over market capitalization) and low volatility. The other characteristics have no premium of their own: quality (ROE), size, one-month reversal and the ESG score, a pure placebo, with two variants (its 1- and 12-month change) like the ones anyone looking for “something” in that column would try. For each characteristic, the monthly rank IC (the Spearman correlation between the characteristic at month-end and the next month's return) and its t:

car = pd.read_parquet("universo2_caracteristicas.parquet")
car = car.sort_values(["ticker", "fecha"]).reset_index(drop=True)
g = car.groupby("ticker")
# the next month's return includes the final row of each
# delisted stock (no price, with the delisting return)
car["rent_sig"] = g["rentabilidad_mes"].shift(-1)
con_precio = car["cierre"].notna()
log_r = np.log1p(car["rentabilidad_mes"].where(con_precio))
mom = log_r.groupby(car["ticker"]).transform(
    lambda s: s.shift(1).rolling(11).sum())
car["momentum_12_1"] = np.expm1(mom)
pn = car["patrimonio_neto"].where(car["patrimonio_neto"] > 0)
car["value"] = np.log(pn / car["capitalizacion"])
car["calidad_roe"] = car["roe"]
car["tamano"] = np.log(car["capitalizacion"])
car["reversion_1m"] = car["rentabilidad_mes"]
car["esg"] = car["puntuacion_esg"]
car["esg_cambio_1m"] = g["puntuacion_esg"].diff()
car["esg_cambio_12m"] = g["puntuacion_esg"].diff(12)

# low volatility: 252-session standard deviation from daily prices
aj = pd.read_parquet("universo2_ajustados.parquet")
aj = aj.sort_values(["ticker", "fecha"])
rd = aj.groupby("ticker")["cierre_ajustado"].pct_change()
aj["vol"] = rd.groupby(aj["ticker"]).transform(
    lambda s: s.rolling(252, min_periods=126).std())
car = car.merge(aj[["fecha", "ticker", "vol"]], how="left",
                on=["fecha", "ticker"])
car["baja_vol"] = -car["vol"]
car = car[car["cierre"].notna()]    # signals only with a price

def ic_mensual(col):
    d = car.dropna(subset=[col, "rent_sig"])
    return d.groupby("fecha")[[col, "rent_sig"]].apply(
        lambda x: x[col].corr(x["rent_sig"], method="spearman"))

senales = ["momentum_12_1", "value", "baja_vol", "calidad_roe",
           "tamano", "reversion_1m", "esg", "esg_cambio_1m",
           "esg_cambio_12m"]
zoo = {}
for s in senales:
    ic = ic_mensual(s)
    t_ic = ic.mean() / ic.std() * np.sqrt(len(ic))
    zoo[s] = {"meses": len(ic), "ic_medio": ic.mean(), "t": t_ic,
              "p": 2 * stats.norm.sf(abs(t_ic))}
zoo = pd.DataFrame(zoo).T
zoo["prima_real"] = ["yes", "yes", "yes", "no", "no", "no", "placebo",
                     "placebo", "placebo"]
zoo.round(4)

                meses  ic_medio       t       p prima_real
momentum_12_1   227.0    0.0282  2.9587  0.0031        yes
value           239.0    0.0254  3.0688  0.0021        yes
baja_vol        233.0    0.0130  1.0602  0.2890        yes
calidad_roe     239.0   -0.0324 -4.5933  0.0000         no
tamano          239.0   -0.0064 -0.8239  0.4100         no
reversion_1m    238.0    0.0246  2.8755  0.0040         no
esg             239.0    0.0103  1.5370  0.1243    placebo
esg_cambio_1m   238.0    0.0058  0.9728  0.3306    placebo
esg_cambio_12m  227.0    0.0065  0.9939  0.3203    placebo

The table has the three errors of the zoo at once:

• Chance. The ESG placebo gives a t of 1.54: it does not reach 1.96, but with a few more variants (another window, another universe, another period) one of them would pass, and its p of 0.12 already “suggests” something in an undemanding report.
• Confounding. The characteristic with the highest t of the nine is quality, at −4.59, and its true premium is zero. ROE is earnings over book equity, and in this universe its rank correlation with value is −0.79: it inherits the value premium with the sign flipped. It clears the t > 3 threshold and any multiple testing correction, because it is not chance.
• The effects that go unseen. Low volatility has a planted premium and gives a t of 1.06: volatile stocks have more beta, and in raw returns the two things offset each other. It only shows up when controlling for beta (lesson 2.3.5). One-month reversal (t = 2.88) has no premium of its own: it inherits the continuation of momentum, of post-earnings drift and of the bankruptcies.

The corrections on these nine tests fix none of that:

for metodo in ["bonferroni", "fdr_bh"]:
    rechaza = multipletests(zoo["p"].astype(float), alpha=0.05,
                            method=metodo)[0]
    print(f"{metodo:10}", list(zoo.index[rechaza]))
print("pass t > 3:", list(zoo.index[zoo["t"].abs() > 3]))

def z_rangos(s):
    r = s.rank()
    return stats.norm.ppf((r - 0.5) / r.count())

def fama_macbeth(cols):
    """t of the mean of each month's cross-sectional slopes."""
    d = car.dropna(subset=cols + ["rent_sig"]).copy()
    for c in cols:
        d[c] = d.groupby("fecha")[c].transform(z_rangos)
    pend = []
    for _, x in d.groupby("fecha"):
        X = np.column_stack([np.ones(len(x)), x[cols]])
        b = np.linalg.lstsq(X, x["rent_sig"], rcond=None)[0]
        pend.append(b[1:])
    pend = pd.DataFrame(pend, columns=cols)
    return pend.mean() / pend.std() * np.sqrt(len(pend))

fm = ["momentum_12_1", "value", "baja_vol", "calidad_roe", "tamano",
      "esg"]
pd.DataFrame({
    "una_a_una": pd.concat([fama_macbeth([c]) for c in fm]),
    "value_y_calidad": fama_macbeth(["value", "calidad_roe"]),
    "todas": fama_macbeth(fm)}).reindex(fm).round(2)

bonferroni ['momentum_12_1', 'value', 'calidad_roe', 'reversion_1m']
fdr_bh     ['momentum_12_1', 'value', 'calidad_roe', 'reversion_1m']
pass t > 3: ['value', 'calidad_roe']
               una_a_una  value_y_calidad  todas
momentum_12_1       4.43              NaN   5.79
value               2.30            -0.15   2.10
baja_vol            1.02              NaN   0.95
calidad_roe        -3.26            -1.80   0.17
tamano              0.05              NaN   0.29
esg                 1.67              NaN   1.63

Bonferroni and Benjamini-Hochberg on the nine leave the same four: momentum, value, quality and one-month reversal. Two of the four have no premium of their own, and low volatility, which does have one, is missing. The t > 3 threshold leaves only value and quality: one real and one false, and momentum, which is real, is out. The Fama-MacBeth regressions (each month, a cross-sectional regression of the next month's return on the rank scores of the signals, and the t of the mean of the slopes) tell the full story. One at a time, quality gives −3.26. Controlling only for value does not fix it: with value and quality together, quality stays at −1.80 and value at −0.15. With two signals correlated at −0.79, the regression attributes the effect to the wrong one. Only with all the signals at once does quality fall to 0.17, while value stays at 2.10 and momentum rises to 5.79, because value, negatively correlated with it, no longer dampens it (part of that momentum, according to the data README, comes from the four bankruptcies, whose prices drift downward for more than a year). Low volatility is still not visible (0.95: it lacks the control for beta) and the ESG placebo gives 1.63, close to the threshold. What separates the real from the false here is the known answer; with real data, a prior economic hypothesis and other markets and periods (lesson 2.2.1).

Common mistakes

• Counting only the tests that worked. Bonferroni with $m = 3$ because “only three are presented” corrects nothing if forty were tried.
• Reading “not significant” as “does not work”. In this sample, the strategy with a true Sharpe of 0.8 comes out with p = 0.004 and the one with 0.6, with p = 0.097: similar edges on either side of the threshold. Absence of evidence is not evidence of absence.
• Choosing the correction after seeing the results. FWER or FDR, and which α, are decided beforehand.
• Correcting and forgetting confounding. Quality passes any correction.
• Trusting the best of a family. Twenty variants of the same idea are twenty tests.
• Repeating the study on the same data. Every return to the same sample is another test, even if it is not in the code. The only clean sample is the one that has not been looked at (lesson 2.5.2).

Exercises

1. Repeat Benjamini-Hochberg on the 1,000 strategies with q = 0.10, 0.20 and 0.30. How many discoveries are there in each case and how many are strategies with an edge?

for q in [0.10, 0.20, 0.30]:
    rechaza = multipletests(p, alpha=q, method="fdr_bh")[0]
    print(f"q = {q:.2f}: {rechaza.sum()} discoveries, "
          f"{rechaza[es_buena].sum()} with an edge")

q = 0.10: 0 discoveries, 0 with an edge
q = 0.20: 0 discoveries, 0 with an edge
q = 0.30: 0 discoveries, 0 with an edge

2. Pick the 20 strategies with the best Sharpe in 2022-2023 and measure their Sharpe in 2024-2025. How many stay above zero? What happens to the ones with an edge?

primera = R.loc[:"2023-12-31"]
segunda = R.loc["2024-01-01":]
sr = lambda x: x.mean() / x.std() * np.sqrt(252)
elegidas = sr(primera).nlargest(20).index
fuera = sr(segunda[elegidas])
print(f"Mean Sharpe of the 20 chosen: in sample "
      f"{sr(primera[elegidas]).mean():.2f}, out {fuera.mean():.2f}")
print("out of sample with Sharpe > 0:", (fuera > 0).sum(), "of 20")
print("with an edge among the chosen:",
      list(elegidas[elegidas.isin(buenas)]))
print("Out-of-sample Sharpe of the five good ones:",
      sr(segunda[buenas]).round(2).to_dict())

Mean Sharpe of the 20 chosen: in sample 2.01, out 0.10
out of sample with Sharpe > 0: 12 of 20
with an edge among the chosen: ['e0133', 'e0518']
Out-of-sample Sharpe of the five good ones: {'e0080': 1.33, 'e0133': -0.04, 'e0258': 0.33, 'e0518': 0.85, 'e0751': 1.01}

3. A colleague has tried 40 variants of an indicator and presents the best with t = 2.6 (p = 0.009). What would you tell them?