--- title: "2.1.8 Multiple testing: Bonferroni, FDR and the “factor zoo”" description: "A thousand synthetic strategies, five with a real edge: how many pass p < 0.05, what Bonferroni, Holm and Benjamini-Hochberg leave, and the “factor zoo”." url: "https://pythonparatrading.com/en/python-trading-course/level-2/multiple-testing-bonferroni-fdr-factor-zoo/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/multiple-testing-bonferroni-fdr-factor-zoo.png" --- # 2.1.8 · The multiple testing problem: Bonferroni, FDR and the “factor zoo” > Whoever tests a thousand strategies finds dozens of “significant” ones even if none of them works. This lesson measures it with a thousand synthetic strategies for which it is known which have an edge, applies the classic corrections and shows their price: with four years of data, not even the good ones survive. Then it does the same with the characteristics of a stock universe. - Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/) - Module: 2.1 · Statistics and econometrics for trading - Duration: 70 min - Prerequisites: [0.7.2 · Probability and distributions, correlation versus causation and simple regression](https://pythonparatrading.com/en/python-trading-course/level-0/probability-distributions-and-regression/), [1.4.4 · Universes without survivorship bias: delistings, mergers and symbol changes](https://pythonparatrading.com/en/python-trading-course/level-1/survivorship-bias-and-point-in-time-data/) ## By the end of this lesson - Explain why a p-value loses its meaning when the best of many tests is chosen and compute how many false positives to expect. - Apply Bonferroni, Holm, Benjamini-Hochberg and Benjamini-Yekutieli with `statsmodels` and tell the family-wise error rate (FWER) from the false discovery rate (FDR). - Relate the significance threshold to the Sharpe ratio and to the years of data needed for a real edge to clear it. - Recognize the three errors in a “factor zoo”: false positives by chance, confounding between correlated characteristics and real effects that go unseen. ## 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. ```mermaid 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] ``` ```python 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()) ``` ```output 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? ```python 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) ``` ```output 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: ```python 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) ``` ```output 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 $$ k = \max \left\{ i : p_{(i)} \le \frac{i}{m}\, q \right\}. $$ 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$). ```mermaid 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`: ```python 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}) ``` ```output 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: ```python 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 ``` ```output 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 ``` *Fig. — The 40 smallest p-values of the 1,000 strategies against the Benjamini-Hochberg step-up threshold with q = 0.05.* Log scale. Forty rising points, from 0.0017 to 0.0436, all below 0.05. The highlighted Benjamini-Hochberg threshold line rises from 0.00005 to 0.002 and stays below every point, by a factor of about 17 to 34 in the first ranks: there are no discoveries. 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$: ```python 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) ``` ```output 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: ```python 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}") ``` ```output 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 ``` *Fig. — Histogram of the p-values of the 1,000 strategies in 0.05 bins: flat, as if none had an edge.* Twenty bars between 39 and 58 strategies around a dashed line at 50, what pure chance would give. The first one, highlighted, p < 0.05, has 46: it does not stand out from the rest. 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. ```python 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)) ``` ```output 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 3[^hlz]. 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: ```python 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) ``` ```output 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 ``` *Fig. — A factor zoo with a known answer: t of the monthly rank IC for nine characteristics of the synthetic universe, 2006–2025.* Nine horizontal bars with lines at ±1.96 and ±3. 12-1 momentum, 2.96; value, 3.07; low volatility, 1.06 (the three with a real premium); quality, highlighted, −4.59 with no premium of its own, through confounding with value; size, −0.82; 1-month reversal, 2.88 with no premium of its own; ESG, 1.54, and its 1- and 12-month changes, 0.97 and 0.99, placebos. 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: ```python 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) ``` ```output 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? **Solution.** ```python 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") ``` ```output 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 ``` Not even with q = 0.30 are there discoveries: the first Benjamini-Hochberg threshold would be 0.30/1,000 = 0.0003, and the smallest p-value is 0.0017; the threshold grows with the rank, but the following p-values grow faster. Raising q is accepting a larger proportion of false ones among the chosen in exchange for power, and even so it is not enough. The problem with this data is not the method: four years do not contain enough information to separate Sharpes of 0.5 to 1.5 from noise among a thousand candidates. 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? **Solution.** ```python 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()) ``` ```output 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} ``` In sample the 20 chosen have a mean Sharpe of 2.01, because they were chosen for that; out of sample, 0.10, and 12 of the 20 stay above zero, little more than chance would give. Only two of the five with an edge make the selection (`e0133` and `e0518`), and out of sample they give −0.04 and 0.85: with two years, even a real edge may not show. The other three, which were not chosen, give between 0.33 and 1.33 in the second half. It is the simplest form of validation (lesson 2.5.2) and the reason why the in-sample Sharpe of a chosen strategy is not an estimate of its future Sharpe. 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? **Solution.** That with 40 tests the p-value of the best is not 0.009. With Bonferroni, the 5% threshold would be 0.05/40 = 0.00125, a t of about 3.2; with Holm, the same for the first one. By chance, with 40 variants with no edge about 2 are expected with p < 0.05, and the best of 40 independent nulls passes p = 0.009 with a probability of $1 - 0.991^{40} \approx 30\%$. If the variants are correlated, the effective number of tests is smaller and the problem somewhat milder, but it is still there. I would ask for the log of the 40 tests, a hypothesis for why it should work, the corrected t or a Deflated Sharpe (lesson 2.5.4) and a validation on data not used to choose. [^hlz]: 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”. ## Self-assessment 1. You test 1,000 strategies with no edge at all with a 5% test. How many do you expect to come out “significant”? - A. None - B. About 50: 5% of the nulls are rejected by chance - C. About 5 - D. It depends on the Sharpe of each one *Answer B: About 50: 5% of the nulls are rejected by chance.* That is the definition of the significance level. In the course data 46 come out with p < 0.05, of which only two have a real edge; the other 44 are chance. 2. What does the Bonferroni correction control, and what does the Benjamini-Hochberg one control? - A. Both control the same thing with different formulas - B. Bonferroni, power; Benjamini-Hochberg, size - C. Bonferroni only works with independent tests and Benjamini-Hochberg with any - D. Bonferroni, the probability of making at least one false positive (FWER); Benjamini-Hochberg, the expected proportion of false positives among the discoveries (FDR) *Answer D: Bonferroni, the probability of making at least one false positive (FWER); Benjamini-Hochberg, the expected proportion of false positives among the discoveries (FDR).* Bonferroni requires p < α/m of each test. Benjamini-Hochberg sorts the p-values and accepts the largest k with p_(k) ≤ k·q/m: it allows some false positives in exchange for more power. Under arbitrary dependence, the safe version is Benjamini-Yekutieli. 3. The best of the 1,000 strategies has a measured Sharpe of 1.58 over four years and p = 0.0017. After Bonferroni it is not significant. Does it mean it has no edge? - A. No: it means that with 1,000 tests and four years there is not enough evidence. In fact it has an edge (true Sharpe 1.5), but it would need t > 4.07, about seven years of data with that Sharpe - B. Yes: Bonferroni proves it is noise - C. Yes, because the Sharpe is below 2 - D. No, because the p-value is below 0.05 *Answer A: No: it means that with 1,000 tests and four years there is not enough evidence. In fact it has an edge (true Sharpe 1.5), but it would need t > 4.07, about seven years of data with that Sharpe.* Failing to reject is not proving the null. With t ≈ Sharpe·√years, clearing 4.07 with a Sharpe of 1.5 requires (4.07/1.5)² ≈ 7.4 years; with a Sharpe of 0.5, 66 years. The corrections protect against false positives at the price of not seeing small edges. 4. In the course's stock universe, quality (ROE) has the highest t of all the characteristics (−4.59) and its true premium is zero. What explains it? - A. Chance, which a Bonferroni correction would remove - B. A data error - C. Confounding: ROE has a correlation of −0.79 with value, which does have a premium, and inherits its effect with the opposite sign. No multiple testing correction fixes it - D. That quality has a negative premium *Answer C: Confounding: ROE has a correlation of −0.79 with value, which does have a premium, and inherits its effect with the opposite sign. No multiple testing correction fixes it.* Its t of −4.59 clears any multiple testing threshold, including the t > 3 of Harvey, Liu and Zhu. Corrections deal with chance, not with confounding. Controlling only for value is not enough (quality stays at −1.80); with all the planted signals in the regression it falls to 0.17. And with real data what is needed, above all, is an economic hypothesis. 5. What do Harvey, Liu and Zhu (2016) propose for a new factor in the cross-section of returns? - A. Use t > 2, as always - B. Do not publish new factors - C. Use only data from outside the US - D. Require a higher threshold, t > 3, because hundreds of factors have been tested on the same data *Answer D: Require a higher threshold, t > 3, because hundreds of factors have been tested on the same data.* Their abstract (NBER w20592; Review of Financial Studies, 2016) argues that with so much data mining the usual criterion of t > 2 makes no sense and estimates that a new factor should clear t > 3. ## Key points - With 1,000 strategies with no edge, about 50 pass p < 0.05. In the course data 46 pass, and only two of them have a real edge. - The best of the noise looks like skill: the second-best of the 1,000 strategies (Sharpe 1.50) has no edge at all. - Bonferroni and Holm control the probability of any false positive; Benjamini-Hochberg, the proportion of false ones among the discoveries. Here none survives, not even the five good ones. - The threshold is demanding because evidence grows with the square root of time: t ≈ Sharpe·√years. A Sharpe of 0.5 needs 66 years to clear Bonferroni with 1,000 tests. - In a factor zoo chance is not the only problem: no correction solves confounding (quality versus value) or effects that cannot be seen without controls (low volatility). ## Resources ### Research - [Introduction to Quant Research](https://pythonparatrading.com/en/quant-research-2/) · Research-0023 · 8 min - [Alpha Factors and Alpha Trading in Python](https://pythonparatrading.com/en/alpha-factors-and-alpha-trading-in-python/) · Research-0007 · 5 min ### Documentation - [statsmodels: multipletests](https://www.statsmodels.org/stable/generated/statsmodels.stats.multitest.multipletests.html) · statsmodels.org ### Books - [Harvey, C. R., Liu, Y. and Zhu, H. (2016): … and the Cross-Section of Expected Returns. Review of Financial Studies 29(1), 5-68](https://doi.org/10.1093/rfs/hhv059) · doi.org - [Harvey, Liu and Zhu: … and the Cross-Section of Expected Returns (NBER working paper w20592, abstract)](https://www.nber.org/papers/w20592) · nber.org - [Benjamini, Y. and Hochberg, Y. (1995): Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. JRSS B 57(1), 289-300](https://doi.org/10.1111/j.2517-6161.1995.tb02031.x) · doi.org ### Data - [Level 2 data: 1,000 fictitious strategies, 2022-2025 (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/ruido2_estrategias.parquet) · Download · ruido2_estrategias.parquet Previous lesson: [2.1.7 · Regime detection: hidden Markov models and structural breaks](https://pythonparatrading.com/en/python-trading-course/level-2/regime-detection-hmm-structural-breaks/) Next lesson: [2.2.1 · Edge hypotheses: risk, behavior, structure and information](https://pythonparatrading.com/en/python-trading-course/level-2/trading-edge-hypothesis/)