SPB Git
spboucher.ai

spb/wp17_uqo

Public

FIRC0002-DATA-02 AI_C-01 ECON-05 ADVE-08 DATA-04 ECON-01 (L2/L3 inputs): second price pass, SEC share re-pull, FinData verification of extreme returns, v2 panel/analysis scripts

hygiene_log.csv: 48 monthly |ret|>300% -> 37 confirmed by FinData, 11 series breaks; 104 firm-year fwd>500% -> 101 confirmed, 3 breaks.
r3k_12: 2,540 tickers with SEC share facts (2,294 dei; 149/95 fallback concepts; 26 none); sum-of-classes at latest end.
r3k_inference.py: CGM two-way, wild-cluster bootstrap, simulated fixed-b cv, Holm/BH/Romano-Wolf, DSR, permutation, sup-Wald.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Simon-Pierre Boucher committed (Nov 5, 2026) parent 6aae05b Browse files

9 changed files +1,351 −0

9 changed files
Unified Split
added results/r3k_v2/hygiene_log.csv +153 −0
@@ -0,0 +1,153 @@
1 +ticker,kind,date,fiscal_year,yahoo_ret,findata_ret,findata_first_bar,verdict,note
2 +ABAT,monthly,2020-12-31,,4.4902,,2023-09-21,series_break,FinData history starts 2023-09-21 (inside the return window)
3 +AMC,monthly,2021-01-31,,5.2547,5.2547,2013-12-19,confirmed,2020-12-31 21.19 -> 2021-01-29 132.5
4 +AR,monthly,2020-04-30,,3.1972,3.1972,2013-10-11,confirmed,2020-03-31 0.71 -> 2020-04-30 2.98
5 +ASST,monthly,2025-05-31,,10.9218,10.9157,2023-02-06,confirmed,2025-04-30 11.77 -> 2025-05-30 140.2
6 +ASTH,monthly,2016-12-31,,4.0,4.0,2015-06-04,confirmed,2016-11-30 1.5 -> 2016-12-30 7.5
7 +AVTX,monthly,2024-03-31,,3.932,3.932,2015-11-16,confirmed,2024-02-29 4.41 -> 2024-03-28 21.75
8 +BBAI,monthly,2023-01-31,,3.8516,3.8516,2021-12-09,confirmed,2022-12-30 0.674 -> 2023-01-31 3.27
9 +BETR,monthly,2023-07-31,,3.4029,,2023-08-24,series_break,FinData history starts 2023-08-24 (inside the return window)
10 +CADL,monthly,2024-04-30,,3.462,3.462,2021-07-28,confirmed,2024-03-28 1.58 -> 2024-04-30 7.05
11 +CAPR,monthly,2014-01-31,,3.1538,3.1538,2008-05-14,confirmed,2013-12-31 32.5 -> 2014-01-31 135
12 +CAPR,monthly,2020-04-30,,4.8333,4.8333,2008-05-14,confirmed,2020-03-31 1.2 -> 2020-04-30 7
13 +CAPR,monthly,2025-12-31,,4.3944,4.3944,2008-05-14,confirmed,2025-11-28 5.35 -> 2025-12-31 28.86
14 +CHRD,monthly,2020-11-30,,309.9092,,2020-11-23,series_break,FinData history starts 2020-11-23 (inside the return window)
15 +CLDX,monthly,2020-06-30,,3.8689,3.8689,2008-10-02,confirmed,2020-05-29 2.67 -> 2020-06-30 13
16 +COGT,monthly,2020-07-31,,4.8478,4.8478,2018-04-02,confirmed,2020-06-30 1.84 -> 2020-07-31 10.76
17 +DEC,monthly,2023-12-31,,16.2612,,2023-12-19,series_break,FinData history starts 2023-12-19 (inside the return window)
18 +DFTX,monthly,2019-08-31,,3.7,,2021-04-28,series_break,FinData history starts 2021-04-28 (inside the return window)
19 +DFTX,monthly,2020-03-31,,10.2,,2021-04-28,series_break,FinData history starts 2021-04-28 (inside the return window)
20 +DJT,monthly,2021-10-31,,5.809,,2021-10-01,series_break,FinData history starts 2021-10-01 (inside the return window)
21 +FCEL,monthly,2020-11-30,,4.1,4.1,2000-06-08,confirmed,2020-10-30 60 -> 2020-11-30 306
22 +FULC,monthly,2021-08-31,,3.0629,3.0629,2019-07-19,confirmed,2021-07-30 7.31 -> 2021-08-31 29.7
23 +GME,monthly,2021-01-31,,16.2505,16.2505,2002-05-13,confirmed,2020-12-31 4.71 -> 2021-01-29 81.25
24 +HYMC,monthly,2022-03-31,,6.4194,6.4194,2018-03-13,confirmed,2022-02-28 3.1 -> 2022-03-31 23
25 +IVT,monthly,2021-09-30,,14.5669,,2021-10-13,series_break,FinData history starts 2021-10-13 (inside the return window)
26 +JANX,monthly,2024-02-29,,4.6519,4.6519,2021-06-14,confirmed,2024-01-31 8.56 -> 2024-02-29 48.38
27 +KNTK,monthly,2020-11-30,,3.0539,3.0538,2020-07-01,confirmed,2020-10-30 3.094 -> 2020-11-30 12.54
28 +KODK,monthly,2020-07-31,,8.7982,8.7982,2013-11-04,confirmed,2020-06-30 2.23 -> 2020-07-31 21.85
29 +LWLG,monthly,2021-06-30,,7.4012,7.4012,2008-03-11,confirmed,2021-05-28 1.72 -> 2021-06-30 14.45
30 +MDGL,monthly,2022-12-31,,3.1399,3.1399,2007-05-17,confirmed,2022-11-30 70.11 -> 2022-12-30 290.2
31 +MXL,monthly,2026-04-30,,3.0684,3.0684,2010-03-25,confirmed,2026-03-31 17.39 -> 2026-04-30 70.75
32 +NGNE,monthly,2015-08-31,,8.2989,8.2989,2014-03-13,confirmed,2015-07-31 36.8 -> 2015-08-31 342.2
33 +OCGN,monthly,2020-12-31,,4.7367,4.7367,2014-12-04,confirmed,2020-11-30 0.319 -> 2020-12-31 1.83
34 +OCGN,monthly,2021-02-28,,5.1864,5.1864,2014-12-04,confirmed,2021-01-29 1.77 -> 2021-02-26 10.95
35 +PAYS,monthly,2010-08-31,,4.2,3.9524,2007-05-18,confirmed,2010-07-30 0.0525 -> 2010-08-31 0.26
36 +PR,monthly,2020-04-30,,3.4867,3.4873,2016-04-25,confirmed,2020-03-31 0.2286 -> 2020-04-30 1.026
37 +QUBT,monthly,2012-01-31,,3.1667,,2018-07-12,series_break,FinData history starts 2018-07-12 (inside the return window)
38 +QUBT,monthly,2013-02-28,,5.6667,,2018-07-12,series_break,FinData history starts 2018-07-12 (inside the return window)
39 +QUBT,monthly,2018-01-31,,6.5,,2018-07-12,series_break,FinData history starts 2018-07-12 (inside the return window)
40 +QUBT,monthly,2024-11-30,,5.3036,5.3036,2018-07-12,confirmed,2024-10-31 1.12 -> 2024-11-29 7.06
41 +QXO,monthly,2023-12-31,,3.6495,3.6483,2011-08-24,confirmed,2023-11-30 29.45 -> 2023-12-29 136.9
42 +RGTI,monthly,2024-12-31,,4.0033,4.0033,2022-03-03,confirmed,2024-11-29 3.05 -> 2024-12-31 15.26
43 +SBET,monthly,2025-05-31,,22.6728,22.2509,2007-05-21,confirmed,2025-04-30 3.299 -> 2025-05-30 76.7
44 +SERV,monthly,2024-07-31,,7.9846,7.9846,2024-04-18,confirmed,2024-06-28 1.95 -> 2024-07-31 17.52
45 +SMMT,monthly,2022-12-31,,4.3125,4.3125,2015-03-06,confirmed,2022-11-30 0.8 -> 2022-12-30 4.25
46 +SOUN,monthly,2024-02-29,,3.4699,3.4699,2022-04-29,confirmed,2024-01-31 1.66 -> 2024-02-29 7.42
47 +TSHA,monthly,2023-08-31,,3.7548,3.7548,2020-09-25,confirmed,2023-07-31 0.673 -> 2023-08-31 3.2
48 +UMAC,monthly,2024-11-30,,9.016,9.016,2024-02-15,confirmed,2024-10-31 1.87 -> 2024-11-29 18.73
49 +VOR,monthly,2025-06-30,,7.1,7.1041,2021-02-08,confirmed,2025-05-30 3.998 -> 2025-06-30 32.4
50 +ABAT,firm_year_fwd,2020-06-30..2021-06-30,2019.0,22.2184,,2023-09-21,series_break,FinData history starts 2023-09-21 (after formation)
51 +ACAD,firm_year_fwd,2012-06-30..2013-06-30,2011.0,9.3125,9.3125,2004-05-28,confirmed,2012-06-29 1.76 -> 2013-06-28 18.15
52 +ACRS,firm_year_fwd,2020-06-30..2021-06-30,2019.0,9.8395,9.8395,2015-10-08,confirmed,2020-06-30 1.62 -> 2021-06-30 17.56
53 +AEHR,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.4292,6.4292,2000-01-04,confirmed,2025-06-30 12.93 -> 2026-06-30 96.06
54 +AEVA,firm_year_fwd,2024-06-30..2025-06-30,2023.0,13.996,13.996,2020-02-28,confirmed,2024-06-28 2.52 -> 2025-06-30 37.79
55 +ALMS,firm_year_fwd,2025-06-30..2026-06-30,2024.0,8.38,8.38,2024-06-28,confirmed,2025-06-30 3 -> 2026-06-30 28.14
56 +AMC,firm_year_fwd,2020-06-30..2021-06-30,2019.0,12.2121,12.2121,2013-12-19,confirmed,2020-06-30 42.88 -> 2021-06-30 566.6
57 +ANRO,firm_year_fwd,2025-06-30..2026-06-30,2024.0,10.9955,10.9955,2024-02-05,confirmed,2025-06-30 2.2 -> 2026-06-30 26.39
58 +APPS,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.0485,5.0485,2007-09-20,confirmed,2020-06-30 12.57 -> 2021-06-30 76.03
59 +ARWR,firm_year_fwd,2013-06-30..2014-06-30,2012.0,6.155,6.155,2004-09-07,confirmed,2013-06-28 2 -> 2014-06-30 14.31
60 +ARWR,firm_year_fwd,2017-06-30..2018-06-30,2016.0,7.3951,7.3951,2004-09-07,confirmed,2017-06-30 1.62 -> 2018-06-29 13.6
61 +AVTX,firm_year_fwd,2017-06-30..2018-06-30,2016.0,6.5742,6.5737,2015-11-16,confirmed,2017-06-30 1650 -> 2018-06-29 1.25e+04
62 +AXSM,firm_year_fwd,2018-06-30..2019-06-30,2017.0,7.0469,7.0469,2015-11-20,confirmed,2018-06-29 3.2 -> 2019-06-28 25.75
63 +AXTI,firm_year_fwd,2025-06-30..2026-06-30,2024.0,33.488,33.488,2000-01-04,confirmed,2025-06-30 2.09 -> 2026-06-30 72.08
64 +BBW,firm_year_fwd,2020-06-30..2021-06-30,2019.0,6.9404,6.9403,2004-10-29,confirmed,2020-06-30 1.804 -> 2021-06-30 14.32
65 +BCRX,firm_year_fwd,2013-06-30..2014-06-30,2012.0,7.2258,7.2258,2000-01-04,confirmed,2013-06-28 1.55 -> 2014-06-30 12.75
66 +BE,firm_year_fwd,2025-06-30..2026-06-30,2024.0,11.6547,11.6547,2018-07-26,confirmed,2025-06-30 23.92 -> 2026-06-30 302.7
67 +BW,firm_year_fwd,2025-06-30..2026-06-30,2024.0,13.657,13.6554,2015-07-01,confirmed,2025-06-30 0.9621 -> 2026-06-30 14.1
68 +CATX,firm_year_fwd,2013-06-30..2014-06-30,2012.0,5.5,5.5,2007-05-17,confirmed,2013-06-28 4.8 -> 2014-06-30 31.2
69 +CDNA,firm_year_fwd,2017-06-30..2018-06-30,2016.0,10.027,10.027,2014-07-18,confirmed,2017-06-30 1.11 -> 2018-06-29 12.24
70 +CELC,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.8367,6.8367,2017-09-21,confirmed,2025-06-30 13.35 -> 2026-06-30 104.6
71 +CELH,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.4647,5.4648,2010-02-11,confirmed,2020-06-30 3.923 -> 2021-06-30 25.36
72 +CHRD,firm_year_fwd,2020-06-30..2021-06-30,2019.0,134.6164,,2020-11-23,series_break,FinData history starts 2020-11-23 (after formation)
73 +CIEN,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.0317,5.0317,2000-01-04,confirmed,2025-06-30 81.33 -> 2026-06-30 490.6
74 +CLSK,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.4247,5.4247,2019-12-11,confirmed,2020-06-30 2.59 -> 2021-06-30 16.64
75 +CLYM,firm_year_fwd,2025-06-30..2026-06-30,2024.0,9.5161,9.5161,2021-08-11,confirmed,2025-06-30 1.24 -> 2026-06-30 13.04
76 +CRMD,firm_year_fwd,2018-06-30..2019-06-30,2017.0,6.475,6.475,2010-05-14,confirmed,2018-06-29 1.2 -> 2019-06-28 8.97
77 +DAVE,firm_year_fwd,2024-06-30..2025-06-30,2023.0,7.8584,7.8584,2022-01-07,confirmed,2024-06-28 30.3 -> 2025-06-30 268.4
78 +DDS,firm_year_fwd,2020-06-30..2021-06-30,2019.0,6.0822,6.0822,2000-01-04,confirmed,2020-06-30 19.83 -> 2021-06-30 140.5
79 +DFTX,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.2481,6.2481,2021-04-28,confirmed,2025-06-30 6.49 -> 2026-06-30 47.04
80 +DMRA,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.6276,6.6276,2020-11-03,confirmed,2025-06-30 3.41 -> 2026-06-30 26.01
81 +ENPH,firm_year_fwd,2017-06-30..2018-06-30,2016.0,6.7356,6.7356,2012-04-02,confirmed,2017-06-30 0.87 -> 2018-06-29 6.73
82 +ERAS,firm_year_fwd,2025-06-30..2026-06-30,2024.0,13.4252,13.4252,2021-07-19,confirmed,2025-06-30 1.27 -> 2026-06-30 18.32
83 +FCEL,firm_year_fwd,2019-06-30..2020-06-30,2018.0,11.5556,11.5555,2000-06-08,confirmed,2019-06-28 5.4 -> 2020-06-30 67.8
84 +FCEL,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.4189,5.4189,2000-06-08,confirmed,2025-06-30 5.61 -> 2026-06-30 36.01
85 +GEVO,firm_year_fwd,2020-06-30..2021-06-30,2019.0,8.2025,8.2025,2011-02-10,confirmed,2020-06-30 0.79 -> 2021-06-30 7.27
86 +GME,firm_year_fwd,2020-06-30..2021-06-30,2019.0,48.341,48.341,2002-05-13,confirmed,2020-06-30 1.085 -> 2021-06-30 53.53
87 +GSAT,firm_year_fwd,2013-06-30..2014-06-30,2012.0,6.7273,6.7272,2006-11-03,confirmed,2013-06-28 8.25 -> 2014-06-30 63.75
88 +HUT,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.2067,5.2067,2021-06-16,confirmed,2025-06-30 18.6 -> 2026-06-30 115.4
89 +HYMC,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.476,6.476,2018-03-13,confirmed,2025-06-30 3.13 -> 2026-06-30 23.4
90 +IBRX,firm_year_fwd,2019-06-30..2020-06-30,2018.0,11.0392,11.0392,2019-06-03,confirmed,2019-06-28 1.02 -> 2020-06-30 12.28
91 +INBX,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.6398,5.6398,2020-08-19,confirmed,2025-06-30 14.27 -> 2026-06-30 94.75
92 +INTC,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.2335,5.2335,2000-01-04,confirmed,2025-06-30 22.4 -> 2026-06-30 139.6
93 +IONQ,firm_year_fwd,2024-06-30..2025-06-30,2023.0,5.1124,5.1124,2021-01-05,confirmed,2024-06-28 7.03 -> 2025-06-30 42.97
94 +KOD,firm_year_fwd,2025-06-30..2026-06-30,2024.0,9.445,9.445,2018-10-05,confirmed,2025-06-30 3.73 -> 2026-06-30 38.96
95 +KOPN,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.1504,5.1504,2000-01-04,confirmed,2020-06-30 1.33 -> 2021-06-30 8.18
96 +LITE,firm_year_fwd,2025-06-30..2026-06-30,2024.0,8.0265,8.0265,2015-07-27,confirmed,2025-06-30 95.06 -> 2026-06-30 858.1
97 +LQDA,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.3989,5.3989,2018-07-27,confirmed,2025-06-30 12.46 -> 2026-06-30 79.73
98 +LWLG,firm_year_fwd,2020-06-30..2021-06-30,2019.0,17.7662,17.7662,2008-03-11,confirmed,2020-06-30 0.77 -> 2021-06-30 14.45
99 +LWLG,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.629,6.629,2008-03-11,confirmed,2025-06-30 1.24 -> 2026-06-30 9.46
100 +MARA,firm_year_fwd,2020-06-30..2021-06-30,2019.0,33.4725,33.4725,2014-07-29,confirmed,2020-06-30 0.91 -> 2021-06-30 31.37
101 +MDGL,firm_year_fwd,2017-06-30..2018-06-30,2016.0,16.2011,16.2011,2007-05-17,confirmed,2017-06-30 16.26 -> 2018-06-29 279.7
102 +MNPR,firm_year_fwd,2024-06-30..2025-06-30,2023.0,8.8027,8.8027,2019-12-20,confirmed,2024-06-28 3.65 -> 2025-06-30 35.78
103 +MU,firm_year_fwd,2025-06-30..2026-06-30,2024.0,8.3879,8.3879,2000-01-04,confirmed,2025-06-30 122.9 -> 2026-06-30 1154
104 +MVST,firm_year_fwd,2024-06-30..2025-06-30,2023.0,6.9605,6.9605,2019-03-28,confirmed,2024-06-28 0.456 -> 2025-06-30 3.63
105 +MXL,firm_year_fwd,2025-06-30..2026-06-30,2024.0,8.0099,8.0099,2010-03-25,confirmed,2025-06-30 14.21 -> 2026-06-30 128
106 +NGNE,firm_year_fwd,2019-06-30..2020-06-30,2018.0,6.0339,6.0339,2014-03-13,confirmed,2019-06-28 47.2 -> 2020-06-30 332
107 +NRGV,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.569,5.5709,2021-09-09,confirmed,2025-06-30 0.7168 -> 2026-06-30 4.71
108 +NTLA,firm_year_fwd,2020-06-30..2021-06-30,2019.0,6.7027,6.7027,2016-05-09,confirmed,2020-06-30 21.02 -> 2021-06-30 161.9
109 +NUTX,firm_year_fwd,2024-06-30..2025-06-30,2023.0,19.9227,19.9262,2022-04-05,confirmed,2024-06-28 5.949 -> 2025-06-30 124.5
110 +NVAX,firm_year_fwd,2019-06-30..2020-06-30,2018.0,13.2235,13.2235,2001-07-11,confirmed,2019-06-28 5.86 -> 2020-06-30 83.35
111 +OCGN,firm_year_fwd,2020-06-30..2021-06-30,2019.0,35.5,35.5,2014-12-04,confirmed,2020-06-30 0.22 -> 2021-06-30 8.03
112 +OKLO,firm_year_fwd,2024-06-30..2025-06-30,2023.0,5.6104,5.6104,2024-05-10,confirmed,2024-06-28 8.47 -> 2025-06-30 55.99
113 +OPEN,firm_year_fwd,2025-06-30..2026-06-30,2024.0,7.6679,7.957,2020-06-23,confirmed,2025-06-30 0.5158 -> 2026-06-30 4.62
114 +ORKA,firm_year_fwd,2025-06-30..2026-06-30,2024.0,7.4897,7.4897,2009-01-29,confirmed,2025-06-30 11.21 -> 2026-06-30 95.17
115 +OVID,firm_year_fwd,2025-06-30..2026-06-30,2024.0,7.1515,7.154,2017-05-08,confirmed,2025-06-30 0.3299 -> 2026-06-30 2.69
116 +PACB,firm_year_fwd,2020-06-30..2021-06-30,2019.0,9.1362,9.1362,2010-10-28,confirmed,2020-06-30 3.45 -> 2021-06-30 34.97
117 +PATK,firm_year_fwd,2011-06-30..2012-06-30,2010.0,5.1005,5.1007,2002-02-06,confirmed,2011-06-30 0.5414 -> 2012-06-29 3.303
118 +PLUG,firm_year_fwd,2013-06-30..2014-06-30,2012.0,11.3158,11.3158,2000-01-04,confirmed,2013-06-28 0.38 -> 2014-06-30 4.68
119 +PR,firm_year_fwd,2020-06-30..2021-06-30,2019.0,6.618,6.6177,2016-04-25,confirmed,2020-06-30 0.7737 -> 2021-06-30 5.894
120 +PRAX,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.9617,6.9617,2020-10-19,confirmed,2025-06-30 42.05 -> 2026-06-30 334.8
121 +PRCH,firm_year_fwd,2024-06-30..2025-06-30,2023.0,6.8079,6.8079,2020-01-14,confirmed,2024-06-28 1.51 -> 2025-06-30 11.79
122 +PSIX,firm_year_fwd,2024-06-30..2025-06-30,2023.0,9.1538,9.1618,2013-05-28,confirmed,2024-06-28 6.365 -> 2025-06-30 64.68
123 +PVLA,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.7799,5.7799,2015-06-29,confirmed,2025-06-30 22.54 -> 2026-06-30 152.8
124 +QBTS,firm_year_fwd,2024-06-30..2025-06-30,2023.0,11.8421,11.8421,2022-08-09,confirmed,2024-06-28 1.14 -> 2025-06-30 14.64
125 +QUBT,firm_year_fwd,2024-06-30..2025-06-30,2023.0,37.1873,37.1873,2018-07-12,confirmed,2024-06-28 0.502 -> 2025-06-30 19.17
126 +QXO,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.8069,5.0547,2011-08-24,confirmed,2020-06-30 16.08 -> 2021-06-30 97.36
127 +RCAT,firm_year_fwd,2024-06-30..2025-06-30,2023.0,5.3304,5.3304,2007-05-23,confirmed,2024-06-28 1.15 -> 2025-06-30 7.28
128 +RGTI,firm_year_fwd,2024-06-30..2025-06-30,2023.0,10.0841,10.0841,2022-03-03,confirmed,2024-06-28 1.07 -> 2025-06-30 11.86
129 +RIOT,firm_year_fwd,2020-06-30..2021-06-30,2019.0,15.9685,15.9685,2016-03-31,confirmed,2020-06-30 2.22 -> 2021-06-30 37.67
130 +RKLB,firm_year_fwd,2024-06-30..2025-06-30,2023.0,6.4521,6.4521,2021-02-17,confirmed,2024-06-28 4.8 -> 2025-06-30 35.77
131 +RLMD,firm_year_fwd,2019-06-30..2020-06-30,2018.0,5.031,5.031,2015-09-09,confirmed,2019-06-28 7.42 -> 2020-06-30 44.75
132 +RLMD,firm_year_fwd,2025-06-30..2026-06-30,2024.0,10.5333,10.5276,2015-09-09,confirmed,2025-06-30 0.6003 -> 2026-06-30 6.92
133 +SEZL,firm_year_fwd,2023-06-30..2024-06-30,2022.0,5.7241,,2023-08-17,series_break,FinData history starts 2023-08-17 (after formation)
134 +SEZL,firm_year_fwd,2024-06-30..2025-06-30,2023.0,11.1911,11.1911,2023-08-17,confirmed,2024-06-28 14.7 -> 2025-06-30 179.2
135 +SIG,firm_year_fwd,2020-06-30..2021-06-30,2019.0,6.8666,6.8666,2004-11-17,confirmed,2020-06-30 9.559 -> 2021-06-30 75.2
136 +SKY,firm_year_fwd,2017-06-30..2018-06-30,2016.0,5.0668,4.9373,2000-01-04,confirmed,2017-06-30 5.902 -> 2018-06-29 35.04
137 +SLS,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.7397,5.7397,2008-03-10,confirmed,2025-06-30 2.19 -> 2026-06-30 14.76
138 +SM,firm_year_fwd,2020-06-30..2021-06-30,2019.0,5.6168,5.6167,2002-11-21,confirmed,2020-06-30 3.364 -> 2021-06-30 22.26
139 +SMCI,firm_year_fwd,2022-06-30..2023-06-30,2021.0,5.1772,5.1772,2007-03-30,confirmed,2022-06-30 4.035 -> 2023-06-30 24.93
140 +SRPT,firm_year_fwd,2012-06-30..2013-06-30,2011.0,9.0635,9.0635,2000-01-04,confirmed,2012-06-29 3.78 -> 2013-06-28 38.04
141 +STTK,firm_year_fwd,2025-06-30..2026-06-30,2024.0,7.7626,7.7648,2020-10-12,confirmed,2025-06-30 0.7918 -> 2026-06-30 6.94
142 +STX,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.7413,5.7413,2002-12-12,confirmed,2025-06-30 143 -> 2026-06-30 964.2
143 +TBCH,firm_year_fwd,2017-06-30..2018-06-30,2016.0,6.2571,6.2571,2012-03-23,confirmed,2017-06-30 2.8 -> 2018-06-29 20.32
144 +TE,firm_year_fwd,2025-06-30..2026-06-30,2024.0,6.7073,6.7073,2025-03-03,confirmed,2025-06-30 1.23 -> 2026-06-30 9.48
145 +TNGX,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.1055,5.1055,2020-12-07,confirmed,2025-06-30 5.12 -> 2026-06-30 31.26
146 +TTI,firm_year_fwd,2020-06-30..2021-06-30,2019.0,7.1887,7.1887,2000-01-04,confirmed,2020-06-30 0.53 -> 2021-06-30 4.34
147 +UAMY,firm_year_fwd,2024-06-30..2025-06-30,2023.0,5.4118,5.4118,2007-05-17,confirmed,2024-06-28 0.34 -> 2025-06-30 2.18
148 +UCTT,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.3177,5.3177,2004-03-26,confirmed,2025-06-30 22.57 -> 2026-06-30 142.6
149 +VICR,firm_year_fwd,2025-06-30..2026-06-30,2024.0,7.3726,7.3726,2000-01-04,confirmed,2025-06-30 45.36 -> 2026-06-30 379.8
150 +VISN,firm_year_fwd,2024-06-30..2025-06-30,2023.0,5.7317,5.7309,2013-10-28,confirmed,2024-06-28 0.3422 -> 2025-06-30 2.303
151 +VKTX,firm_year_fwd,2017-06-30..2018-06-30,2016.0,7.787,7.787,2015-04-30,confirmed,2017-06-30 1.08 -> 2018-06-29 9.49
152 +VSAT,firm_year_fwd,2025-06-30..2026-06-30,2024.0,5.1514,5.1514,2000-01-04,confirmed,2025-06-30 14.6 -> 2026-06-30 89.81
153 +WDC,firm_year_fwd,2025-06-30..2026-06-30,2024.0,9.0086,9.0086,2000-01-04,confirmed,2025-06-30 63.8 -> 2026-06-30 638.5
added results/r3k_v2/hygiene_summary.json +16 −0
@@ -0,0 +1,16 @@
1 +{
2 + "n_flagged_monthly": 48,
3 + "n_flagged_firm_years": 104,
4 + "verdicts_monthly": {
5 + "confirmed": 37,
6 + "series_break": 11
7 + },
8 + "verdicts_firm_year": {
9 + "confirmed": 101,
10 + "series_break": 3
11 + },
12 + "n_queries": 261,
13 + "seconds": 221,
14 + "rule": "series_break if FinData history starts inside the window or has no bar at a window end; discrepancy if |yahoo-findata| > max(0.5, 0.25|findata|); else confirmed",
15 + "second_source": "FinData OS market_bars_stock_1d_split_div (FirstRate), read-only, 2026-11-05T14:38:34.235663"
16 +}
\ No newline at end of file
added scripts/r3k_10_download_v2.py +85 −0
@@ -0,0 +1,85 @@
1 +"""r3k_10_download_v2.py — Stage D (FIRC-0002, package L3-i / L5): second price pass for the review.
2 +Writes ONLY to data/processed_r3k_v2/ (the staged author inputs in data/processed_r3k/ are read-only):
3 + prices_monthly_unadj.parquet ticker, date (month-end), close_unadj, adj_close (yfinance auto_adjust=False, 2009-06..2026-10)
4 + splits.parquet ticker, date, ratio (yf.Ticker(t).splits)
5 + prices_daily_v2.parquet ticker, date, adj_close, close_unadj (daily 2015-01-01..2024-06-30, placebo launches + 2022 event)
6 + benchmarks_daily_v2.parquet SPY, IWV, IWM daily adj_close 2015-01-01..2024-06-30
7 + benchmarks_monthly_v2.parquet SPY, IWV, IWM monthly adj_close + close_unadj
8 + ff_daily.parquet Ken French daily FF5 + MOM (public download)
9 +Price vintage = download date (recorded in download_v2.log); the author's vintage is 2026-10-28.
10 +Run: OMP_NUM_THREADS=4 nice -n 10 .venv-review/bin/python -I scripts/r3k_10_download_v2.py
11 +"""
12 +import io, time, zipfile, warnings, datetime as dt, requests, pandas as pd, numpy as np, yfinance as yf
13 +warnings.filterwarnings("ignore")
14 +OUT = "data/processed_r3k_v2"; LOG = open(f"{OUT}/download_v2.log", "a")
15 +def log(*a):
16 + s = " ".join(str(x) for x in a); print(s, flush=True); LOG.write(s + "\n"); LOG.flush()
17 +log("START", dt.datetime.now().isoformat(), "yfinance", yf.__version__)
18 +tickers = sorted(pd.read_parquet("data/raw_russell3000/annual_final").ticker.dropna().unique().tolist())
19 +log(len(tickers), "tickers")
20 +START, END = "2009-06-01", "2026-10-28"
21 +
22 +def dl(tks, start, interval):
23 + for a in range(3):
24 + try:
25 + return yf.download(tks, start=start, end=END, interval=interval, auto_adjust=False, progress=False, group_by="ticker", threads=True)
26 + except Exception as e:
27 + log(" retry", a, repr(e)[:200]); time.sleep(5)
28 + return None
29 +def collect(df, batch, cols):
30 + out = []
31 + for t in batch:
32 + try:
33 + sub = df[t] if isinstance(df.columns, pd.MultiIndex) else df
34 + s = sub[cols].dropna(how="all")
35 + if len(s):
36 + tmp = s.reset_index(); tmp.columns = ["date"] + [c.lower().replace("adj close", "adj_close").replace("close", "close_unadj") if c == "Close" else ("adj_close" if c == "Adj Close" else c.lower()) for c in cols]; tmp["ticker"] = t; out.append(tmp)
37 + except Exception as e:
38 + log(" collect fail", t, repr(e)[:100])
39 + return out
40 +
41 +BATCH = 100; monthly, daily = [], []
42 +for i in range(0, len(tickers), BATCH):
43 + b = tickers[i:i + BATCH]; log(f"batch {i // BATCH + 1}/{(len(tickers) + BATCH - 1) // BATCH} {b[0]}..{b[-1]}")
44 + m = dl(b, START, "1mo")
45 + if m is not None: monthly += collect(m, b, ["Close", "Adj Close"])
46 + d = dl(b, "2015-01-01", "1d")
47 + if d is not None: daily += collect(d, b, ["Close", "Adj Close"])
48 + time.sleep(1)
49 +m = pd.concat(monthly, ignore_index=True); m["date"] = pd.to_datetime(m["date"]).dt.tz_localize(None); m = m.sort_values(["ticker", "date"])
50 +m.to_parquet(f"{OUT}/prices_monthly_unadj.parquet", index=False); log("monthly", len(m), "rows", m.ticker.nunique(), "tickers")
51 +d = pd.concat(daily, ignore_index=True); d["date"] = pd.to_datetime(d["date"]).dt.tz_localize(None); d = d[d.date <= "2024-06-30"].sort_values(["ticker", "date"])
52 +d.to_parquet(f"{OUT}/prices_daily_v2.parquet", index=False); log("daily", len(d), "rows", d.ticker.nunique(), "tickers")
53 +
54 +bmk = ["SPY", "IWV", "IWM"]
55 +bm = dl(bmk, START, "1mo"); rows = collect(bm, bmk, ["Close", "Adj Close"]); x = pd.concat(rows); x["date"] = pd.to_datetime(x["date"]).dt.tz_localize(None); x.to_parquet(f"{OUT}/benchmarks_monthly_v2.parquet", index=False)
56 +bd = dl(bmk, "2015-01-01", "1d"); rows = collect(bd, bmk, ["Close", "Adj Close"]); x = pd.concat(rows); x["date"] = pd.to_datetime(x["date"]).dt.tz_localize(None); x = x[x.date <= "2024-06-30"]; x.to_parquet(f"{OUT}/benchmarks_daily_v2.parquet", index=False)
57 +log("benchmarks ok", len(x))
58 +
59 +# splits (one request per ticker; ~2,500 x ~0.3 s)
60 +sp = []
61 +for i, t in enumerate(tickers):
62 + try:
63 + s = yf.Ticker(t).splits
64 + if s is not None and len(s):
65 + ss = s.reset_index(); ss.columns = ["date", "ratio"]; ss["ticker"] = t; sp.append(ss)
66 + except Exception as e:
67 + log(" splits fail", t, repr(e)[:80])
68 + if (i + 1) % 250 == 0: log(f" splits {i + 1}/{len(tickers)}"); time.sleep(1)
69 +s = pd.concat(sp, ignore_index=True); s["date"] = pd.to_datetime(s["date"]).dt.tz_localize(None); s.to_parquet(f"{OUT}/splits.parquet", index=False); log("splits", len(s), "events", s.ticker.nunique(), "tickers")
70 +
71 +# Ken French daily factors (FF5 2x3 daily + momentum daily)
72 +BASE = "https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/"
73 +def fetch_daily(name):
74 + r = requests.get(BASE + name, timeout=120); r.raise_for_status(); z = zipfile.ZipFile(io.BytesIO(r.content)); txt = z.read(z.namelist()[0]).decode("latin1").splitlines()
75 + start = next(i for i, l in enumerate(txt) if len(l.split(",")[0].strip()) == 8 and l.split(",")[0].strip().isdigit()); hdr = txt[start - 1]; rows = []
76 + for l in txt[start:]:
77 + p = [x.strip() for x in l.split(",")]
78 + if len(p[0]) != 8 or not p[0].isdigit(): break
79 + rows.append(p)
80 + cols = ["date"] + [c.strip() for c in hdr.split(",")[1:]]; df = pd.DataFrame(rows, columns=cols[:len(rows[0])])
81 + for c in df.columns[1:]: df[c] = pd.to_numeric(df[c], errors="coerce") / 100.0
82 + df["date"] = pd.to_datetime(df["date"], format="%Y%m%d"); return df
83 +ff5 = fetch_daily("F-F_Research_Data_5_Factors_2x3_daily_CSV.zip").rename(columns={"Mkt-RF": "MktRF"}); mom = fetch_daily("F-F_Momentum_Factor_daily_CSV.zip"); mom.columns = ["date", "MOM"]
84 +ff = ff5.merge(mom, on="date", how="left"); ff.to_parquet(f"{OUT}/ff_daily.parquet", index=False); log("ff daily", len(ff), ff.date.min().date(), ff.date.max().date())
85 +log("DONE", dt.datetime.now().isoformat())
added scripts/r3k_12_sec_shares_v2.py +68 −0
@@ -0,0 +1,68 @@
1 +"""r3k_12_sec_shares_v2.py — Stage D (FIRC-0002, package L3-ii, cluster A: DATA-04, ECON-01, EDIT-10).
2 +Point-in-time shares outstanding from SEC XBRL company facts, re-pulled with the raw per-class rows kept
3 +(the author's r3k_01c kept only max() per filing and overwrote marketcap.parquet):
4 + 1. dei:EntityCommonStockSharesOutstanding (cover-page count, by class) for every CIK of the panel;
5 + 2. fallback for CIKs without any dei fact: us-gaap:CommonStockSharesOutstanding, then
6 + us-gaap:WeightedAverageNumberOfSharesOutstandingBasic;
7 + 3. aggregation per (cik, filed): classes = distinct values reported for the latest `end` date of that filing ->
8 + shares_sum (sum of classes, canonical), shares_max (author's rule, for comparison), n_classes; unit = 'shares' only.
9 +Writes data/processed_r3k_v2/sec_shares_raw.parquet (every row, with concept) and shares_history_v2.parquet
10 +(ticker, cik, date = filing date, shares_sum, shares_max, n_classes, shares_source). Never touches data/processed_r3k/.
11 +Rate: <= 8 requests/s (SEC fair-use limit is 10/s). Run: nice -n 10 .venv-review/bin/python -I scripts/r3k_12_sec_shares_v2.py
12 +"""
13 +import time, threading, datetime as dt, pandas as pd, numpy as np, httpx
14 +from concurrent.futures import ThreadPoolExecutor, as_completed
15 +OUT = "data/processed_r3k_v2"; R = "data/raw_russell3000"
16 +comp = pd.read_parquet(f"{R}/companies")[["cik", "ticker"]]; panel = pd.read_parquet(f"{R}/annual_final")[["cik", "ticker"]].drop_duplicates()
17 +m = pd.concat([comp, panel]).drop_duplicates(); m["cik"] = m["cik"].astype(str).str.zfill(10)
18 +ciks = sorted(m.cik.unique()); print(len(ciks), "CIKs", dt.datetime.now().isoformat(), flush=True)
19 +H = {"User-Agent": "UQO-Research simon-pierre.boucher@uqo.ca", "Accept": "application/json"}
20 +cli = httpx.Client(timeout=30, headers=H, http2=False); lock = threading.Lock(); last = [0.0]
21 +def rate():
22 + with lock:
23 + now = time.monotonic(); wait = max(0, last[0] + 0.125 - now)
24 + if wait: time.sleep(wait)
25 + last[0] = time.monotonic()
26 +CONCEPTS = [("dei", "EntityCommonStockSharesOutstanding"), ("us-gaap", "CommonStockSharesOutstanding"), ("us-gaap", "WeightedAverageNumberOfSharesOutstandingBasic")]
27 +def get(cik, tax, concept):
28 + for a in range(4):
29 + rate()
30 + try:
31 + r = cli.get(f"https://data.sec.gov/api/xbrl/companyconcept/CIK{cik}/{tax}/{concept}.json")
32 + if r.status_code == 404: return []
33 + if r.status_code in (403, 429): time.sleep(30); continue
34 + r.raise_for_status(); d = r.json(); units = d.get("units", {})
35 + rows = [(cik, tax + ":" + concept, x.get("filed"), x.get("end"), x.get("val"), x.get("form"), x.get("fy"), x.get("fp"), u) for u, lst in units.items() for x in lst]
36 + return rows
37 + except Exception:
38 + time.sleep(2)
39 + return None
40 +def one(cik):
41 + out = []; src = None
42 + for tax, concept in CONCEPTS:
43 + rows = get(cik, tax, concept)
44 + if rows is None: return cik, None, None
45 + rows = [r for r in rows if r[8] == "shares" and r[4] is not None and r[4] > 0]
46 + if rows: out = rows; src = tax + ":" + concept; break
47 + return cik, out, src
48 +rows, fail, done, srcs = [], 0, 0, {}
49 +with ThreadPoolExecutor(4) as ex:
50 + for f in as_completed([ex.submit(one, c) for c in ciks]):
51 + cik, rr, src = f.result(); done += 1
52 + if rr is None: fail += 1
53 + else: rows += rr; srcs[cik] = src
54 + if done % 250 == 0: print(f"{done}/{len(ciks)} fail={fail}", dt.datetime.now().strftime("%H:%M:%S"), flush=True)
55 +raw = pd.DataFrame(rows, columns=["cik", "concept", "filed", "end", "shares", "form", "fy", "fp", "unit"]); raw["filed"] = pd.to_datetime(raw["filed"]); raw["end"] = pd.to_datetime(raw["end"])
56 +raw.to_parquet(f"{OUT}/sec_shares_raw.parquet", index=False)
57 +# aggregation: per filing, classes at the latest `end`
58 +raw = raw.sort_values(["cik", "filed", "end"])
59 +latest_end = raw.groupby(["cik", "filed"])["end"].transform("max"); le = raw[raw["end"] == latest_end]
60 +agg = le.groupby(["cik", "filed"]).agg(shares_sum=("shares", lambda s: float(np.unique(s.values).sum())), n_classes=("shares", lambda s: int(len(np.unique(s.values)))), concept=("concept", "first")).reset_index()
61 +mx = raw.groupby(["cik", "filed"])["shares"].max().rename("shares_max").reset_index()
62 +agg = agg.merge(mx, on=["cik", "filed"]).merge(m, on="cik", how="inner").rename(columns={"filed": "date"})
63 +agg["shares_source"] = agg["concept"]; agg = agg[["ticker", "cik", "date", "shares_sum", "shares_max", "n_classes", "shares_source"]].sort_values(["ticker", "date"])
64 +agg.to_parquet(f"{OUT}/shares_history_v2.parquet", index=False)
65 +multi = agg[agg.n_classes > 1].ticker.nunique()
66 +msg = (f"DONE {dt.datetime.now().isoformat()} CIKs {len(ciks)} failures {fail} | raw rows {len(raw)} | filings {len(agg)} | tickers {agg.ticker.nunique()} | "
67 + f"source counts {pd.Series(srcs).value_counts().to_dict()} | tickers without any fact {len(set(m.ticker) - set(agg.ticker))} | tickers with >1 class on some filing {multi}")
68 +print(msg, flush=True); open(f"{OUT}/sec_shares_v2.log", "a").write(msg + "\n")
added scripts/r3k_13_hygiene_findata.py +73 −0
@@ -0,0 +1,73 @@
1 +"""r3k_13_hygiene_findata.py — Stage D (FIRC-0002, package L2, cluster B: DATA-02, AI_C-01, ECON-05, ADVE-08, EDIT-08).
2 +Second-source verification of extreme Yahoo returns against FinData OS (M3U96a, read-only SQL, symbol-scoped):
3 + (a) every monthly |ret| > 300 % in the staged monthly_returns.parquet (Yahoo adjusted close, author vintage 2026-10-28);
4 + (b) every firm-year forward 12-month return > 500 % in the staged analysis panel (the holding window July Y+1 .. June Y+2).
5 +For each flagged observation the FinData split- and dividend-adjusted daily bars (market_bars_stock_1d_split_div) give the
6 +month-end close path; the verdict rule (pre-registered here, applied mechanically):
7 + series_break FinData has no bar on or before the start of the window (history begins inside it, e.g. CHRD 2020-11-23
8 + = chapter-11 emergence) or no bars at all -> Yahoo value cannot be confirmed -> set to missing;
9 + discrepancy |yahoo - findata| > max(0.50, 0.25 x |findata|) -> set to missing;
10 + confirmed otherwise (the extreme return is real: both sources agree) -> kept.
11 +Writes results/r3k_v2/hygiene_log.csv (ticker, kind, date/window, yahoo_ret, findata_ret, findata_first_bar, verdict) and
12 +results/r3k_v2/hygiene_summary.json. Queries are logged (count, wall-clock). Run: nice -n 10 .venv-review/bin/python -I scripts/r3k_13_hygiene_findata.py
13 +"""
14 +import os, json, time, urllib.request, datetime as dt, numpy as np, pandas as pd
15 +D = "data/processed_r3k"; T = "results/r3k_v2"
16 +BASE = os.environ.get("FINDATA_URL", "https://findata.spboucher.ai"); KEY = open(os.path.expanduser("~/Desktop/Cluster/secrets/findata/api-key-claude-code.txt")).read().strip()
17 +NQ = [0]
18 +def sql(q, max_rows=5000):
19 + req = urllib.request.Request(BASE + "/api/v1/query", data=json.dumps({"sql": q, "max_rows": max_rows}).encode(), headers={"Authorization": "Bearer " + KEY, "Content-Type": "application/json"})
20 + for a in range(3):
21 + try:
22 + with urllib.request.urlopen(req, timeout=120) as r: d = json.load(r)
23 + NQ[0] += 1; cols = [c["name"] for c in d["columns"]]; return pd.DataFrame(d["rows"], columns=cols)
24 + except Exception as e:
25 + time.sleep(2 + 3 * a); err = e
26 + raise err
27 +def fd_sym(t): return t # FinData uses the dot convention (BRK.B), as does the panel
28 +def closes(sym, start, end):
29 + q = f"SELECT session_date, close FROM market_bars_stock_1d_split_div WHERE symbol='{fd_sym(sym)}' AND session_date BETWEEN '{start}' AND '{end}' ORDER BY session_date"
30 + df = sql(q);
31 + if df.empty: return df
32 + df["session_date"] = pd.to_datetime(df["session_date"]); return df
33 +def first_bar(sym):
34 + df = sql(f"SELECT min(session_date) AS d FROM market_bars_stock_1d_split_div WHERE symbol='{fd_sym(sym)}'"); v = df.iloc[0, 0] if len(df) else None
35 + return pd.Timestamp(v) if v else None
36 +def last_close_on_or_before(df, date):
37 + s = df[df.session_date <= date]; return (float(s.close.iloc[-1]), s.session_date.iloc[-1]) if len(s) else (np.nan, None)
38 +
39 +mr = pd.read_parquet(f"{D}/monthly_returns.parquet"); mr["date"] = pd.to_datetime(mr["date"])
40 +ap = pd.read_parquet(f"{D}/analysis_panel.parquet")
41 +flag_m = mr[mr.ret.abs() > 3.0].copy(); flag_y = ap[ap.fwd_ret > 5.0][["ticker", "fiscal_year", "fwd_ret", "mktcap", "size_group"]].copy()
42 +print("flagged monthly", len(flag_m), "| flagged firm-years", len(flag_y), flush=True)
43 +rows = []; t0 = time.time(); cache_first = {}
44 +def fb(sym):
45 + if sym not in cache_first: cache_first[sym] = first_bar(sym)
46 + return cache_first[sym]
47 +for r in flag_m.itertuples():
48 + end = r.date; start = (end - pd.offsets.MonthEnd(1)); fbar = fb(r.ticker); verdict = ""; fret = np.nan
49 + if fbar is None: verdict = "series_break"; note = "no FinData history"
50 + elif fbar > start: verdict = "series_break"; note = f"FinData history starts {fbar.date()} (inside the return window)"
51 + else:
52 + c = closes(r.ticker, (start - pd.Timedelta(days=10)).date(), end.date())
53 + p0, d0 = last_close_on_or_before(c, start); p1, d1 = last_close_on_or_before(c, end)
54 + if np.isnan(p0) or np.isnan(p1) or d0 == d1: verdict = "series_break"; note = "no bar at one window end"
55 + else:
56 + fret = p1 / p0 - 1; note = f"{d0.date()} {p0:.4g} -> {d1.date()} {p1:.4g}"
57 + verdict = "discrepancy" if abs(r.ret - fret) > max(0.5, 0.25 * abs(fret)) else "confirmed"
58 + rows.append({"ticker": r.ticker, "kind": "monthly", "date": end.date().isoformat(), "fiscal_year": None, "yahoo_ret": round(float(r.ret), 4), "findata_ret": None if np.isnan(fret) else round(float(fret), 4), "findata_first_bar": fbar.date().isoformat() if fbar is not None else None, "verdict": verdict, "note": note})
59 +print("monthly checks done", NQ[0], "queries", round(time.time() - t0), "s", flush=True)
60 +for r in flag_y.itertuples():
61 + Y = int(r.fiscal_year); form = pd.Timestamp(Y + 1, 6, 30); end1 = form + pd.offsets.MonthEnd(12); fbar = fb(r.ticker); fret = np.nan
62 + if fbar is None: verdict = "series_break"; note = "no FinData history"
63 + elif fbar > form: verdict = "series_break"; note = f"FinData history starts {fbar.date()} (after formation)"
64 + else:
65 + c = closes(r.ticker, (form - pd.Timedelta(days=10)).date(), end1.date()); p0, d0 = last_close_on_or_before(c, form); p1, d1 = last_close_on_or_before(c, end1)
66 + if np.isnan(p0) or np.isnan(p1) or d0 == d1: verdict = "series_break"; note = "no bar at one window end"
67 + else:
68 + fret = p1 / p0 - 1; note = f"{d0.date()} {p0:.4g} -> {d1.date()} {p1:.4g}"; verdict = "discrepancy" if abs(r.fwd_ret - fret) > max(0.5, 0.25 * abs(fret)) else "confirmed"
69 + rows.append({"ticker": r.ticker, "kind": "firm_year_fwd", "date": f"{form.date()}..{end1.date()}", "fiscal_year": Y, "yahoo_ret": round(float(r.fwd_ret), 4), "findata_ret": None if np.isnan(fret) else round(float(fret), 4), "findata_first_bar": fbar.date().isoformat() if fbar is not None else None, "verdict": verdict, "note": note})
70 +log = pd.DataFrame(rows); log.to_csv(f"{T}/hygiene_log.csv", index=False)
71 +summ = {"n_flagged_monthly": int(len(flag_m)), "n_flagged_firm_years": int(len(flag_y)), "verdicts_monthly": log[log.kind == "monthly"].verdict.value_counts().to_dict(), "verdicts_firm_year": log[log.kind == "firm_year_fwd"].verdict.value_counts().to_dict(),
72 + "n_queries": NQ[0], "seconds": round(time.time() - t0), "rule": "series_break if FinData history starts inside the window or has no bar at a window end; discrepancy if |yahoo-findata| > max(0.5, 0.25|findata|); else confirmed", "second_source": "FinData OS market_bars_stock_1d_split_div (FirstRate), read-only, " + dt.datetime.now().isoformat()}
73 +json.dump(summ, open(f"{T}/hygiene_summary.json", "w"), indent=2); print(json.dumps(summ, indent=1)); print(log.groupby(["kind", "verdict"]).size())
added scripts/r3k_14_build_panel_v2.py +133 −0
@@ -0,0 +1,133 @@
1 +"""r3k_14_build_panel_v2.py — Stage D (FIRC-0002): rebuilt Russell 3000 firm-year panel.
2 +Packages L3 (cluster A: ECON-01, EDIT-01, ECON-02, ADVE-09, THEO-08, EDIT-10, DATA-04), L2 (cluster B: DATA-02, AI_C-01,
3 +ECON-05, ADVE-08, THEO-09, EDIT-08, AI_C-11, THEO-16, DATA-12) and the size-partition part of L6 (THEO-06, AI_C-12).
4 +Changes relative to the committed r3k_02_build_panel.py (every one is tagged in the code):
5 + [A1] market cap = UNADJUSTED month-end close (yfinance auto_adjust=False, r3k_10) x point-in-time SEC shares (sum of share
6 + classes at the filing's latest `end`, r3k_12; fallback concepts for filers without dei facts); the author's
7 + split/dividend-adjusted close x max-class count is kept as mktcap_v1 for the before/after.
8 + [A2] cap screen: cap > 10x or < 0.1x the firm's trailing 12-month median cap, or > 5e12 USD before 2023 -> NaN + mktcap_flag.
9 + [A3] size groups: (i) top-1,000 by cap within cohort (author); (ii) fixed within-cohort percentile (top 42 %); N by group.
10 + [B1] return hygiene: besides the author's filter (price <= 0, non-finite, ret < -99 %), every monthly |ret| > 300 % and every
11 + firm-year forward return > 500 % whose second-source verdict (FinData OS, r3k_13 hygiene_log.csv) is series_break or
12 + discrepancy is set to missing; a firm-year whose holding or momentum window contains a removed month has that
13 + return set to missing (splice_in_window).
14 + [B2] fwd_ret_12m_strict (all 12 months valid) and n_valid_fwd beside the author's 6-of-12 zero-filled fwd_ret.
15 + [B3] panel keyed on (cik, fiscal_year): exact duplicate rows dropped; one ticker per CIK (dual-class filers: the ticker
16 + with the longest price history); dropped rows counted in panel_meta_v2.json.
17 +Reads only staged author inputs (data/raw_russell3000, data/processed_r3k) + review downloads (data/processed_r3k_v2) +
18 +results/r3k_v2/hygiene_log.csv. Writes data/processed_r3k_v2/{monthly_returns_v2,merged_panel_v2}.parquet, panel_meta_v2.json,
19 +results/r3k_v2/{mktcap_diag.json, missing_cap_table.csv, size_group_changes.csv}.
20 +Run: OMP_NUM_THREADS=4 nice -n 10 .venv-review/bin/python -I scripts/r3k_14_build_panel_v2.py
21 +"""
22 +import json, os, numpy as np, pandas as pd
23 +D = "data/processed_r3k"; V = "data/processed_r3k_v2"; R = "data/raw_russell3000"; T = "results/r3k_v2"
24 +panel = pd.read_parquet(f"{R}/annual_final")
25 +comp = pd.read_parquet(f"{R}/companies")[["cik", "ticker", "in_sp500", "in_russell3000", "sic_code", "exchange", "russell_weight_pct", "date_added"]]
26 +px = pd.read_parquet(f"{D}/prices_monthly.parquet"); bm = pd.read_parquet(f"{D}/benchmarks_monthly.parquet")
27 +unadj = pd.read_parquet(f"{V}/prices_monthly_unadj.parquet"); sh = pd.read_parquet(f"{V}/shares_history_v2.parquet")
28 +mc_v1 = pd.read_parquet(f"{D}/marketcap.parquet"); sh_v1 = pd.read_parquet(f"{D}/shares_history.parquet")
29 +hyg = pd.read_csv(f"{T}/hygiene_log.csv"); spid = pd.read_csv(f"{R}/spid_firm_year_signals.csv"); v1 = pd.read_parquet(f"{D}/merged_panel.parquet")
30 +META = {}
31 +def mend(s): return pd.to_datetime(s).dt.to_period("M").dt.to_timestamp("M")
32 +
33 +# ---------------- monthly returns (author's construction) + hygiene [B1] ----------------
34 +px = px.sort_values(["ticker", "date"]).copy(); px["date"] = mend(px["date"]); px = px.drop_duplicates(["ticker", "date"], keep="last")
35 +px["ret"] = px.groupby("ticker")["adj_close"].pct_change()
36 +bad = (px["adj_close"] <= 0) | ~np.isfinite(px["ret"]) | (px["ret"] < -0.99); px.loc[bad, "ret"] = np.nan; META["n_bad_author_filter"] = int(bad.sum())
37 +hm = hyg[(hyg.kind == "monthly") & hyg.verdict.isin(["series_break", "discrepancy"])].copy(); hm["date"] = pd.to_datetime(hm["date"])
38 +px["splice"] = False; key = set(zip(hm.ticker, hm.date)); px.loc[[(t, d) in key for t, d in zip(px.ticker, px.date)], "splice"] = True
39 +px.loc[px.splice, "ret"] = np.nan; META["n_monthly_removed_second_source"] = int(px.splice.sum()); META["n_monthly_flagged_gt300"] = int((hyg.kind == "monthly").sum()); META["n_monthly_confirmed_gt300"] = int(((hyg.kind == "monthly") & (hyg.verdict == "confirmed")).sum())
40 +# unadjusted close [A1]
41 +unadj["date"] = mend(unadj["date"]); unadj = unadj.drop_duplicates(["ticker", "date"], keep="last")[["ticker", "date", "close_unadj"]]
42 +px = px.merge(unadj, on=["ticker", "date"], how="left"); META["share_months_with_unadj_close"] = round(float(px.close_unadj.notna().mean()), 4)
43 +# shares [A1]: point-in-time, backward as-of on filing date
44 +sh = sh.sort_values(["ticker", "date"]); sh["date"] = pd.to_datetime(sh["date"]); sh = sh.drop_duplicates(["ticker", "date"], keep="last")
45 +px = pd.merge_asof(px.sort_values("date"), sh[["ticker", "date", "shares_sum", "shares_max", "shares_source"]].rename(columns={"date": "shares_date"}).sort_values("shares_date"), left_on="date", right_on="shares_date", by="ticker", direction="backward").sort_values(["ticker", "date"])
46 +px["shares_age_m"] = (px["date"] - px["shares_date"]).dt.days / 30.44
47 +px["mktcap"] = px["shares_sum"] * px["close_unadj"]
48 +# author's variable for the before/after: adjusted close x max-class SEC shares (history + constant), exactly as r3k_02
49 +cur = mc_v1.set_index("ticker")["shares"].to_dict(); sh_v1 = sh_v1.sort_values(["ticker", "date"]); sh_v1["date"] = pd.to_datetime(sh_v1["date"]); sh_v1 = sh_v1.drop_duplicates(["ticker", "date"], keep="last")
50 +tmp = pd.merge_asof(px[["ticker", "date", "adj_close"]].sort_values("date"), sh_v1.rename(columns={"shares": "sh1", "date": "d1"}).sort_values("d1"), left_on="date", right_on="d1", by="ticker", direction="backward")
51 +tmp["sh1"] = tmp["sh1"].fillna(tmp["ticker"].map(cur)); px = px.merge(tmp[["ticker", "date", "sh1"]], on=["ticker", "date"], how="left"); px["mktcap_v1"] = px["sh1"] * px["adj_close"]
52 +# cap screen [A2]
53 +px["cap_med12"] = px.groupby("ticker")["mktcap"].transform(lambda s: s.rolling(12, min_periods=3).median().shift(1))
54 +flag = ((px["mktcap"] > 10 * px["cap_med12"]) | (px["mktcap"] < 0.1 * px["cap_med12"])) & px["cap_med12"].notna()
55 +flag |= (px["mktcap"] > 5e12) & (px["date"] < "2023-01-01")
56 +px["mktcap_flag"] = flag; px.loc[flag, "mktcap"] = np.nan; META["n_months_cap_flagged"] = int(flag.sum()); META["n_tickers_cap_flagged"] = int(px[flag].ticker.nunique())
57 +px = px.sort_values(["ticker", "date"]).reset_index(drop=True); px.to_parquet(f"{V}/monthly_returns_v2.parquet", index=False)
58 +print("monthly rows", len(px), "tickers", px.ticker.nunique(), "| caps flagged", int(flag.sum()), "| splices removed", int(px.splice.sum()), flush=True)
59 +
60 +# ---------------- benchmarks ----------------
61 +bm["date"] = mend(bm["date"]); bm = bm.sort_values(["ticker", "date"]); bm["ret"] = bm.groupby("ticker")["adj_close"].pct_change(); bix = {t: g.set_index("date")["ret"] for t, g in bm.groupby("ticker")}
62 +def cum_b(t, start, end):
63 + s = bix[t]; m = s[(s.index > start) & (s.index <= end)]; return (1 + m.fillna(0)).prod() - 1 if len(m) else np.nan
64 +
65 +# ---------------- firm-year panel ----------------
66 +def cum_ret(m, minn):
67 + if m.empty or m["ret"].notna().sum() < minn: return np.nan
68 + return (1 + m["ret"].fillna(0)).prod() - 1
69 +rows = []; px_by = {t: g for t, g in px.groupby("ticker")}
70 +hy = hyg[(hyg.kind == "firm_year_fwd") & hyg.verdict.isin(["series_break", "discrepancy"])]; hy_key = set(zip(hy.ticker, hy.fiscal_year.astype(int)))
71 +for r in panel[["ticker", "fiscal_year"]].drop_duplicates().itertuples(index=False):
72 + t, Y = r.ticker, int(r.fiscal_year); d = px_by.get(t)
73 + if d is None: continue
74 + form = pd.Timestamp(Y + 1, 6, 30) + pd.offsets.MonthEnd(0); end1 = form + pd.offsets.MonthEnd(12); mom_start = form - pd.offsets.MonthEnd(12)
75 + w1 = d[(d.date > form) & (d.date <= end1)]; w6 = d[(d.date > form) & (d.date <= form + pd.offsets.MonthEnd(6))]; w24 = d[(d.date > form) & (d.date <= form + pd.offsets.MonthEnd(24))]
76 + wm = d[(d.date > mom_start) & (d.date <= form - pd.offsets.MonthEnd(1))]; win = d[(d.date > mom_start) & (d.date <= form)]
77 + fwd = cum_ret(w1, 6); fwd6 = cum_ret(w6, 3); fwd24 = cum_ret(w24, 12); mom = cum_ret(wm, 6)
78 + nvalid = int(w1["ret"].notna().sum()); strict = (1 + w1["ret"]).prod() - 1 if (nvalid == 12 and len(w1) == 12) else np.nan
79 + splice_hold = bool(w24["splice"].any()); splice_mom = bool(win["splice"].any())
80 + if splice_hold or (t, Y) in hy_key: fwd = fwd6 = fwd24 = strict = np.nan
81 + if splice_mom: mom = np.nan
82 + vol = win["ret"].std() * np.sqrt(12) if win["ret"].notna().sum() >= 6 else np.nan; fvol = w1["ret"].std() * np.sqrt(12) if w1["ret"].notna().sum() >= 6 else np.nan
83 + capd = d.loc[(d.date <= form) & d.mktcap.notna()]; cap = capd.mktcap.iloc[-1] if len(capd) else np.nan; cap_age = (form - capd.date.iloc[-1]).days / 30.44 if len(capd) else np.nan
84 + sh_age = capd.shares_age_m.iloc[-1] if len(capd) else np.nan
85 + cap1 = d.loc[(d.date <= form) & d.mktcap_v1.notna(), "mktcap_v1"]; cap1 = cap1.iloc[-1] if len(cap1) else np.nan
86 + rows.append({"ticker": t, "fiscal_year": Y, "fwd_ret": fwd, "fwd_ret_6m": fwd6, "fwd_ret_24m": fwd24, "fwd_ret_12m_strict": strict, "n_valid_fwd": nvalid,
87 + "spy_fwd": cum_b("SPY", form, end1), "iwv_fwd": cum_b("IWV", form, end1), "iwm_fwd": cum_b("IWM", form, end1), "mom_12m": mom, "vol_12m": vol, "fwd_vol": fvol,
88 + "mktcap": cap, "log_mktcap": np.log(cap) if pd.notna(cap) and cap > 0 else np.nan, "cap_age_m": cap_age, "shares_age_m": sh_age, "mktcap_v1": cap1,
89 + "splice_in_window": splice_hold or (t, Y) in hy_key, "first_price": d.date.min(), "months_listed": (form - d.date.min()).days / 30.44})
90 +rp = pd.DataFrame(rows); rp["fwd_ret_mktadj"] = rp["fwd_ret"] - rp["spy_fwd"]; rp["fwd_ret_iwvadj"] = rp["fwd_ret"] - rp["iwv_fwd"]
91 +META["n_firm_years_fwd_removed_hygiene"] = int(rp.splice_in_window.sum())
92 +
93 +# ---------------- key on (cik, fiscal_year) [B3] ----------------
94 +panel["cik"] = panel["cik"].astype(str).str.zfill(10); n0 = len(panel)
95 +panel = panel.drop_duplicates(); n1 = len(panel); META["n_exact_duplicate_rows_dropped"] = n0 - n1
96 +hist_len = px.groupby("ticker").size(); panel["_hist"] = panel.ticker.map(hist_len).fillna(0)
97 +dual = panel.groupby("cik").ticker.nunique(); dual = dual[dual > 1].index.tolist(); META["dual_class_ciks"] = {c: sorted(panel[panel.cik == c].ticker.unique().tolist()) for c in dual}
98 +keep_t = panel[panel.cik.isin(dual)].groupby("cik").apply(lambda g: g.sort_values(["_hist", "ticker"], ascending=[False, True]).ticker.iloc[0]).to_dict()
99 +panel = panel[~panel.cik.isin(dual) | (panel.ticker == panel.cik.map(keep_t))].drop(columns="_hist"); n2 = len(panel); META["n_dual_class_rows_dropped"] = n1 - n2
100 +assert not panel.duplicated(["cik", "fiscal_year"]).any(), "panel not unique on (cik, fiscal_year)"
101 +META["kept_ticker_per_dual_class_cik"] = keep_t; META["n_rows_before_dedup"] = n0; META["n_rows_after_dedup"] = n2
102 +merged = panel.merge(rp, on=["ticker", "fiscal_year"], how="left").merge(comp.drop(columns=["ticker"]).assign(cik=lambda x: x.cik.astype(str).str.zfill(10)), on="cik", how="left")
103 +spid = spid.drop(columns=["ticker"]).rename(columns={c: f"spid_{c}" for c in spid.columns if c not in ("cik", "fiscal_year")}); spid["cik"] = spid["cik"].astype(str).str.zfill(10)
104 +merged = merged.merge(spid, on=["cik", "fiscal_year"], how="left")
105 +for c in [c for c in merged.columns if c.startswith("spid_")]: merged[c] = merged[c].fillna(0)
106 +# ---------------- size groups [A3] ----------------
107 +merged["cap_rank"] = merged.groupby("fiscal_year")["mktcap"].rank(ascending=False, method="first")
108 +merged["size_group"] = np.where(merged["cap_rank"] <= 1000, "large", np.where(merged["cap_rank"].notna(), "small", None))
109 +merged["cap_pct"] = merged.groupby("fiscal_year")["mktcap"].rank(ascending=False, pct=True)
110 +merged["size_group_pct42"] = np.where(merged["cap_pct"] <= 0.42, "large", np.where(merged["cap_pct"].notna(), "small", None))
111 +merged["cap_rank_v1"] = merged.groupby("fiscal_year")["mktcap_v1"].rank(ascending=False, method="first"); merged["size_group_v1"] = np.where(merged["cap_rank_v1"] <= 1000, "large", np.where(merged["cap_rank_v1"].notna(), "small", None))
112 +merged = merged.sort_values(["cik", "fiscal_year"]).reset_index(drop=True); merged.to_parquet(f"{V}/merged_panel_v2.parquet", index=False)
113 +est = merged[merged.fwd_ret.notna() & (merged.fiscal_year <= 2024)]
114 +META.update({"n_firm_years": int(len(merged)), "n_with_fwd": int(merged.fwd_ret.notna().sum()), "n_estimation_sample": int(len(est)), "n_estimation_sample_v1": int(len(v1[v1.fwd_ret.notna() & (v1.fiscal_year <= 2024)])),
115 + "n_est_without_cap": int(est.mktcap.isna().sum()), "n_est_without_cap_v1": int(v1[v1.fwd_ret.notna() & (v1.fiscal_year <= 2024)].mktcap.isna().sum()), "n_est_strict12": int(est.fwd_ret_12m_strict.notna().sum()), "n_est_lt12_valid": int((est.n_valid_fwd < 12).sum()),
116 + "size_groups": est.size_group.value_counts(dropna=False).to_dict(), "size_groups_pct42": est.size_group_pct42.value_counts(dropna=False).to_dict(), "price_vintage_unadj": "r3k_10 download (see data/processed_r3k_v2/download_v2.log)", "shares_source": "SEC companyfacts, sum of classes at latest end (r3k_12)"})
117 +json.dump(META, open(f"{V}/panel_meta_v2.json", "w"), indent=2, default=str)
118 +# ---------------- diagnostics for the before/after [A: acceptance] ----------------
119 +diag = {}
120 +e = est.copy(); diag["log_cap_min_max_v2"] = [round(float(e.log_mktcap.min()), 2), round(float(e.log_mktcap.max()), 2)]; diag["log_cap_min_max_v1"] = [round(float(np.log(e.mktcap_v1).min()), 2), round(float(np.log(e.mktcap_v1).max()), 2)]
121 +a11 = e[(e.ticker == "AAPL") & (e.fiscal_year == 2010)]; diag["AAPL_cohort2010_cap_v2_usd"] = float(a11.mktcap.iloc[0]) if len(a11) else None; diag["AAPL_cohort2010_rank_v2"] = float(a11.cap_rank.iloc[0]) if len(a11) else None; diag["AAPL_cohort2010_cap_v1_usd"] = float(a11.mktcap_v1.iloc[0]) if len(a11) else None; diag["AAPL_cohort2010_rank_v1"] = float(a11.cap_rank_v1.iloc[0]) if len(a11) else None
122 +for t, y in [("AEP", 2010), ("STX", 2014), ("AA", 2018), ("SLS", 2014), ("ON", 2012), ("XOM", 2011)]:
123 + s = e[(e.ticker == t) & (e.fiscal_year == y)]; diag[f"{t}_{y}"] = {"cap_v1": float(s.mktcap_v1.iloc[0]) if len(s) and pd.notna(s.mktcap_v1.iloc[0]) else None, "cap_v2": float(s.mktcap.iloc[0]) if len(s) and pd.notna(s.mktcap.iloc[0]) else None}
124 +diag["n_caps_gt_3e12"] = int((e.mktcap > 3e12).sum()); diag["n_caps_gt_3e12_v1"] = int((e.mktcap_v1 > 3e12).sum())
125 +both = e.dropna(subset=["size_group", "size_group_v1"]); diag["share_firm_years_changing_size_group"] = round(float((both.size_group != both.size_group_v1).mean()), 4); diag["n_firm_years_changing_size_group"] = int((both.size_group != both.size_group_v1).sum())
126 +diag["share_est_without_cap_v2"] = round(float(e.mktcap.isna().mean()), 4); diag["share_est_without_cap_v1"] = round(float(e.mktcap_v1.isna().mean()), 4)
127 +diag["sum_vs_max_shares_share_filings_differ"] = round(float((sh.shares_sum != sh.shares_max).mean()), 4)
128 +json.dump(diag, open(f"{T}/mktcap_diag.json", "w"), indent=2); print(json.dumps(diag, indent=1))
129 +# missing-cap table by year x forward-return tercile (EDIT-10)
130 +e["fwd_terc"] = e.groupby("fiscal_year")["fwd_ret"].transform(lambda s: pd.qcut(s.rank(method="first"), 3, labels=["T1 low", "T2", "T3 high"]))
131 +mt = e.assign(miss_v2=e.mktcap.isna(), miss_v1=e.mktcap_v1.isna()).groupby(["fiscal_year", "fwd_terc"], observed=True).agg(n=("ticker", "size"), missing_cap_v1=("miss_v1", "sum"), missing_cap_v2=("miss_v2", "sum")).reset_index(); mt.to_csv(f"{T}/missing_cap_table.csv", index=False)
132 +both.assign(changed=both.size_group != both.size_group_v1).groupby("fiscal_year").agg(n=("ticker", "size"), changed=("changed", "sum"), large_v1=("size_group_v1", lambda s: (s == "large").sum()), large_v2=("size_group", lambda s: (s == "large").sum())).reset_index().to_csv(f"{T}/size_group_changes.csv", index=False)
133 +print(json.dumps(META, indent=1, default=str))
added scripts/r3k_15_analysis_v2.py +307 −0
@@ -0,0 +1,307 @@
1 +"""r3k_15_analysis_v2.py — Stage D (FIRC-0002): core asset-pricing tests on the rebuilt panel (v2 of r3k_03_analysis.py).
2 +Packages L2 (cluster B), L3 (cluster A), L4 (clusters E/M/N, inference), L7 (facets / forms, cluster H), L8 (SPID, cluster I).
3 +Reads data/processed_r3k_v2/merged_panel_v2.parquet + monthly_returns_v2.parquet, the staged factors and filing_scores.
4 +Writes results/r3k_v2/results.json, data/processed_r3k_v2/analysis_panel_v2.parquet, ls_ts_v2_*.parquet, vw_weights_diag.csv.
5 +Differences from r3k_03 (tagged): [B] winsorization WITHIN formation year (fwd_w) with the author's pooled version kept (fwd_w_pooled);
6 +every sort/spread/FM reported on raw | winsorized | median | strict-12 outcomes with plain t, NW t (n beside), 95 % CI (t and cohort
7 +bootstrap) and the top-1 % contribution; [A] z_size from within-year winsorized log cap, VW weights on the corrected cap with the
8 +single-firm weight diagnostic; [E] two-way (firm x year) clustered and wild-cluster-bootstrap inference for every pooled regression;
9 +[H] AI Risk, PC1, equal-form and 10-K composites carried as signals; [I] SPID any-tech / any-AI binary contrasts with CI and MDE.
10 +Run: OMP_NUM_THREADS=4 nice -n 10 .venv-review/bin/python -I scripts/r3k_15_analysis_v2.py
11 +"""
12 +import os, sys, json, time, warnings, numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf
13 +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))); import r3k_inference as I
14 +from linearmodels.panel import PanelOLS
15 +warnings.filterwarnings("ignore"); T0 = time.time()
16 +D = "data/processed_r3k"; V = "data/processed_r3k_v2"; R = "data/raw_russell3000"; T = "results/r3k_v2"
17 +m = pd.read_parquet(f"{V}/merged_panel_v2.parquet"); ret = pd.read_parquet(f"{V}/monthly_returns_v2.parquet")
18 +ff = pd.read_parquet(f"{D}/ff5_monthly.parquet").merge(pd.read_parquet(f"{D}/mom_monthly.parquet"), on="date", how="left")
19 +RES = {"script": "r3k_15_analysis_v2.py", "panel": "merged_panel_v2.parquet (r3k_14)"}
20 +SIGNALS = {"AI_ADOPT": "ai_adoption_wmean", "AI_INFRA": "ai_infrastructure_wmean", "AI_RISK": "ai_disruption_risk_wmean"}
21 +PRE, POST = range(2010, 2022), range(2022, 2025)
22 +df = m[m.fwd_ret.notna() & (m.fiscal_year <= 2024)].copy()
23 +def zscore(g): return (g - g.mean()) / g.std(ddof=0)
24 +def wins_within(s): return s.clip(s.quantile(.01), s.quantile(.99))
25 +def wins_pooled(s): return s.clip(s.quantile(.01), s.quantile(.99))
26 +# [B] outcomes and controls
27 +for c, w in [("fwd_ret", "fwd_w"), ("fwd_ret_mktadj", "fwd_madj_w"), ("mom_12m", "mom_w"), ("vol_12m", "vol_w"), ("log_mktcap", "logcap_w"), ("fwd_ret_12m_strict", "fwd_strict_w"), ("fwd_ret_6m", "fwd6_w"), ("fwd_ret_24m", "fwd24_w")]:
28 + df[w] = df.groupby("fiscal_year")[c].transform(wins_within)
29 +df["fwd_w_pooled"] = wins_pooled(df["fwd_ret"])
30 +for k, col in SIGNALS.items(): df[f"z_{k}"] = df.groupby("fiscal_year")[col].transform(zscore)
31 +df["z_AI_COMPOS"] = df[["z_AI_ADOPT", "z_AI_INFRA"]].mean(axis=1)
32 +# [H] PC1 of the three facets, computed before the sorts (cluster H), signed with adoption, re-standardised within year
33 +Z = df[["z_AI_ADOPT", "z_AI_INFRA", "z_AI_RISK"]].dropna(); Zc = Z - Z.mean(); U, S, Vt = np.linalg.svd(Zc.values, full_matrices=False); pc1 = Zc.values @ Vt[0]
34 +if np.corrcoef(pc1, Z["z_AI_ADOPT"])[0, 1] < 0: pc1 = -pc1
35 +df.loc[Z.index, "z_AI_PC1"] = pc1; df["z_AI_PC1"] = df.groupby("fiscal_year")["z_AI_PC1"].transform(zscore); RES["pca"] = {"evr": [round(float(x), 3) for x in (S ** 2) / np.sum(S ** 2)], "loadings": [round(float(x), 3) for x in Vt[0]]}
36 +df["z_CYBER"] = df.groupby("fiscal_year")["cybersecurity_wmean"].transform(zscore)
37 +# [H] form-specific composites (10-K / 10-Q / 8-K) from the staged filing-level scores; equal-form composite
38 +fs = pd.read_parquet(f"{R}/filing_scores"); fl = pd.read_parquet(f"{R}/filings")[["accession_number", "form_type"]]
39 +fs = fs.drop(columns=[c for c in ["form_type"] if c in fs.columns]).merge(fl, on="accession_number", how="left"); fs["fiscal_year"] = pd.to_datetime(fs["filing_date"]).dt.year
40 +by = fs.groupby(["ticker", "fiscal_year", "form_type"]).agg(ad=("ai_adoption_wmean", "mean"), inf=("ai_infrastructure_wmean", "mean"), n=("accession_number", "count")).reset_index()
41 +for form, tag in [("10-K", "10k"), ("10-Q", "10q"), ("8-K", "8k")]:
42 + b = by[by.form_type == form][["ticker", "fiscal_year", "ad", "inf", "n"]].rename(columns={"ad": f"ad_{tag}", "inf": f"inf_{tag}", "n": f"n_{tag}"})
43 + df = df.drop(columns=[c for c in df.columns if c in b.columns and c not in ("ticker", "fiscal_year")]).merge(b, on=["ticker", "fiscal_year"], how="left")
44 + df[f"z_AI_{tag}"] = (df.groupby("fiscal_year")[f"ad_{tag}"].transform(zscore) + df.groupby("fiscal_year")[f"inf_{tag}"].transform(zscore)) / 2
45 +df["z_AI_formw"] = df[["z_AI_10k", "z_AI_10q", "z_AI_8k"]].mean(axis=1, skipna=False); df["z_AI_formw"] = df.groupby("fiscal_year")["z_AI_formw"].transform(zscore)
46 +# SPID [I]
47 +df["spid_tech_log"] = np.log1p(df["spid_mentions_tech"]); df["spid_ai_any"] = (df["spid_mentions_ai"] + df["spid_mentions_ai_tech_tag"] > 0).astype(float); df["spid_tech_any"] = (df["spid_mentions_tech"] > 0).astype(float)
48 +df["z_SPID_TECH"] = df.groupby("fiscal_year")["spid_tech_log"].transform(zscore)
49 +# [A] controls: z_size from within-year winsorized log cap
50 +df["z_size"] = df.groupby("fiscal_year")["logcap_w"].transform(zscore); df["z_mom"] = df.groupby("fiscal_year")["mom_w"].transform(zscore); df["z_vol"] = df.groupby("fiscal_year")["vol_w"].transform(zscore)
51 +df["post"] = (df.fiscal_year >= 2022).astype(int); df = df.sort_values(["cik", "fiscal_year"]).reset_index(drop=True)
52 +RES["n_obs"] = int(len(df)); RES["n_firms"] = int(df.ticker.nunique()); RES["yrs"] = f"{int(df.fiscal_year.min())}-{int(df.fiscal_year.max())}"; RES["n_by_year"] = df.groupby("fiscal_year").size().to_dict()
53 +RES["n_sp500_firm_years"] = int(df.in_sp500.sum()); RES["n_large"] = int((df.size_group == "large").sum()); RES["n_without_cap"] = int(df.mktcap.isna().sum()); RES["n_strict12"] = int(df.fwd_ret_12m_strict.notna().sum())
54 +print("N", len(df), "firms", df.ticker.nunique(), RES["yrs"], "| no cap", RES["n_without_cap"], flush=True)
55 +summ_cols = list(SIGNALS.values()) + ["spid_mentions_tech", "spid_ai_any", "spid_tech_any", "fwd_ret", "fwd_w", "fwd_ret_mktadj", "mom_12m", "vol_12m", "log_mktcap", "mktcap", "total_chunks"]
56 +df[summ_cols].describe(percentiles=[.01, .05, .25, .5, .75, .95, .99]).T.to_csv(f"{T}/t1_summary_stats.csv")
57 +df[["z_AI_ADOPT", "z_AI_INFRA", "z_AI_RISK", "z_AI_PC1", "z_SPID_TECH", "z_CYBER", "z_size", "z_mom", "fwd_ret"]].corr().to_csv(f"{T}/t1b_correlations.csv")
58 +RES["corr"] = df[["z_AI_ADOPT", "z_AI_INFRA", "z_AI_RISK", "z_SPID_TECH", "z_CYBER", "z_AI_PC1"]].corr().round(3).to_dict()
59 +
60 +# ---------------- sort machinery ----------------
61 +VWD = []
62 +def spread_stats(ls, by_cohort=None):
63 + mu, t = I.nw_t(ls); out = {"Q5_minus_Q1": mu, "t_stat": t, "t_plain": I.plain_t(ls), "ci95": I.ci95(ls), "ci95_boot": I.block_boot_ci(ls), "n_years": len(ls)}
64 + if by_cohort is not None: out["by_cohort"] = {int(k): round(v, 4) for k, v in by_cohort.items()}
65 + return out
66 +def sort_table(d, z, weight=None, oc="fwd_ret", stat="mean", minn=40, yrs=None, diag_label=None):
67 + ann = {q: [] for q in range(1, 6)}; ls = []; ls_trim = []; byc = {}
68 + for yr, g in d.groupby("fiscal_year"):
69 + if yrs is not None and yr not in yrs: continue
70 + g = g.dropna(subset=[z, oc]).copy()
71 + if len(g) < minn: continue
72 + g["q"] = pd.qcut(g[z].rank(method="first"), 5, labels=False) + 1; p99 = g[oc].quantile(.99); vals = {}; vals_trim = {}
73 + for q in range(1, 6):
74 + s = g[g.q == q]
75 + if s.empty: continue
76 + if weight is None:
77 + vals[q] = s[oc].median() if stat == "median" else s[oc].mean(); st = s[s[oc] <= p99]; vals_trim[q] = st[oc].mean() if len(st) else np.nan
78 + else:
79 + w = s[weight].clip(lower=0).fillna(0); vals[q] = np.average(s[oc], weights=w) if w.sum() > 0 else s[oc].mean(); vals_trim[q] = vals[q]
80 + if diag_label and w.sum() > 0:
81 + ws = (w / w.sum()).sort_values(ascending=False); VWD.append({"sort": diag_label, "cohort": int(yr), "q": q, "n": len(s), "max_weight": round(float(ws.iloc[0]), 4), "max_ticker": s.loc[ws.index[0], "ticker"], "top10_weight": round(float(ws.iloc[:10].sum()), 4)})
82 + ann[q].append(vals[q])
83 + if 5 in vals and 1 in vals: ls.append(vals[5] - vals[1]); ls_trim.append(vals_trim[5] - vals_trim[1]); byc[int(yr)] = vals[5] - vals[1]
84 + out = {f"Q{q}": float(np.mean(ann[q])) if ann[q] else np.nan for q in range(1, 6)}; out.update(spread_stats(ls, byc)); out["spread_ex_top1pct"] = float(np.nanmean(ls_trim)) if ls_trim else np.nan; out["top1pct_contribution"] = out["Q5_minus_Q1"] - out["spread_ex_top1pct"]; return out
85 +OUTS = [("", "fwd_ret", "mean"), ("_w", "fwd_w", "mean"), ("_wp", "fwd_w_pooled", "mean"), ("_med", "fwd_ret", "median"), ("_strict", "fwd_ret_12m_strict", "mean"), ("_strict_w", "fwd_strict_w", "mean")]
86 +SIGL = ["AI_ADOPT", "AI_INFRA", "AI_RISK", "AI_COMPOS", "AI_PC1", "AI_formw", "AI_10k", "SPID_TECH", "CYBER"]
87 +sorts = {}
88 +for k in SIGL:
89 + for sfx, oc, st in OUTS:
90 + sorts[f"{k}_EW{sfx}"] = sort_table(df, f"z_{k}", None, oc, st)
91 + if sfx in ("", "_w"): sorts[f"{k}_VW{sfx}"] = sort_table(df, f"z_{k}", "mktcap", oc, st, diag_label=f"{k}_VW{sfx}"); sorts[f"{k}_VW_v1{sfx}"] = sort_table(df, f"z_{k}", "mktcap_v1", oc, st, diag_label=f"{k}_VW_v1{sfx}")
92 + sorts[f"{k}_EW_madj"] = sort_table(df, f"z_{k}", None, "fwd_ret_mktadj"); sorts[f"{k}_EW_madj_w"] = sort_table(df, f"z_{k}", None, "fwd_madj_w")
93 +RES["sorts"] = {k: {kk: (round(vv, 5) if isinstance(vv, float) else vv) for kk, vv in v.items()} for k, v in sorts.items()}
94 +pd.DataFrame(sorts).T.drop(columns=["by_cohort"]).to_csv(f"{T}/t2_portfolio_sorts.csv")
95 +vwd = pd.DataFrame(VWD); vwd.to_csv(f"{T}/vw_weights_diag.csv", index=False)
96 +RES["vw_diag"] = {"max_single_firm_weight_v2": float(vwd[vwd["sort"].str.contains("_VW$|_VW_w$", regex=True)].max_weight.max()), "max_single_firm_weight_v1": float(vwd[vwd["sort"].str.contains("_VW_v1")].max_weight.max()),
97 + "n_cells_above_25pct_v2": int((vwd[vwd["sort"].str.contains("_VW$|_VW_w$", regex=True)].max_weight > 0.25).sum()), "n_cells_above_25pct_v1": int((vwd[vwd["sort"].str.contains("_VW_v1")].max_weight > 0.25).sum()),
98 + "worst_cells_v2": vwd[vwd["sort"] == "AI_COMPOS_VW"].nlargest(5, "max_weight")[["cohort", "q", "max_weight", "max_ticker"]].to_dict("records"), "worst_cells_v1": vwd[vwd["sort"] == "AI_COMPOS_VW_v1"].nlargest(5, "max_weight")[["cohort", "q", "max_weight", "max_ticker"]].to_dict("records")}
99 +for k in ["AI_COMPOS_EW", "AI_COMPOS_EW_w", "AI_COMPOS_EW_med", "AI_COMPOS_VW", "AI_COMPOS_VW_v1", "AI_COMPOS_EW_strict"]:
100 + s = sorts[k]; print(f"{k:22s} Q1 {s['Q1']:+.4f} Q5 {s['Q5']:+.4f} spread {s['Q5_minus_Q1']:+.4f} NW t {s['t_stat']:.2f} plain t {s['t_plain']:.2f} CI {np.round(s['ci95'], 4)} n {s['n_years']}", flush=True)
101 +print("VW diag", {k: v for k, v in RES["vw_diag"].items() if not k.startswith("worst")}, flush=True)
102 +
103 +# ---------------- monthly long-short + alphas ----------------
104 +ff["date"] = pd.to_datetime(ff["date"]).dt.to_period("M").dt.to_timestamp("M"); ff_ix = ff.set_index("date")
105 +ret["date"] = pd.to_datetime(ret["date"]).dt.to_period("M").dt.to_timestamp("M"); ret_by = {t: g.set_index("date")["ret"] for t, g in ret.groupby("ticker")}
106 +def ls_monthly(d, z, n=5, minn=40):
107 + monthly = {}
108 + for yr, g in d.groupby("fiscal_year"):
109 + g = g.dropna(subset=[z]).copy()
110 + if len(g) < minn: continue
111 + g["q"] = pd.qcut(g[z].rank(method="first"), n, labels=False) + 1; hold = pd.date_range(pd.Timestamp(yr + 1, 7, 31), pd.Timestamp(yr + 2, 6, 30), freq="ME")
112 + for q in (1, n):
113 + tks = [t for t in g[g.q == q].ticker if t in ret_by]
114 + for dt in hold:
115 + r = [ret_by[t].get(dt, np.nan) for t in tks]; r = [x for x in r if pd.notna(x)]
116 + if len(r) >= 10: monthly.setdefault(dt, {}).setdefault(q, []).append(np.mean(r))
117 + recs = [{"date": dt, "lo": np.mean(v[1]), "hi": np.mean(v[n])} for dt, v in sorted(monthly.items()) if 1 in v and n in v]; ts = pd.DataFrame(recs); ts["ls"] = ts["hi"] - ts["lo"]; return ts
118 +SR_TRIALS = []
119 +def factor_alpha(ts, record=True):
120 + s = ts.merge(ff_ix.reset_index(), on="date", how="inner").dropna(subset=["ls", "MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]); y = s["ls"]; out = {}
121 + if len(s) < 12: return {"n_months": int(len(s))}
122 + r1 = sm.OLS(y, sm.add_constant(s[["MktRF"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 3}); r2 = sm.OLS(y, sm.add_constant(s[["MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 3})
123 + mu, t = I.nw_t(y.values); sr = mu / y.std() * np.sqrt(12)
124 + out.update(mean_ls_ann_pct=round(mu * 1200, 2), hac_t=round(t, 2), plain_t=round(I.plain_t(y.values), 2), sharpe=round(sr, 2), skew=round(float(y.skew()), 3), kurt=round(float(y.kurt() + 3), 3), capm_alpha_ann_pct=round(r1.params["const"] * 1200, 2), capm_alpha_t=round(r1.tvalues["const"], 2),
125 + ff6_alpha_ann_pct=round(r2.params["const"] * 1200, 2), ff6_alpha_t=round(r2.tvalues["const"], 2), n_months=int(len(s)), betas={k: round(r2.params[k], 3) for k in ["MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]})
126 + if record: SR_TRIALS.append(sr)
127 + return out
128 +RES["ls_alpha"] = {}; POST_START = pd.Timestamp("2023-07-31")
129 +for k in ["AI_ADOPT", "AI_INFRA", "AI_COMPOS", "AI_RISK", "AI_PC1", "AI_formw", "AI_10k", "SPID_TECH", "CYBER"]:
130 + ts = ls_monthly(df, f"z_{k}"); ts.to_parquet(f"{V}/ls_ts_v2_{k}.parquet", index=False)
131 + RES["ls_alpha"][k] = {"full": factor_alpha(ts), "pre": factor_alpha(ts[ts.date < POST_START]), "post": factor_alpha(ts[ts.date >= POST_START])}
132 + print(k, "LS full", RES["ls_alpha"][k]["full"].get("mean_ls_ann_pct"), "ff6a", RES["ls_alpha"][k]["full"].get("ff6_alpha_ann_pct"), "(t", RES["ls_alpha"][k]["full"].get("ff6_alpha_t"), ") | post ff6a", RES["ls_alpha"][k]["post"].get("ff6_alpha_ann_pct"), "(t", RES["ls_alpha"][k]["post"].get("ff6_alpha_t"), ")", flush=True)
133 +
134 +# ---------------- two-way FE panel (fwd_w within-year) ----------------
135 +pdf = df.dropna(subset=["fwd_w", "z_AI_ADOPT", "z_AI_INFRA", "z_AI_RISK", "z_size", "z_mom", "z_vol"]).set_index(["ticker", "fiscal_year"])
136 +specs = {"M1_adopt": "fwd_w ~ z_AI_ADOPT", "M2_adopt_ctrl": "fwd_w ~ z_AI_ADOPT + z_size + z_mom + z_vol", "M3_infra_ctrl": "fwd_w ~ z_AI_INFRA + z_size + z_mom + z_vol", "M4_risk_ctrl": "fwd_w ~ z_AI_RISK + z_size + z_mom + z_vol",
137 + "M5_all": "fwd_w ~ z_AI_ADOPT + z_AI_INFRA + z_AI_RISK + z_size + z_mom + z_vol", "M6_spid": "fwd_w ~ z_SPID_TECH + z_size + z_mom + z_vol", "M7_all_spid": "fwd_w ~ z_AI_COMPOS + z_SPID_TECH + z_size + z_mom + z_vol", "M8_compos": "fwd_w ~ z_AI_COMPOS + z_size + z_mom + z_vol"}
138 +RES["panel_fe"] = {}
139 +for lab, f in specs.items():
140 + r = PanelOLS.from_formula(f + " + EntityEffects + TimeEffects", data=pdf, drop_absorbed=True).fit(cov_type="clustered", cluster_entity=True, cluster_time=True)
141 + RES["panel_fe"][lab] = {"coefs": {v: {"coef": round(r.params[v], 4), "t": round(r.tstats[v], 2)} for v in r.params.index}, "n": int(r.nobs), "r2_within": round(r.rsquared_within, 4)}
142 +print("FE M8", RES["panel_fe"]["M8_compos"]["coefs"]["z_AI_COMPOS"], flush=True)
143 +
144 +# ---------------- Fama-MacBeth: raw | within-year winsorized | pooled-winsorized | strict ----------------
145 +FB15 = I.fixed_b_cv(15, 3)
146 +def fama_macbeth(d, xvars, minn=40, oc="fwd_ret"):
147 + coefs = {v: [] for v in ["const"] + xvars}; ns = []
148 + for yr, g in d.groupby("fiscal_year"):
149 + g = g.dropna(subset=[oc] + xvars)
150 + if len(g) < minn: continue
151 + r = sm.OLS(g[oc], sm.add_constant(g[xvars])).fit(); ns.append(int(r.nobs))
152 + for v in coefs: coefs[v].append(r.params.get(v, np.nan))
153 + out = {}
154 + for v in coefs:
155 + mu, t = I.nw_t(coefs[v]); out[v] = {"mean_coef": round(mu, 4), "t": round(t, 2), "t_plain": round(I.plain_t(coefs[v]), 2), "ci95": [round(x, 4) for x in I.ci95(coefs[v])], "ci95_boot": [round(x, 4) for x in I.block_boot_ci(coefs[v])], "fixed_b_cv05": FB15["cv"], "by_year": [round(x, 4) for x in coefs[v]]}
156 + out["n_years"] = len(coefs["const"]); out["n_mean_per_year"] = int(np.mean(ns)) if ns else 0; return out
157 +RES["fama_macbeth"] = {}
158 +for sfx, oc in [("", "fwd_ret"), ("_w", "fwd_w"), ("_wp", "fwd_w_pooled"), ("_strict", "fwd_ret_12m_strict")]:
159 + RES["fama_macbeth"]["compos" + sfx] = fama_macbeth(df, ["z_AI_COMPOS", "z_size", "z_mom", "z_vol"], oc=oc)
160 + if sfx in ("", "_w"):
161 + for lab, xv in [("adopt", ["z_AI_ADOPT"]), ("infra", ["z_AI_INFRA"]), ("risk", ["z_AI_RISK"]), ("all", ["z_AI_ADOPT", "z_AI_INFRA", "z_AI_RISK"]), ("pc1", ["z_AI_PC1"]), ("formw", ["z_AI_formw"]), ("tenk", ["z_AI_10k"]), ("spid", ["z_SPID_TECH"]), ("compos_spid", ["z_AI_COMPOS", "z_SPID_TECH"]), ("cyber", ["z_CYBER"]), ("compos_cyber", ["z_AI_COMPOS", "z_CYBER"]), ("spid_any_ai", ["spid_ai_any"]), ("spid_any_tech", ["spid_tech_any"]), ("compos_spid_any", ["z_AI_COMPOS", "spid_ai_any"]), ("compos_nocontrols", [])]:
162 + RES["fama_macbeth"][lab + sfx] = fama_macbeth(df, (xv if lab != "compos_nocontrols" else ["z_AI_COMPOS"]) + (["z_size", "z_mom", "z_vol"] if lab != "compos_nocontrols" else []), oc=oc)
163 +for k in ["compos", "compos_w", "compos_wp", "compos_strict", "all_w", "pc1_w", "formw_w", "tenk_w"]:
164 + v = RES["fama_macbeth"][k]; key = [x for x in v if x.startswith("z_AI")][0]; print("FM", k, key, v[key]["mean_coef"], "NW t", v[key]["t"], "plain t", v[key]["t_plain"], "CI", v[key]["ci95"], flush=True)
165 +
166 +# ---------------- structural break regressions: firm / year / two-way / wild bootstrap ----------------
167 +def break_reg(d, z, oc="fwd_w", boot=True):
168 + d2 = d.dropna(subset=[oc, z, "z_size", "z_mom"]).copy(); d2["zp"] = d2[z] * d2["post"]
169 + X = sm.add_constant(d2[[z, "post", "zp", "z_size", "z_mom"]].values); y = d2[oc].values; j = 3
170 + r = sm.OLS(y, X).fit(); gf = pd.factorize(d2.ticker)[0]; gy = pd.factorize(d2.fiscal_year)[0]
171 + se_f = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params())); se_y = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gy).cov_params())); se_tw = np.sqrt(np.diag(I.cgm_cov(r, gf, gy)))
172 + out = {"pre_slope": round(float(r.params[1]), 4), "interaction": round(float(r.params[j]), 4), "post_slope": round(float(r.params[1] + r.params[j]), 4), "interaction_t_firm": round(float(r.params[j] / se_f[j]), 2), "interaction_t_year": round(float(r.params[j] / se_y[j]), 2), "interaction_t_twoway": round(float(r.params[j] / se_tw[j]), 2),
173 + "pre_slope_t_firm": round(float(r.params[1] / se_f[1]), 2), "pre_slope_t_twoway": round(float(r.params[1] / se_tw[1]), 2), "interaction_ci95_firm": [round(float(r.params[j] - 1.96 * se_f[j]), 4), round(float(r.params[j] + 1.96 * se_f[j]), 4)], "n": int(r.nobs), "G_years": int(gy.max() + 1)}
174 + if boot:
175 + for wname in ["rademacher", "webb"]:
176 + b = I.wild_cluster_boot(y, X, gy, j, B=9999, weights=wname); out[f"interaction_p_wild_{wname}"] = b["p"]
177 + return out
178 +RES["structural_break"] = {}
179 +for k in ["AI_ADOPT", "AI_INFRA", "AI_COMPOS", "AI_RISK", "AI_PC1", "AI_formw", "SPID_TECH"]:
180 + RES["structural_break"][k] = break_reg(df, f"z_{k}"); RES["structural_break"][k + "_raw"] = break_reg(df, f"z_{k}", oc="fwd_ret", boot=False)
181 + print("break", k, {x: RES["structural_break"][k][x] for x in ["interaction", "interaction_t_firm", "interaction_t_twoway", "interaction_p_wild_rademacher", "interaction_p_wild_webb"]}, flush=True)
182 +
183 +# ---------------- sub-period spreads (raw | winsorized | median), n beside every t ----------------
184 +def spread_sub(d, z, yrs, oc="fwd_ret", stat="mean", minn=40):
185 + s = sort_table(d, z, None, oc, stat, minn, yrs=set(yrs)); return {"mean": round(s["Q5_minus_Q1"], 4) if s["n_years"] else None, "t": round(s["t_stat"], 2) if s["t_stat"] == s["t_stat"] else None, "t_plain": round(s["t_plain"], 2) if s["t_plain"] == s["t_plain"] else None, "ci95": s["ci95"], "ci95_boot": s["ci95_boot"], "n": s["n_years"], "by_cohort": s["by_cohort"], "top1pct_contribution": round(s["top1pct_contribution"], 4) if s["n_years"] else None}
186 +RES["subperiod_spread"] = {}
187 +for k in SIGL:
188 + RES["subperiod_spread"][k] = {}
189 + for sfx, oc, st in [("", "fwd_ret", "mean"), ("_w", "fwd_w", "mean"), ("_med", "fwd_ret", "median"), ("_strict", "fwd_ret_12m_strict", "mean")]:
190 + RES["subperiod_spread"][k]["pre" + sfx] = spread_sub(df, f"z_{k}", PRE, oc, st); RES["subperiod_spread"][k]["post" + sfx] = spread_sub(df, f"z_{k}", POST, oc, st)
191 +v = RES["subperiod_spread"]["AI_COMPOS"]; print("COMPOS pre/post raw", v["pre"]["mean"], v["post"]["mean"], "| wins", v["pre_w"]["mean"], v["post_w"]["mean"], "| median", v["pre_med"]["mean"], v["post_med"]["mean"], "| post cohorts wins", v["post_w"]["by_cohort"], flush=True)
192 +
193 +# ---------------- pooled Q1/Q5 regressions with firm / two-way / wild inference [E] ----------------
194 +def pooled(d, yrs, dv="fwd_ret", z="z_AI_COMPOS", boot=True):
195 + d = d[d.fiscal_year.isin(list(yrs))].dropna(subset=[z, dv]).copy(); d["q"] = d.groupby("fiscal_year")[z].transform(lambda s: pd.qcut(s.rank(method="first"), 5, labels=False) + 1)
196 + d = d[d.q.isin([1, 5])].copy(); d["hi"] = (d.q == 5).astype(float); yd = pd.get_dummies(d.fiscal_year.astype(int), prefix="y", drop_first=True).astype(float)
197 + X = np.c_[np.ones(len(d)), d["hi"].values, yd.values]; y = d[dv].values; r = sm.OLS(y, X).fit(); gf = pd.factorize(d.ticker)[0]; gy = pd.factorize(d.fiscal_year)[0]
198 + se_f = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params())); G = gy.max() + 1
199 + out = {"coef_pct": round(float(r.params[1]) * 100, 2), "t_firm": round(float(r.params[1] / se_f[1]), 2), "ci95_firm_pct": [round(float(r.params[1] - 1.96 * se_f[1]) * 100, 2), round(float(r.params[1] + 1.96 * se_f[1]) * 100, 2)], "n": int(r.nobs), "G_years": int(G)}
200 + if G >= 2:
201 + se_tw = np.sqrt(np.diag(I.cgm_cov(r, gf, gy))); se_y = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gy).cov_params())); out["t_twoway"] = round(float(r.params[1] / se_tw[1]), 2); out["t_year"] = round(float(r.params[1] / se_y[1]), 2)
202 + if boot:
203 + for wname in ["rademacher", "webb"]: out[f"p_wild_{wname}"] = I.wild_cluster_boot(y, X, gy, 1, B=9999, weights=wname)["p"]
204 + return out, d
205 +RES["pooled"] = {}
206 +for sfx, dv in [("", "fwd_ret"), ("_w", "fwd_w"), ("_med", None)]:
207 + if dv is None: continue
208 + RES["pooled"]["post" + sfx], dpost = pooled(df, POST, dv); RES["pooled"]["pre" + sfx], dpre = pooled(df, PRE, dv); RES["pooled"]["post_excl_2024" + sfx], _ = pooled(df, range(2022, 2024), dv); RES["pooled"]["post_2023_2024" + sfx], _ = pooled(df, range(2023, 2025), dv)
209 + dall = pd.concat([dpre.assign(post=0), dpost.assign(post=1)]); dall["hp"] = dall.hi * dall.post; yd = pd.get_dummies(dall.fiscal_year.astype(int), prefix="y", drop_first=True).astype(float)
210 + X = np.c_[np.ones(len(dall)), dall["hi"].values, dall["hp"].values, yd.values]; y = dall[dv].values; r = sm.OLS(y, X).fit(); gf = pd.factorize(dall.ticker)[0]; gy = pd.factorize(dall.fiscal_year)[0]
211 + se_f = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params())); se_tw = np.sqrt(np.diag(I.cgm_cov(r, gf, gy)))
212 + RES["pooled"]["diff" + sfx] = {"coef_pct": round(float(r.params[2]) * 100, 2), "t_firm": round(float(r.params[2] / se_f[2]), 2), "t_twoway": round(float(r.params[2] / se_tw[2]), 2), "ci95_firm_pct": [round(float(r.params[2] - 1.96 * se_f[2]) * 100, 2), round(float(r.params[2] + 1.96 * se_f[2]) * 100, 2)], "n": int(r.nobs),
213 + "p_wild_rademacher": I.wild_cluster_boot(y, X, gy, 2, B=9999, weights="rademacher")["p"], "p_wild_webb": I.wild_cluster_boot(y, X, gy, 2, B=9999, weights="webb")["p"], "mde_pct_firm": round(I.mde(se_f[2]) * 100, 2)}
214 +print("pooled post", RES["pooled"]["post"], "\npooled post wins", RES["pooled"]["post_w"], "\ndiff", RES["pooled"]["diff"], flush=True)
215 +
216 +# ---------------- size groups (both partitions) and index membership ----------------
217 +RES["groups"] = {}
218 +def group_block(d, gname):
219 + blk = {"n": int(len(d)), "firms": int(d.ticker.nunique())}
220 + for sfx, oc, st in [("", "fwd_ret", "mean"), ("_w", "fwd_w", "mean"), ("_med", "fwd_ret", "median")]:
221 + blk["sort_EW" + sfx] = sort_table(d, "z_AI_COMPOS", None, oc, st); blk["pre" + sfx] = spread_sub(d, "z_AI_COMPOS", PRE, oc, st); blk["post" + sfx] = spread_sub(d, "z_AI_COMPOS", POST, oc, st)
222 + cs = blk["sort_EW" + sfx]["by_cohort"]; post = {k: v for k, v in cs.items() if k >= 2022}
223 + blk["post_loo" + sfx] = {k: round(float(np.mean([v for kk, v in post.items() if kk != k])), 4) for k in post}
224 + blk["pooled_post" + sfx], _ = pooled(d, POST, oc, boot=False) if st == "mean" else (None, None); blk["pooled_post_excl_2024" + sfx], _ = pooled(d, range(2022, 2024), oc, boot=False) if st == "mean" else (None, None)
225 + blk["fm_compos"] = fama_macbeth(d, ["z_AI_COMPOS", "z_size", "z_mom", "z_vol"]); blk["fm_compos_w"] = fama_macbeth(d, ["z_AI_COMPOS", "z_size", "z_mom", "z_vol"], oc="fwd_w")
226 + ts = ls_monthly(d, "z_AI_COMPOS"); ts.to_parquet(f"{V}/ls_ts_v2_AI_COMPOS_{gname}.parquet", index=False); blk["ls"] = {"full": factor_alpha(ts), "pre": factor_alpha(ts[ts.date < POST_START]), "post": factor_alpha(ts[ts.date >= POST_START])}
227 + return blk
228 +for gname, mask in [("large_top1000", df.size_group == "large"), ("small_rest", df.size_group == "small"), ("large_pct42", df.size_group_pct42 == "large"), ("small_pct42", df.size_group_pct42 == "small"), ("large_top1000_v1", df.size_group_v1 == "large"), ("small_rest_v1", df.size_group_v1 == "small"), ("sp500", df.in_sp500 == True), ("non_sp500", df.in_sp500 == False)]:
229 + RES["groups"][gname] = group_block(df[mask], gname); g = RES["groups"][gname]
230 + print(gname, "n", g["n"], "full", round(g["sort_EW"]["Q5_minus_Q1"], 4), "| pre/post raw", g["pre"]["mean"], g["post"]["mean"], "| wins", g["pre_w"]["mean"], g["post_w"]["mean"], "| LOO wins", g["post_loo_w"], "| ff6 post", g["ls"]["post"].get("ff6_alpha_ann_pct"), flush=True)
231 +df["z_AI_COMPOS_sizeneu"] = df.groupby(["fiscal_year", "size_group"])["z_AI_COMPOS"].transform(lambda g: zscore(g) if g.std(ddof=0) > 0 else g * 0)
232 +RES["size_neutral"] = {"full": sort_table(df, "z_AI_COMPOS_sizeneu"), "full_w": sort_table(df, "z_AI_COMPOS_sizeneu", None, "fwd_w"), "pre": spread_sub(df, "z_AI_COMPOS_sizeneu", PRE), "post": spread_sub(df, "z_AI_COMPOS_sizeneu", POST), "pre_w": spread_sub(df, "z_AI_COMPOS_sizeneu", PRE, "fwd_w"), "post_w": spread_sub(df, "z_AI_COMPOS_sizeneu", POST, "fwd_w")}
233 +
234 +# ---------------- contrasts: hi x post x large, hi x post x member [E] ----------------
235 +def contrast(d, grp, dv="fwd_w"):
236 + d = d.dropna(subset=["z_AI_COMPOS", dv, grp]).copy(); d["q"] = d.groupby("fiscal_year")["z_AI_COMPOS"].transform(lambda s: pd.qcut(s.rank(method="first"), 5, labels=False) + 1); d = d[d.q.isin([1, 5])].copy()
237 + d["hi"] = (d.q == 5).astype(float); d["g"] = d[grp].astype(float); d["hp"] = d.hi * d.post; d["hg"] = d.hi * d.g; d["pg"] = d.post * d.g; d["hpg"] = d.hi * d.post * d.g
238 + yd = pd.get_dummies(d.fiscal_year.astype(int), prefix="y", drop_first=True).astype(float); X = np.c_[np.ones(len(d)), d[["hi", "hp", "g", "hg", "pg", "hpg"]].values, yd.values]; y = d[dv].values
239 + r = sm.OLS(y, X).fit(); gf = pd.factorize(d.ticker)[0]; gy = pd.factorize(d.fiscal_year)[0]; se_f = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params())); se_tw = np.sqrt(np.diag(I.cgm_cov(r, gf, gy)))
240 + out = {}
241 + for j, nm in [(2, "hi_x_post"), (4, "hi_x_group"), (6, "hi_x_post_x_group")]:
242 + out[nm] = {"coef_pct": round(float(r.params[j]) * 100, 2), "t_firm": round(float(r.params[j] / se_f[j]), 2), "t_twoway": round(float(r.params[j] / se_tw[j]), 2), "ci95_firm_pct": [round(float(r.params[j] - 1.96 * se_f[j]) * 100, 2), round(float(r.params[j] + 1.96 * se_f[j]) * 100, 2)], "mde_pct_firm": round(I.mde(se_f[j]) * 100, 2)}
243 + out["hi_x_post_x_group"]["p_wild_rademacher"] = I.wild_cluster_boot(y, X, gy, 6, B=9999, weights="rademacher")["p"]; out["hi_x_post_x_group"]["p_wild_webb"] = I.wild_cluster_boot(y, X, gy, 6, B=9999, weights="webb")["p"]; out["n"] = int(r.nobs); return out
244 +df["is_large"] = (df.size_group == "large").astype(float).where(df.size_group.notna()); df["is_member"] = df.in_sp500.astype(float)
245 +RES["contrasts"] = {"large_vs_small_w": contrast(df, "is_large"), "large_vs_small_raw": contrast(df, "is_large", "fwd_ret"), "member_vs_non_w": contrast(df, "is_member"), "member_vs_non_raw": contrast(df, "is_member", "fwd_ret")}
246 +# cohort-block bootstrap of the difference in (post - pre) spreads between groups
247 +def group_diff_boot(ga, gb, sfx="_w", B=9999):
248 + a = RES["groups"][ga]["sort_EW" + sfx]["by_cohort"]; b = RES["groups"][gb]["sort_EW" + sfx]["by_cohort"]; yrs = sorted(set(a) & set(b)); pre = [y for y in yrs if y < 2022]; post = [y for y in yrs if y >= 2022]
249 + def stat(pre_i, post_i): return (np.mean([a[y] for y in post_i]) - np.mean([a[y] for y in pre_i])) - (np.mean([b[y] for y in post_i]) - np.mean([b[y] for y in pre_i]))
250 + obs = stat(pre, post); rng = np.random.default_rng(I.RNG_SEED); bs = [stat(rng.choice(pre, len(pre)), rng.choice(post, len(post))) for _ in range(B)]
251 + return {"diff_in_diff_pct": round(obs * 100, 2), "ci95_boot_pct": [round(float(np.quantile(bs, .025)) * 100, 2), round(float(np.quantile(bs, .975)) * 100, 2)], "n_pre": len(pre), "n_post": len(post)}
252 +RES["contrasts"]["large_minus_small_postpre_boot_w"] = group_diff_boot("large_top1000", "small_rest"); RES["contrasts"]["member_minus_non_postpre_boot_w"] = group_diff_boot("sp500", "non_sp500"); RES["contrasts"]["large_minus_small_postpre_boot_raw"] = group_diff_boot("large_top1000", "small_rest", "")
253 +print("contrast large", RES["contrasts"]["large_vs_small_w"]["hi_x_post_x_group"], "| member", RES["contrasts"]["member_vs_non_w"]["hi_x_post_x_group"], flush=True)
254 +
255 +# ---------------- SPID binary contrasts, CI, MDE [I] ----------------
256 +def binary_contrast(d, flag, yrs, oc="fwd_ret", min_treated=15):
257 + diffs = {}; ns = {}
258 + for yr, g in d[d.fiscal_year.isin(list(yrs))].groupby("fiscal_year"):
259 + a = g[g[flag] == 1][oc].dropna(); b = g[g[flag] == 0][oc].dropna()
260 + if len(a) >= min_treated and len(b) >= min_treated: diffs[int(yr)] = float(a.mean() - b.mean()); ns[int(yr)] = int(len(a))
261 + ls = list(diffs.values()); mu, t = I.nw_t(ls); se = (np.std(ls, ddof=1) / np.sqrt(len(ls))) if len(ls) > 1 else np.nan
262 + return {"mean_pct": round(mu * 100, 2) if ls else None, "t": round(t, 2) if t == t else None, "t_plain": round(I.plain_t(ls), 2) if len(ls) > 1 else None, "ci95_pct": [round(x * 100, 2) for x in I.ci95(ls)] if len(ls) > 1 else [None, None], "n_cohorts": len(ls), "n_treated_by_cohort": ns, "mde_pct": round(I.mde(se) * 100, 2) if se == se else None, "by_cohort_pct": {k: round(v * 100, 2) for k, v in diffs.items()}}
263 +def firm_level_contrast(d, flag, yrs, oc="fwd_w"):
264 + d = d[d.fiscal_year.isin(list(yrs))].dropna(subset=[oc, flag, "z_size", "z_mom"]).copy(); yd = pd.get_dummies(d.fiscal_year.astype(int), prefix="y", drop_first=True).astype(float)
265 + X = np.c_[np.ones(len(d)), d[flag].values, d[["z_size", "z_mom"]].values, yd.values]; y = d[oc].values; r = sm.OLS(y, X).fit(); gf = pd.factorize(d.ticker)[0]; se = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params()))
266 + return {"coef_pct": round(float(r.params[1]) * 100, 2), "t_firm": round(float(r.params[1] / se[1]), 2), "ci95_pct": [round(float(r.params[1] - 1.96 * se[1]) * 100, 2), round(float(r.params[1] + 1.96 * se[1]) * 100, 2)], "mde_pct": round(I.mde(se[1]) * 100, 2), "n": int(r.nobs), "n_treated": int(d[flag].sum())}
267 +RES["spid"] = {"share_firm_years_any_ai": round(float(df.spid_ai_any.mean()), 4), "share_firm_years_any_tech": round(float(df.spid_tech_any.mean()), 4), "share_post_any_ai": round(float(df[df.post == 1].spid_ai_any.mean()), 4), "any_ai_by_year": df.groupby("fiscal_year").spid_ai_any.mean().round(4).to_dict(), "n_any_ai_by_year": df.groupby("fiscal_year").spid_ai_any.sum().astype(int).to_dict(),
268 + "spid_tech_quantiles": {q: float(df.spid_mentions_tech.quantile(q)) for q in [.5, .75, .9, .95]}}
269 +for flag in ["spid_ai_any", "spid_tech_any"]:
270 + RES["spid"][flag] = {}
271 + for per, yrs in [("full", range(2010, 2025)), ("pre", PRE), ("post", POST)]:
272 + RES["spid"][flag][per] = {"cohort_raw": binary_contrast(df, flag, yrs), "cohort_w": binary_contrast(df, flag, yrs, "fwd_w"), "firm_w": firm_level_contrast(df, flag, yrs), "firm_raw": firm_level_contrast(df, flag, yrs, "fwd_ret")}
273 +# terciles among non-zero tech-project firms (>= 15 treated per cohort)
274 +nz = df[df.spid_mentions_tech > 0].copy(); nz["spid_terc"] = nz.groupby("fiscal_year")["spid_mentions_tech"].transform(lambda s: pd.qcut(s.rank(method="first"), 3, labels=False) + 1 if len(s) >= 45 else np.nan)
275 +tsp = {}
276 +for yr, g in nz.dropna(subset=["spid_terc"]).groupby("fiscal_year"):
277 + a = g[g.spid_terc == 3].fwd_w; b = g[g.spid_terc == 1].fwd_w
278 + if len(a) >= 15 and len(b) >= 15: tsp[int(yr)] = float(a.mean() - b.mean())
279 +mu, t = I.nw_t(list(tsp.values())); RES["spid"]["tercile_T3_minus_T1_nonzero_w"] = {"mean_pct": round(mu * 100, 2), "t": round(t, 2) if t == t else None, "t_plain": round(I.plain_t(list(tsp.values())), 2), "ci95_pct": [round(x * 100, 2) for x in I.ci95(list(tsp.values()))], "n_cohorts": len(tsp), "by_cohort_pct": {k: round(v * 100, 2) for k, v in tsp.items()}}
280 +# legacy quintile spread (SPID_TECH) kept only for the before/after, with its MDE
281 +sp_leg = RES["subperiod_spread"]["SPID_TECH"]; ls_leg = list(RES["sorts"]["SPID_TECH_EW"]["by_cohort"].values()); RES["spid"]["legacy_quintile_full_mde_pct"] = round(I.mde(np.std(ls_leg, ddof=1) / np.sqrt(len(ls_leg))) * 100, 2)
282 +print("SPID any-AI post cohort", RES["spid"]["spid_ai_any"]["post"]["cohort_raw"], "\n firm", RES["spid"]["spid_ai_any"]["post"]["firm_w"], flush=True)
283 +
284 +# ---------------- side-by-side samples: Russell 3000 vs S&P 500 members ----------------
285 +RES["samples"] = {}
286 +for sname, d in [("r3k", df), ("sp500", df[df.in_sp500 == True].copy())]:
287 + for k, col in SIGNALS.items(): d[f"zz_{k}"] = d.groupby("fiscal_year")[col].transform(zscore)
288 + d["zz_AI_COMPOS"] = d[["zz_AI_ADOPT", "zz_AI_INFRA"]].mean(axis=1); d["zz_SPID_TECH"] = d.groupby("fiscal_year")["spid_tech_log"].transform(zscore); d["zz_AI_PC1"] = d.groupby("fiscal_year")["z_AI_PC1"].transform(zscore); d["zz_AI_formw"] = d.groupby("fiscal_year")["z_AI_formw"].transform(zscore)
289 + for c in ["size", "mom", "vol"]: d[f"zz_{c}"] = d.groupby("fiscal_year")[{"size": "logcap_w", "mom": "mom_w", "vol": "vol_w"}[c]].transform(zscore)
290 + for c, w in [("fwd_ret", "fwd_w_s"), ("fwd_ret_mktadj", "fwd_madj_w_s")]: d[w] = d.groupby("fiscal_year")[c].transform(wins_within)
291 + blk = {"n": int(len(d)), "firms": int(d.ticker.nunique()), "sorts": {}, "ls": {}, "fm": {}, "break": {}, "sub": {}}
292 + for k in ["AI_ADOPT", "AI_INFRA", "AI_RISK", "AI_COMPOS", "AI_PC1", "AI_formw", "SPID_TECH"]:
293 + for sfx, oc, st in [("", "fwd_ret", "mean"), ("_w", "fwd_w_s", "mean"), ("_med", "fwd_ret", "median")]:
294 + blk["sorts"][f"{k}_EW{sfx}"] = sort_table(d, f"zz_{k}", None, oc, st)
295 + if sfx != "_med": blk["sorts"][f"{k}_VW{sfx}"] = sort_table(d, f"zz_{k}", "mktcap", oc, st, diag_label=f"{sname}_{k}_VW{sfx}")
296 + blk["sub"].setdefault(k, {})["pre" + sfx] = spread_sub(d, f"zz_{k}", PRE, oc, st); blk["sub"][k]["post" + sfx] = spread_sub(d, f"zz_{k}", POST, oc, st)
297 + ts = ls_monthly(d, f"zz_{k}"); ts.to_parquet(f"{V}/ls_ts_v2_{k}_{sname}.parquet", index=False); blk["ls"][k] = {"full": factor_alpha(ts), "pre": factor_alpha(ts[ts.date < POST_START]), "post": factor_alpha(ts[ts.date >= POST_START])}
298 + if k not in ("AI_RISK",): blk["break"][k] = break_reg(d.rename(columns={f"zz_{k}": f"zzz_{k}", "zz_size": "z_size_s", "zz_mom": "z_mom_s"}).assign(z_size=lambda x: x.z_size_s, z_mom=lambda x: x.z_mom_s), f"zzz_{k}", oc="fwd_w_s", boot=(k == "AI_COMPOS"))
299 + for sfx, oc in [("", "fwd_ret"), ("_w", "fwd_w_s")]:
300 + blk["fm"]["compos" + sfx] = fama_macbeth(d, ["zz_AI_COMPOS", "zz_size", "zz_mom", "zz_vol"], oc=oc); blk["fm"]["all" + sfx] = fama_macbeth(d, ["zz_AI_ADOPT", "zz_AI_INFRA", "zz_AI_RISK", "zz_size", "zz_mom", "zz_vol"], oc=oc); blk["fm"]["spid" + sfx] = fama_macbeth(d, ["zz_SPID_TECH", "zz_size", "zz_mom", "zz_vol"], oc=oc); blk["fm"]["pc1" + sfx] = fama_macbeth(d, ["zz_AI_PC1", "zz_size", "zz_mom", "zz_vol"], oc=oc); blk["fm"]["formw" + sfx] = fama_macbeth(d, ["zz_AI_formw", "zz_size", "zz_mom", "zz_vol"], oc=oc)
301 + blk["cohorts_pct"] = {k_: round(v * 100, 1) for k_, v in blk["sorts"]["AI_COMPOS_EW"]["by_cohort"].items()}; blk["cohorts_w_pct"] = {k_: round(v * 100, 1) for k_, v in blk["sorts"]["AI_COMPOS_EW_w"]["by_cohort"].items()}
302 + blk["pooled_post"], _ = pooled(d, POST, "fwd_ret", "zz_AI_COMPOS"); blk["pooled_post_w"], _ = pooled(d, POST, "fwd_w_s", "zz_AI_COMPOS"); blk["pooled_post_excl_2024_w"], _ = pooled(d, range(2022, 2024), "fwd_w_s", "zz_AI_COMPOS", boot=False)
303 + RES["samples"][sname] = blk
304 + print("SAMPLE", sname, "n", blk["n"], "compos EW raw", round(blk["sorts"]["AI_COMPOS_EW"]["Q5_minus_Q1"], 4), "wins", round(blk["sorts"]["AI_COMPOS_EW_w"]["Q5_minus_Q1"], 4), "(NW t", round(blk["sorts"]["AI_COMPOS_EW_w"]["t_stat"], 2), ") pre/post wins", blk["sub"]["AI_COMPOS"]["pre_w"]["mean"], blk["sub"]["AI_COMPOS"]["post_w"]["mean"], "| pooled post wins", blk["pooled_post_w"], flush=True)
305 +pd.DataFrame(VWD).to_csv(f"{T}/vw_weights_diag.csv", index=False)
306 +RES["sr_trials_ann"] = [round(float(x), 4) for x in SR_TRIALS]; RES["runtime_s"] = round(time.time() - T0)
307 +json.dump(RES, open(f"{T}/results.json", "w"), indent=2, default=str); df.to_parquet(f"{V}/analysis_panel_v2.parquet", index=False); print("saved results/r3k_v2/results.json in", RES["runtime_s"], "s", flush=True)
added scripts/r3k_16_extended_v2.py +370 −0
@@ -0,0 +1,370 @@
1 +"""r3k_16_extended_v2.py — Stage D (FIRC-0002): extended tests on the rebuilt panel (v2 of r3k_04_extended.py; idempotent: reads
2 +analysis_panel_v2.parquet written by r3k_15 and never re-merges into it).
3 +Blocks (tags = review clusters): A placebo factors (+ two-way / wild-bootstrap inference, MHT family) · B sector-neutral, double sorts
4 +(both size partitions), residual sort, FM +/- controls with and without Energy [H] · C mega-cap exclusion incl. the ex-top-500 PRE
5 +cell [O] · D alt signals / horizons / benchmark (raw & winsorized) · E year-by-year FM (raw & winsorized) · G DiD / event-time with
6 +the main effect restored (differences from 2021 with firm-clustered CIs), trailing-exposure and cohorts >= 2019 variants [G] ·
7 +H forward vol · I cohorts, LOO, pooled tests (from r3k_15) · J alternative splits with pooled firm-clustered tests [M] · K entrants
8 +· L sector controls with plain / NW / fixed-b t and two-way + wild inference [E, H] · M TECH spanning · N filing-type composition,
9 +HAC Wald test of form-slope equality, post-launch share of 2022-cohort chunks and the 2022-cohort re-run on pre-launch filings [H, M]
10 +· O placebo break dates + exact permutation p + sup-Wald [M] · P SPID interactions (explicit terms) [C, I] · Q reversal and cyber-
11 +neutral tests [K] · R multiple-testing families (Holm / BH / Romano-Wolf on shared bootstrap draws), deflated Sharpe [N] ·
12 +S specification curve [N] · coverage.
13 +Writes results/r3k_v2/results2.json (+ CSVs). Run: OMP_NUM_THREADS=4 nice -n 10 .venv-review/bin/python -I scripts/r3k_16_extended_v2.py
14 +"""
15 +import os, sys, json, time, itertools, warnings, numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf
16 +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))); import r3k_inference as I
17 +from linearmodels.panel import PanelOLS
18 +warnings.filterwarnings("ignore"); T0 = time.time()
19 +D = "data/processed_r3k"; V = "data/processed_r3k_v2"; R = "data/raw_russell3000"; T = "results/r3k_v2"
20 +df = pd.read_parquet(f"{V}/analysis_panel_v2.parquet"); RES = json.load(open(f"{T}/results.json")); mret = pd.read_parquet(f"{V}/monthly_returns_v2.parquet")
21 +mret["date"] = pd.to_datetime(mret["date"]).dt.to_period("M").dt.to_timestamp("M")
22 +ff = pd.read_parquet(f"{D}/ff5_monthly.parquet").merge(pd.read_parquet(f"{D}/mom_monthly.parquet"), on="date", how="left"); ff["date"] = pd.to_datetime(ff["date"]).dt.to_period("M").dt.to_timestamp("M"); ff = ff.set_index("date")
23 +R2 = {"script": "r3k_16_extended_v2.py"}; PRE, POST, ALL = range(2010, 2022), range(2022, 2025), range(2010, 2025); POST_START = pd.Timestamp("2023-07-31")
24 +def zc(g): return (g - g.mean()) / g.std(ddof=0)
25 +def wins_within(s): return s.clip(s.quantile(.01), s.quantile(.99))
26 +def cohort_spreads(d, z, yrs, oc="fwd_ret", stat="mean", minn=40):
27 + out = {}
28 + for yr, g in d[d.fiscal_year.isin(list(yrs))].groupby("fiscal_year"):
29 + g = g.dropna(subset=[z, oc]).copy()
30 + if len(g) < minn: continue
31 + g["q"] = pd.qcut(g[z].rank(method="first"), 5, labels=False) + 1; f = (lambda s: s.median()) if stat == "median" else (lambda s: s.mean()); out[int(yr)] = float(f(g[g.q == 5][oc]) - f(g[g.q == 1][oc]))
32 + return out
33 +def ss(d, z, yrs, oc="fwd_ret", stat="mean"):
34 + sp = cohort_spreads(d, z, yrs, oc, stat); ls = list(sp.values()); mu, t = I.nw_t(ls)
35 + return {"mean": round(mu, 4) if ls else None, "t": round(t, 2) if t == t else None, "t_plain": round(I.plain_t(ls), 2) if len(ls) > 1 else None, "ci95": [round(x, 4) for x in I.ci95(ls)] if len(ls) > 1 else [None, None], "n": len(sp), "by_cohort": {k: round(v, 4) for k, v in sp.items()}}
36 +def ss3(d, z, yrs): return {"raw": ss(d, z, yrs), "w": ss(d, z, yrs, "fwd_w"), "med": ss(d, z, yrs, "fwd_ret", "median")}
37 +def clustered(formula, d, boot_j=None, boot=False):
38 + r = smf.ols(formula, data=d).fit(); X = r.model.exog; y = r.model.endog; gf = pd.factorize(d.ticker)[0]; gy = pd.factorize(d.fiscal_year)[0]
39 + se_f = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gf).cov_params())); se_y = np.sqrt(np.diag(r.get_robustcov_results(cov_type="cluster", groups=gy).cov_params())); se_tw = np.sqrt(np.diag(I.cgm_cov(r, gf, gy)))
40 + out = {"n": int(r.nobs), "coef": {}}
41 + for i, v in enumerate(r.model.exog_names):
42 + if v.startswith("C(") or v == "Intercept": continue
43 + out["coef"][v] = {"coef_pct": round(float(r.params.iloc[i]) * 100, 2), "t_firm": round(float(r.params.iloc[i] / se_f[i]), 2), "t_year": round(float(r.params.iloc[i] / se_y[i]), 2), "t_twoway": round(float(r.params.iloc[i] / se_tw[i]), 2), "ci95_firm_pct": [round(float(r.params.iloc[i] - 1.96 * se_f[i]) * 100, 2), round(float(r.params.iloc[i] + 1.96 * se_f[i]) * 100, 2)]}
44 + if boot and v == boot_j:
45 + for wname in ["rademacher", "webb"]: b = I.wild_cluster_boot(y, X, gy, i, B=9999, weights=wname); out["coef"][v][f"p_wild_{wname}"] = b["p"]; out["coef"][v][f"_tboot_{wname}"] = b["t_boot"]
46 + return out
47 +
48 +# ---------------- A. placebo factors ----------------
49 +FACT = {"AI adoption": "ai_adoption_wmean", "AI infrastructure": "ai_infrastructure_wmean", "AI disruption risk": "ai_disruption_risk_wmean", "ESG commitment": "esg_commitment_wmean", "Climate risk": "climate_risk_wmean", "Cybersecurity": "cybersecurity_wmean", "Automation/labor": "automation_labor_wmean", "Macro risk": "macro_risk_wmean"}
50 +PL = {}; tboot_pl = {}
51 +for nm, col in FACT.items():
52 + z = f"_z_{col}"; df[z] = df.groupby("fiscal_year")[col].transform(zc); d2 = df.dropna(subset=["fwd_w", z, "z_size", "z_mom"]).copy(); d2["zp"] = d2[z] * d2["post"]
53 + c = clustered(f"fwd_w ~ {z} + post + zp + z_size + z_mom", d2, boot_j="zp", boot=True); tboot_pl[nm] = {w: c["coef"]["zp"].pop(f"_tboot_{w}") for w in ["rademacher", "webb"]}
54 + PL[nm] = {"pre": ss(df, z, PRE), "post": ss(df, z, POST), "full": ss(df, z, ALL), "pre_w": ss(df, z, PRE, "fwd_w"), "post_w": ss(df, z, POST, "fwd_w"), "full_w": ss(df, z, ALL, "fwd_w"), "post_med": ss(df, z, POST, "fwd_ret", "median"), "interaction": c["coef"]["zp"], "is_ai": nm.startswith("AI")}
55 +R2["placebo"] = PL; print("A", {k: (v["post_w"]["mean"], v["interaction"]["t_firm"], v["interaction"]["t_twoway"], v["interaction"]["p_wild_webb"]) for k, v in PL.items()}, flush=True)
56 +
57 +# ---------------- B. sector-neutral, double sorts, residual sort, FM +/- Energy ----------------
58 +ta = df.groupby(["fiscal_year", "sector"])["ai_adoption_wmean"].transform(lambda g: zc(g) if g.std(ddof=0) > 0 else g * 0); ti = df.groupby(["fiscal_year", "sector"])["ai_infrastructure_wmean"].transform(lambda g: zc(g) if g.std(ddof=0) > 0 else g * 0); df["z_AI_COMPOS_sn"] = (ta + ti) / 2
59 +R2["sector_neutral"] = {"full": ss(df, "z_AI_COMPOS_sn", ALL), "pre": ss(df, "z_AI_COMPOS_sn", PRE), "post": ss(df, "z_AI_COMPOS_sn", POST), "full_w": ss(df, "z_AI_COMPOS_sn", ALL, "fwd_w"), "pre_w": ss(df, "z_AI_COMPOS_sn", PRE, "fwd_w"), "post_w": ss(df, "z_AI_COMPOS_sn", POST, "fwd_w"), "full_med": ss(df, "z_AI_COMPOS_sn", ALL, "fwd_ret", "median"),
60 + "full_ex_energy_w": ss(df[df.sector != "Energy"], "z_AI_COMPOS_sn", ALL, "fwd_w"), "full_ex_energy": ss(df[df.sector != "Energy"], "z_AI_COMPOS_sn", ALL)}
61 +# residual sort: fwd_w orthogonalised to sector x year FE + controls, then quintiles on the composite
62 +g = df.dropna(subset=["fwd_w", "z_AI_COMPOS", "z_size", "z_mom", "z_vol", "sector"]).copy(); g["sy"] = g.sector + "_" + g.fiscal_year.astype(str)
63 +g["resid"] = smf.ols("fwd_w ~ z_size + z_mom + z_vol + C(sy)", data=g).fit().resid; g["resid_raw"] = smf.ols("fwd_ret ~ z_size + z_mom + z_vol + C(sy)", data=g).fit().resid
64 +R2["residual_sort"] = {"full_w": ss(g, "z_AI_COMPOS", ALL, "resid"), "pre_w": ss(g, "z_AI_COMPOS", PRE, "resid"), "post_w": ss(g, "z_AI_COMPOS", POST, "resid"), "full_raw": ss(g, "z_AI_COMPOS", ALL, "resid_raw"), "n": int(len(g))}
65 +heat = []
66 +for per, msk in [("pre", df.fiscal_year < 2022), ("post", df.fiscal_year >= 2022)]:
67 + for sec, gg in df[msk].groupby("sector"):
68 + gg = gg.dropna(subset=["z_AI_COMPOS", "fwd_w"])
69 + if len(gg) < 30: continue
70 + med = gg.z_AI_COMPOS.median(); heat.append({"period": per, "sector": sec, "hi_minus_lo_pct": round((gg[gg.z_AI_COMPOS > med].fwd_ret.mean() - gg[gg.z_AI_COMPOS <= med].fwd_ret.mean()) * 100, 2), "hi_minus_lo_w_pct": round((gg[gg.z_AI_COMPOS > med].fwd_w.mean() - gg[gg.z_AI_COMPOS <= med].fwd_w.mean()) * 100, 2), "n": len(gg)})
71 +pd.DataFrame(heat).to_csv(f"{T}/t_sector_effect.csv", index=False); R2["sector_effect"] = heat
72 +def double_sort(d, z, by, oc="fwd_ret"):
73 + d = d.dropna(subset=[z, by, oc]).copy(); d["sz"] = d.groupby("fiscal_year")[by].transform(lambda g: pd.qcut(g.rank(method="first"), 3, labels=["Small", "Mid", "Large"]))
74 + d["ai"] = d.groupby(["fiscal_year", "sz"], observed=True)[z].transform(lambda g: pd.qcut(g.rank(method="first"), 3, labels=["Lo", "Mid", "Hi"]))
75 + tab = d.groupby(["sz", "ai"], observed=True)[oc].mean().mul(100).round(2).unstack(); tab["Hi-Lo"] = tab["Hi"] - tab["Lo"]; tab["N"] = d.groupby("sz", observed=True).size(); return tab
76 +for per, msk in [("post", df.fiscal_year >= 2022), ("pre", df.fiscal_year < 2022), ("full", df.fiscal_year >= 2010)]:
77 + R2[f"double_sort_{per}"] = double_sort(df[msk], "z_AI_COMPOS", "log_mktcap").to_dict(); R2[f"double_sort_{per}_w"] = double_sort(df[msk], "z_AI_COMPOS", "log_mktcap", "fwd_w").to_dict(); R2[f"double_sort_{per}_v1cap"] = double_sort(df[msk].assign(lc1=np.log(df[msk].mktcap_v1)), "z_AI_COMPOS", "lc1").to_dict()
78 +print("B sector-neutral full raw/w/med", R2["sector_neutral"]["full"]["mean"], R2["sector_neutral"]["full_w"]["mean"], R2["sector_neutral"]["full_med"]["mean"], "| residual sort full", R2["residual_sort"]["full_w"]["mean"], flush=True)
79 +
80 +# ---------------- C. mega-cap exclusion (ex-top-500 PRE cell computed) ----------------
81 +MAG7 = ["AAPL", "MSFT", "NVDA", "AMZN", "GOOGL", "META", "TSLA", "GOOG", "AVGO"]
82 +R2["mag7"] = {"baseline_post": ss(df, "z_AI_COMPOS", POST), "baseline_pre": ss(df, "z_AI_COMPOS", PRE), "ex_mag7_post": ss(df[~df.ticker.isin(MAG7)], "z_AI_COMPOS", POST), "ex_mag7_pre": ss(df[~df.ticker.isin(MAG7)], "z_AI_COMPOS", PRE),
83 + "ex_top50_post": ss(df[df.cap_rank > 50], "z_AI_COMPOS", POST), "ex_top50_pre": ss(df[df.cap_rank > 50], "z_AI_COMPOS", PRE), "ex_top500_post": ss(df[df.cap_rank > 500], "z_AI_COMPOS", POST), "ex_top500_pre": ss(df[df.cap_rank > 500], "z_AI_COMPOS", PRE),
84 + "ex_mag7_post_w": ss(df[~df.ticker.isin(MAG7)], "z_AI_COMPOS", POST, "fwd_w"), "ex_mag7_pre_w": ss(df[~df.ticker.isin(MAG7)], "z_AI_COMPOS", PRE, "fwd_w"), "ex_top500_post_w": ss(df[df.cap_rank > 500], "z_AI_COMPOS", POST, "fwd_w"), "ex_top500_pre_w": ss(df[df.cap_rank > 500], "z_AI_COMPOS", PRE, "fwd_w")}
85 +
86 +# ---------------- D. alt signals / horizons / benchmark ----------------
87 +for sfx in ["mean", "top5pct"]: df[f"z_AI_{sfx}"] = (df.groupby("fiscal_year")[f"ai_adoption_{sfx}"].transform(zc) + df.groupby("fiscal_year")[f"ai_infrastructure_{sfx}"].transform(zc)) / 2
88 +R2["alt_signal"] = {s: {"pre": ss(df, c, PRE), "post": ss(df, c, POST), "pre_w": ss(df, c, PRE, "fwd_w"), "post_w": ss(df, c, POST, "fwd_w")} for s, c in [("wmean(base)", "z_AI_COMPOS"), ("mean", "z_AI_mean"), ("top5pct", "z_AI_top5pct")]}
89 +R2["alt_horizon"] = {f"{H}m": {"pre": ss(df, "z_AI_COMPOS", PRE, oc=f"fwd_ret_{H}m"), "post": ss(df, "z_AI_COMPOS", range(2022, 2024 if H == 24 else 2025), oc=f"fwd_ret_{H}m"), "pre_w": ss(df, "z_AI_COMPOS", PRE, oc=f"fwd{H}_w"), "post_w": ss(df, "z_AI_COMPOS", range(2022, 2024 if H == 24 else 2025), oc=f"fwd{H}_w")} for H in (6, 24)}
90 +R2["alt_benchmark"] = {"iwv_adj": {"pre": ss(df, "z_AI_COMPOS", PRE, oc="fwd_ret_iwvadj"), "post": ss(df, "z_AI_COMPOS", POST, oc="fwd_ret_iwvadj"), "pre_w": ss(df, "z_AI_COMPOS", PRE, oc="fwd_madj_w"), "post_w": ss(df, "z_AI_COMPOS", POST, oc="fwd_madj_w")}}
91 +
92 +# ---------------- E. year-by-year FM (raw & winsorized) ----------------
93 +fm = []
94 +for yr, g in df.groupby("fiscal_year"):
95 + g = g.dropna(subset=["fwd_ret", "fwd_w", "z_AI_COMPOS", "z_size", "z_mom", "z_vol"])
96 + if len(g) < 40: continue
97 + r = sm.OLS(g["fwd_ret"], sm.add_constant(g[["z_AI_COMPOS", "z_size", "z_mom", "z_vol"]])).fit(cov_type="HC1"); rw = sm.OLS(g["fwd_w"], sm.add_constant(g[["z_AI_COMPOS", "z_size", "z_mom", "z_vol"]])).fit(cov_type="HC1")
98 + fm.append({"year": int(yr), "coef": r.params["z_AI_COMPOS"], "se": r.bse["z_AI_COMPOS"], "t": r.tvalues["z_AI_COMPOS"], "coef_w": rw.params["z_AI_COMPOS"], "se_w": rw.bse["z_AI_COMPOS"], "n": int(r.nobs)})
99 +pd.DataFrame(fm).to_csv(f"{T}/t_fm_byyear.csv", index=False); R2["fm_byyear"] = fm
100 +
101 +# ---------------- G. DiD / event-time with the main effect restored [G] ----------------
102 +exp_c = df[df.fiscal_year.isin([2019, 2020, 2021])].groupby("ticker")[["ai_adoption_wmean", "ai_infrastructure_wmean"]].mean().mean(axis=1)
103 +df["pre_AI_cont"] = df.ticker.map((exp_c - exp_c.mean()) / exp_c.std()); df["treated"] = df.ticker.map((exp_c > exp_c.median()).astype(float))
104 +did = df.dropna(subset=["fwd_w", "pre_AI_cont", "z_size", "z_mom"]).copy(); did["tp"] = did.treated * did.post; did["cp"] = did.pre_AI_cont * did.post; didx = did.set_index(["ticker", "fiscal_year"])
105 +mb = PanelOLS.from_formula("fwd_w ~ tp + z_size + z_mom + EntityEffects + TimeEffects", data=didx, drop_absorbed=True).fit(cov_type="clustered", cluster_entity=True, cluster_time=True)
106 +mc = PanelOLS.from_formula("fwd_w ~ cp + z_size + z_mom + EntityEffects + TimeEffects", data=didx, drop_absorbed=True).fit(cov_type="clustered", cluster_entity=True, cluster_time=True)
107 +R2["did"] = {"binary": {"coef": round(float(mb.params["tp"]), 4), "t": round(float(mb.tstats["tp"]), 2), "n": int(mb.nobs)}, "continuous": {"coef": round(float(mc.params["cp"]), 4), "t": round(float(mc.tstats["cp"]), 2), "n": int(mc.nobs)}}
108 +def event_time(d, expo, years, ref=2021, label=""):
109 + ev = d.dropna(subset=["fwd_w", expo, "z_size", "z_mom"]).copy(); ev = ev[ev.fiscal_year.between(years[0], years[-1])]; ev["yr"] = ev.fiscal_year.astype(int).astype(str)
110 + # (i) author's coding (no main effect): absolute slope per year
111 + me0 = smf.ols(f"fwd_w ~ C(yr, Treatment('{ref}')):{expo} + C(yr) + z_size + z_mom", data=ev).fit(cov_type="cluster", cov_kwds={"groups": ev["ticker"]})
112 + # (ii) main effect restored: coefficients are DIFFERENCES from the reference year
113 + me1 = smf.ols(f"fwd_w ~ {expo} + C(yr, Treatment('{ref}')):{expo} + C(yr) + z_size + z_mom", data=ev).fit(cov_type="cluster", cov_kwds={"groups": ev["ticker"]})
114 + rows = []
115 + for yr in years:
116 + k0 = next((k for k in me0.params.index if k.startswith(f"C(yr, Treatment('{ref}'))[") and k.endswith(f"{yr}]:{expo}")), None); k1 = next((k for k in me1.params.index if k.startswith(f"C(yr, Treatment('{ref}'))[T.") and k.endswith(f"{yr}]:{expo}")), None)
117 + rows.append({"year": int(yr), "abs_slope": float(me0.params[k0]) if k0 else None, "abs_se": float(me0.bse[k0]) if k0 else None, "diff_from_ref": float(me1.params[k1]) if k1 else 0.0, "diff_se": float(me1.bse[k1]) if k1 else 0.0, "n": int(me1.nobs)})
118 + base = float(me1.params[expo]); return {"ref": ref, "base_slope_ref_year": round(base, 4), "base_t": round(float(me1.tvalues[expo]), 2), "rows": rows, "n": int(me1.nobs), "label": label}
119 +R2["did_eventtime"] = event_time(df, "pre_AI_cont", list(range(2014, 2025)), label="exposure = 2019-21 mean composite (author), outcomes 2014-2024: pre-2019 rows are look-back, not pre-trends")
120 +pd.DataFrame(R2["did_eventtime"]["rows"]).to_csv(f"{T}/t_did_eventtime.csv", index=False)
121 +# trailing exposure: mean raw composite over Y-3..Y-1, standardised within year (strictly pre-determined for every outcome year)
122 +comp_raw = df.groupby(["ticker", "fiscal_year"])[["ai_adoption_wmean", "ai_infrastructure_wmean"]].mean().mean(axis=1).rename("c").reset_index()
123 +tr = []
124 +for lag in (1, 2, 3): x = comp_raw.copy(); x["fiscal_year"] = x["fiscal_year"] + lag; tr.append(x.rename(columns={"c": f"c_l{lag}"}))
125 +tmp = comp_raw[["ticker", "fiscal_year"]].copy()
126 +for x in tr: tmp = tmp.merge(x, on=["ticker", "fiscal_year"], how="left")
127 +tmp["trail"] = tmp[["c_l1", "c_l2", "c_l3"]].mean(axis=1); df = df.merge(tmp[["ticker", "fiscal_year", "trail"]], on=["ticker", "fiscal_year"], how="left"); df["trail_z"] = df.groupby("fiscal_year")["trail"].transform(zc)
128 +R2["did_eventtime_trailing"] = event_time(df, "trail_z", list(range(2014, 2025)), label="exposure = trailing 3-year mean composite ending Y-1 (pre-determined)")
129 +R2["did_eventtime_post2019"] = event_time(df, "pre_AI_cont", list(range(2019, 2025)), label="author exposure, cohorts >= 2019 only")
130 +pd.DataFrame(R2["did_eventtime_trailing"]["rows"]).to_csv(f"{T}/t_did_eventtime_trailing.csv", index=False)
131 +print("G", R2["did"], [(r["year"], round(r["diff_from_ref"] * 100, 1)) for r in R2["did_eventtime"]["rows"] if r["year"] >= 2021], flush=True)
132 +
133 +# ---------------- H. forward vol ----------------
134 +R2["fwd_vol"] = {"pre": ss(df, "z_AI_COMPOS", PRE, oc="fwd_vol"), "post": ss(df, "z_AI_COMPOS", POST, oc="fwd_vol"), "full": ss(df, "z_AI_COMPOS", ALL, oc="fwd_vol")}
135 +
136 +# ---------------- I. cohorts, LOO, pooled (pooled blocks from r3k_15) ----------------
137 +def cohort_block(oc="fwd_ret", stat="mean"):
138 + cs = cohort_spreads(df, "z_AI_COMPOS", ALL, oc, stat); post = {k: v for k, v in cs.items() if k >= 2022}; pre = [v for k, v in cs.items() if k < 2022]
139 + return {"cohort_spreads_pct": {k: round(v * 100, 1) for k, v in cs.items()}, "post_mean_pct": round(np.mean(list(post.values())) * 100, 1), "post_loo_pct": {k: round(np.mean([v for kk, v in post.items() if kk != k]) * 100, 1) for k in post}, "pre_mean_pct": round(np.mean(pre) * 100, 1),
140 + "pre_ci95_pct": [round(x * 100, 1) for x in I.ci95(pre)], "post_ci95_pct": [round(x * 100, 1) for x in I.ci95(list(post.values()))], "pre_cohorts_positive_share": round(float(np.mean([v > 0 for v in pre])), 2), "n_pre_above_2022": int(sum(v > cs[2022] for k, v in cs.items() if k < 2022)), "pre_min_pct": round(min(pre) * 100, 1)}
141 +R2["cohorts"] = cohort_block(); R2["cohorts_w"] = cohort_block("fwd_w"); R2["cohorts_med"] = cohort_block("fwd_ret", "median")
142 +for k in ["post", "post_w", "pre", "pre_w", "post_excl_2024", "post_excl_2024_w", "diff", "diff_w"]: R2["cohorts"]["pooled_" + k] = RES["pooled"][k]
143 +print("I raw", R2["cohorts"]["cohort_spreads_pct"], "| wins", R2["cohorts_w"]["cohort_spreads_pct"], flush=True)
144 +
145 +# ---------------- J. alternative splits with pooled firm-clustered tests [M] ----------------
146 +def pooled_split(d, pre_yrs, post_yrs, dv="fwd_w"):
147 + d = d[d.fiscal_year.isin(list(pre_yrs) + list(post_yrs))].dropna(subset=["z_AI_COMPOS", dv]).copy(); d["q"] = d.groupby("fiscal_year")["z_AI_COMPOS"].transform(lambda s: pd.qcut(s.rank(method="first"), 5, labels=False) + 1); d = d[d.q.isin([1, 5])].copy()
148 + d["hi"] = (d.q == 5).astype(float); d["p2"] = d.fiscal_year.isin(list(post_yrs)).astype(float); d["hp"] = d.hi * d.p2; d["yr"] = d.fiscal_year.astype(str); c = clustered(f"{dv} ~ hi + hp + C(yr)", d); return c["coef"]["hp"] | {"n": c["n"]}
149 +R2["breaks"] = {}
150 +for lab, pre_y, post_y in [("split_2022", range(2010, 2022), range(2022, 2025)), ("split_2021", range(2010, 2021), range(2021, 2025))]:
151 + R2["breaks"][lab] = {"pre": ss(df, "z_AI_COMPOS", pre_y), "post": ss(df, "z_AI_COMPOS", post_y), "pre_w": ss(df, "z_AI_COMPOS", pre_y, "fwd_w"), "post_w": ss(df, "z_AI_COMPOS", post_y, "fwd_w"), "pre_med": ss(df, "z_AI_COMPOS", pre_y, "fwd_ret", "median"), "post_med": ss(df, "z_AI_COMPOS", post_y, "fwd_ret", "median"),
152 + "pooled_diff_w": pooled_split(df, pre_y, post_y), "pooled_diff_raw": pooled_split(df, pre_y, post_y, "fwd_ret")}
153 +R2["breaks"]["post_2022_2023"] = ss(df, "z_AI_COMPOS", range(2022, 2024)); R2["breaks"]["post_2024"] = ss(df, "z_AI_COMPOS", range(2024, 2025)); R2["breaks"]["cohort_2021_w"] = ss(df, "z_AI_COMPOS", range(2021, 2022), "fwd_w")
154 +print("J", {k: (v["pre_w"]["mean"], v["post_w"]["mean"], v["pooled_diff_w"]) for k, v in R2["breaks"].items() if k.startswith("split")}, flush=True)
155 +
156 +# ---------------- K. entrants / seasoned ----------------
157 +seasoned = df[df.first_price <= "2009-07-31"]; l36 = df[df.months_listed >= 36]
158 +df["q"] = df.groupby("fiscal_year")["z_AI_COMPOS"].transform(lambda s: pd.qcut(s.rank(method="first"), 5, labels=False) + 1); df["entrant"] = (df.first_price > "2009-07-31").astype(int)
159 +R2["entrants"] = {"n_firms": int(df.ticker.nunique()), "n_seasoned": int(seasoned.ticker.nunique()), "seasoned": {"pre": ss(seasoned, "z_AI_COMPOS", PRE), "post": ss(seasoned, "z_AI_COMPOS", POST), "pre_w": ss(seasoned, "z_AI_COMPOS", PRE, "fwd_w"), "post_w": ss(seasoned, "z_AI_COMPOS", POST, "fwd_w"), "full_w": ss(seasoned, "z_AI_COMPOS", ALL, "fwd_w")},
160 + "listed36": {"pre": ss(l36, "z_AI_COMPOS", PRE), "post": ss(l36, "z_AI_COMPOS", POST), "pre_w": ss(l36, "z_AI_COMPOS", PRE, "fwd_w"), "post_w": ss(l36, "z_AI_COMPOS", POST, "fwd_w")}, "entrant_share_by_q": df.groupby("q").entrant.mean().round(3).to_dict(),
161 + "entrant_minus_seasoned_fwd_pct": round((df[df.entrant == 1].fwd_ret.mean() - df[df.entrant == 0].fwd_ret.mean()) * 100, 1), "entrant_minus_seasoned_fwd_w_pct": round((df[df.entrant == 1].fwd_w.mean() - df[df.entrant == 0].fwd_w.mean()) * 100, 1)}
162 +
163 +# ---------------- L. sector controls: FM +/- sector dummies, +/- Energy, raw & winsorized; pooled with two-way + wild ----------------
164 +FB = I.fixed_b_cv(15, 3)["cv"]
165 +def fm_series(d, oc, with_sector):
166 + out = {}
167 + for yr, g in d.groupby("fiscal_year"):
168 + g = g.dropna(subset=[oc, "z_AI_COMPOS", "z_size", "z_mom", "z_vol", "sector"])
169 + if len(g) < 40: continue
170 + out[int(yr)] = float(smf.ols(f"{oc} ~ z_AI_COMPOS + z_size + z_mom + z_vol" + (" + C(sector)" if with_sector else ""), data=g).fit().params["z_AI_COMPOS"])
171 + return out
172 +L = {}; fm_boot = {}
173 +for samp, dd in [("all", df), ("ex_energy", df[df.sector != "Energy"])]:
174 + for oc, tag in [("fwd_ret", "raw"), ("fwd_w", "w")]:
175 + for ws, stag in [(False, "nosector"), (True, "sector")]:
176 + s = fm_series(dd, oc, ws); fm_boot[f"{samp}_{tag}_{stag}"] = s
177 + for lab, yrs in [("pre", PRE), ("post", POST), ("full", ALL)]:
178 + v = [s[y] for y in s if y in yrs]; mu, t = I.nw_t(v); L[f"fm_{samp}_{tag}_{stag}_{lab}"] = {"slope_pct": round(mu * 100, 2), "t_nw": round(t, 2) if t == t else None, "t_plain": round(I.plain_t(v), 2) if len(v) > 1 else None, "fixed_b_cv05": FB if lab == "full" else None, "ci95_pct": [round(x * 100, 2) for x in I.ci95(v)] if len(v) > 1 else [None, None], "n_years": len(v)}
179 + for lab, yrs in [("pre", PRE), ("post", POST), ("full", ALL)]:
180 + g = dd[dd.fiscal_year.isin(list(yrs))].dropna(subset=["fwd_w", "z_AI_COMPOS", "z_size", "z_mom", "z_vol", "sector"]).copy(); g["sy"] = g.sector + "_" + g.fiscal_year.astype(str); g["yr"] = g.fiscal_year.astype(str)
181 + c1 = clustered("fwd_w ~ z_AI_COMPOS + z_size + z_mom + z_vol + C(sy)", g, boot_j="z_AI_COMPOS", boot=(samp == "all")); c0 = clustered("fwd_w ~ z_AI_COMPOS + z_size + z_mom + z_vol + C(yr)", g, boot_j="z_AI_COMPOS", boot=(samp == "all"))
182 + for c in (c0, c1):
183 + for w in ["rademacher", "webb"]: c["coef"]["z_AI_COMPOS"].pop(f"_tboot_{w}", None)
184 + L[f"pooled_{samp}_{lab}"] = {"year_fe": c0["coef"]["z_AI_COMPOS"], "sector_year_fe": c1["coef"]["z_AI_COMPOS"], "n": c1["n"]}
185 +R2["sector_controls"] = L; print("L full all: FM raw sector", L["fm_all_raw_sector_full"], "| FM w sector", L["fm_all_w_sector_full"], "| pooled sy", L["pooled_all_full"]["sector_year_fe"], flush=True)
186 +
187 +# ---------------- M. TECH-factor spanning ----------------
188 +sect = df.drop_duplicates("ticker").set_index("ticker")["sector"]; mret["sector"] = mret.ticker.map(sect)
189 +tech = (mret[mret.sector == "Information Technology"].groupby("date").ret.mean() - mret.groupby("date").ret.mean()).rename("TECH").replace([np.inf, -np.inf], np.nan)
190 +M = {}
191 +for k in ["AI_COMPOS", "AI_ADOPT", "AI_INFRA", "AI_RISK", "AI_PC1", "AI_formw", "SPID_TECH"]:
192 + ts = pd.read_parquet(f"{V}/ls_ts_v2_{k}.parquet"); ts["date"] = pd.to_datetime(ts["date"]); ts = ts.set_index("date").join(ff, how="left").join(tech, how="left").dropna(subset=["ls"]); M[k] = {}
193 + for lab, msk in [("full", ts.index >= ts.index.min()), ("pre", ts.index < POST_START), ("post", ts.index >= POST_START)]:
194 + s = ts[msk]; mu, t = I.nw_t(s["ls"].values); blk = {"n_months": int(len(s)), "mean_ann_pct": round(mu * 1200, 2), "sharpe": round(mu / s["ls"].std() * np.sqrt(12), 2), "hac_t": round(t, 2), "plain_t": round(I.plain_t(s["ls"].values), 2)}
195 + for nm, cols in [("capm", ["MktRF"]), ("ff6", ["MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]), ("ff6_tech", ["MktRF", "SMB", "HML", "RMW", "CMA", "MOM", "TECH"])]:
196 + s2 = s.dropna(subset=cols)
197 + if len(s2) < 12: continue
198 + r = sm.OLS(s2["ls"], sm.add_constant(s2[cols])).fit(cov_type="HAC", cov_kwds={"maxlags": 3}); blk[f"{nm}_alpha_pct"] = round(r.params["const"] * 1200, 2); blk[f"{nm}_t"] = round(r.tvalues["const"], 2)
199 + if "TECH" in cols: blk["tech_beta"] = round(r.params["TECH"], 3); blk["tech_beta_t"] = round(r.tvalues["TECH"], 2)
200 + M[k][lab] = blk
201 +R2["monthly_ls"] = M; print("M", M["AI_COMPOS"]["post"], flush=True)
202 +
203 +# ---------------- N. filing-type composition, Wald equality, 2022-cohort post-launch share and re-run [H, M] ----------------
204 +fs = pd.read_parquet(f"{R}/filing_scores"); fl = pd.read_parquet(f"{R}/filings")[["accession_number", "form_type", "filing_date"]].rename(columns={"filing_date": "fdate"})
205 +fs = fs.drop(columns=[c for c in ["form_type"] if c in fs.columns]).merge(fl, on="accession_number", how="left"); fs["fiscal_year"] = pd.to_datetime(fs["filing_date"]).dt.year
206 +N = {"filings_scored_by_form": fs.form_type.value_counts().to_dict()}
207 +for tag, lab in [("10k", "10-K"), ("10q", "10-Q"), ("8k", "8-K"), ("formw", "equal-form")]:
208 + z = f"z_AI_{tag}"; N[tag] = {"share_firm_years": round(float(df[z].notna().mean()), 3), "corr_with_baseline": round(float(df[["z_AI_COMPOS", z]].corr().iloc[0, 1]), 3), "pre": ss(df, z, PRE), "post": ss(df, z, POST), "full": ss(df, z, ALL), "pre_w": ss(df, z, PRE, "fwd_w"), "post_w": ss(df, z, POST, "fwd_w"), "full_w": ss(df, z, ALL, "fwd_w"),
209 + "fm_full": {"slope_pct": round(RES["fama_macbeth"]["tenk_w" if tag == "10k" else "formw_w" if tag == "formw" else "compos_w"][z if tag in ("10k", "formw") else "z_AI_COMPOS"]["mean_coef"] * 100, 2)} if tag in ("10k", "formw") else None}
210 +# FM per form (raw and winsorized) + joint horse race + HAC Wald test of equality of the three slope series
211 +def fm_form(oc):
212 + coefs = {k: [] for k in ["z_AI_10k", "z_AI_10q", "z_AI_8k"]}; single = {k: [] for k in coefs}
213 + for yr, g in df.groupby("fiscal_year"):
214 + g = g.dropna(subset=[oc, "z_AI_10k", "z_AI_10q", "z_AI_8k", "z_size", "z_mom", "z_vol"])
215 + if len(g) < 40: continue
216 + r = sm.OLS(g[oc], sm.add_constant(g[list(coefs) + ["z_size", "z_mom", "z_vol"]])).fit()
217 + for k in coefs: coefs[k].append(r.params[k]); single[k].append(sm.OLS(g[oc], sm.add_constant(g[[k, "z_size", "z_mom", "z_vol"]])).fit().params[k])
218 + Cm = np.array([coefs[k] for k in coefs]).T; Tn = Cm.shape[0]; mu = Cm.mean(0)
219 + # HAC (Bartlett, 3 lags) long-run covariance of the 3 slope series -> Wald test of H0: b_10k = b_10q = b_8k
220 + Xc = Cm - mu; S = Xc.T @ Xc / Tn
221 + for l in range(1, 4): G_ = Xc[l:].T @ Xc[:-l] / Tn; S += (1 - l / 4) * (G_ + G_.T)
222 + Rm = np.array([[1, -1, 0], [1, 0, -1]]); W = float((Rm @ mu).T @ np.linalg.inv(Rm @ (S / Tn) @ Rm.T) @ (Rm @ mu)); from scipy import stats as st
223 + return {"horse_race": {k: {"slope_pct": round(I.nw_t(v)[0] * 100, 2), "t": round(I.nw_t(v)[1], 2), "t_plain": round(I.plain_t(v), 2)} for k, v in coefs.items()}, "single": {k: {"slope_pct": round(I.nw_t(v)[0] * 100, 2), "t": round(I.nw_t(v)[1], 2), "t_plain": round(I.plain_t(v), 2)} for k, v in single.items()},
224 + "wald_equality": {"stat": round(W, 3), "df": 2, "p": round(float(1 - st.chi2.cdf(W, 2)), 4), "n_years": Tn}}
225 +N["fm_forms_raw"] = fm_form("fwd_ret"); N["fm_forms_w"] = fm_form("fwd_w"); N["horse_race_fm"] = N["fm_forms_raw"]["horse_race"]
226 +N["corr_matrix"] = df[["z_AI_COMPOS", "z_AI_10k", "z_AI_10q", "z_AI_8k", "z_AI_formw"]].corr().round(3).to_dict()
227 +# post-launch share of the 2022 cohort's text and the pre-launch-only re-run (eq. (3) = unweighted mean of filing wmeans, verified in r3k_16 log)
228 +f22 = fs[fs.fiscal_year == 2022].copy(); f22["post_launch"] = pd.to_datetime(f22.filing_date) > pd.Timestamp("2022-11-30")
229 +N["cohort2022"] = {"share_filings_post_launch": round(float(f22.post_launch.mean()), 4), "share_chunks_post_launch": round(float(f22.loc[f22.post_launch, "chunk_count"].sum() / f22.chunk_count.sum()), 4), "share_filings_post_launch_by_form": f22.groupby("form_type").post_launch.mean().round(4).to_dict()}
230 +chk = fs.groupby(["ticker", "fiscal_year"])[["ai_adoption_wmean", "ai_infrastructure_wmean"]].mean().reset_index().merge(df[["ticker", "fiscal_year", "ai_adoption_wmean"]].rename(columns={"ai_adoption_wmean": "panel"}), on=["ticker", "fiscal_year"])
231 +N["cohort2022"]["aggregation_reproduced_share_1e6"] = round(float(((chk.ai_adoption_wmean - chk.panel).abs() < 1e-6).mean()), 5)
232 +pre_only = f22[~f22.post_launch].groupby("ticker")[["ai_adoption_wmean", "ai_infrastructure_wmean"]].mean().reset_index(); d22 = df[df.fiscal_year == 2022].drop(columns=[c for c in ["ad22", "in22"] if c in df.columns]).merge(pre_only.rename(columns={"ai_adoption_wmean": "ad22", "ai_infrastructure_wmean": "in22"}), on="ticker", how="inner")
233 +d22["z22"] = (zc(d22.ad22) + zc(d22.in22)) / 2; cs22 = cohort_spreads(d22, "z22", [2022]); cs22w = cohort_spreads(d22, "z22", [2022], "fwd_w"); cs22b = cohort_spreads(d22, "z_AI_COMPOS", [2022]); cs22bw = cohort_spreads(d22, "z_AI_COMPOS", [2022], "fwd_w")
234 +N["cohort2022"]["spread_prelaunch_text_only_pct"] = round(cs22.get(2022, np.nan) * 100, 2); N["cohort2022"]["spread_prelaunch_text_only_w_pct"] = round(cs22w.get(2022, np.nan) * 100, 2); N["cohort2022"]["spread_full_text_same_firms_pct"] = round(cs22b.get(2022, np.nan) * 100, 2); N["cohort2022"]["spread_full_text_same_firms_w_pct"] = round(cs22bw.get(2022, np.nan) * 100, 2); N["cohort2022"]["n_firms_with_prelaunch_filings"] = int(len(d22)); N["cohort2022"]["rank_corr_prelaunch_vs_full"] = round(float(d22[["z22", "z_AI_COMPOS"]].corr(method="spearman").iloc[0, 1]), 3)
235 +R2["filing_type"] = N; print("N", {k: (v["pre_w"]["mean"], v["post_w"]["mean"], v["corr_with_baseline"]) for k, v in N.items() if isinstance(v, dict) and "pre_w" in v}, "| Wald", N["fm_forms_w"]["wald_equality"], "| 2022", N["cohort2022"], flush=True)
236 +
237 +# ---------------- O. placebo break dates, permutation, sup-Wald [M] ----------------
238 +O = {}
239 +for tag, oc, st in [("raw", "fwd_ret", "mean"), ("w", "fwd_w", "mean"), ("med", "fwd_ret", "median")]:
240 + cs_all = cohort_spreads(df, "z_AI_COMPOS", ALL, oc, st); diffs = {}
241 + for b in range(2013, 2023):
242 + pre_ = [v for k, v in cs_all.items() if k < b]; post_ = [v for k, v in cs_all.items() if k >= b]; diffs[b] = round((np.mean(post_) - np.mean(pre_)) * 100, 1)
243 + O[tag] = {"diff_post_minus_pre_pct_by_break_year": diffs, "rank_of_2022": int(sorted(diffs.values(), reverse=True).index(diffs[2022]) + 1), "n_candidates": len(diffs), "permutation": I.permutation_three_cohorts(cs_all, [2022, 2023, 2024]), "sup_wald": I.sup_wald([cs_all[y] for y in sorted(cs_all)])}
244 +R2["placebo_breaks"] = O["raw"] | {"winsorized": O["w"], "median": O["med"]}; print("O raw perm p", O["raw"]["permutation"]["p_one_sided"], "rank", O["raw"]["permutation"]["rank_among_subsets"], "| wins", O["w"]["permutation"]["p_one_sided"], "| sup-Wald raw", O["raw"]["sup_wald"]["sup_wald"], O["raw"]["sup_wald"]["p_simulated"], flush=True)
245 +
246 +# ---------------- P. SPID interactions on the rebuilt panel (explicit terms) [C, I] ----------------
247 +P = {}
248 +for k in ["z_SPID_TECH", "spid_ai_any", "spid_tech_any"]:
249 + d2 = df.dropna(subset=["fwd_w", k, "z_AI_COMPOS", "z_size", "z_mom"]).copy(); d2["kp"] = d2[k] * d2.post; d2["cp"] = d2.z_AI_COMPOS * d2.post
250 + c = clustered(f"fwd_w ~ {k} + post + kp + z_AI_COMPOS + cp + z_size + z_mom", d2); X = smf.ols(f"fwd_w ~ {k} + post + kp + z_AI_COMPOS + cp + z_size + z_mom", data=d2).fit().model.exog
251 + P[k] = {"terms": {{"kp": f"{k}:post", "cp": "z_AI_COMPOS:post"}.get(v, v): vals for v, vals in c["coef"].items()}, "n": c["n"], "rank": int(np.linalg.matrix_rank(X)), "cond": round(float(np.linalg.cond(X)), 2)}
252 +P["share_firm_years_any_ai_project"] = RES["spid"]["share_firm_years_any_ai"]; P["share_post_any_ai_project"] = RES["spid"]["share_post_any_ai"]; P["ai_any_by_quintile_post"] = df[df.post == 1].groupby("q").spid_ai_any.mean().round(3).to_dict()
253 +P["binary_contrasts"] = {k: RES["spid"][k] for k in ["spid_ai_any", "spid_tech_any"]}; P["tercile_nonzero_w"] = RES["spid"]["tercile_T3_minus_T1_nonzero_w"]; R2["spid"] = P; print("P", {k: P[k]["terms"][f"{k}:post"] for k in ["z_SPID_TECH", "spid_ai_any"]}, flush=True)
254 +
255 +# ---------------- Q. reversal (AI-washing) and cyber-neutral sort [K] ----------------
256 +rv = df[(df.fiscal_year <= 2023)].dropna(subset=["fwd_ret", "fwd_ret_24m", "z_AI_COMPOS"]).copy(); rv["r13_24"] = (1 + rv.fwd_ret_24m) / (1 + rv.fwd_ret) - 1; rv["r13_24_w"] = rv.groupby("fiscal_year")["r13_24"].transform(wins_within)
257 +rv["hi_lang"] = rv.groupby("fiscal_year")["z_AI_COMPOS"].transform(lambda s: s >= s.quantile(0.8)); hl = rv[rv.hi_lang].copy(); hl["proj"] = (hl.spid_mentions_tech > 0).astype(int)
258 +def contrast_by_cohort(d, flag, oc, min_n=15):
259 + out = {}
260 + for yr, g in d.groupby("fiscal_year"):
261 + a = g[g[flag] == 1][oc].dropna(); b = g[g[flag] == 0][oc].dropna()
262 + if len(a) >= min_n and len(b) >= min_n: out[int(yr)] = float(a.mean() - b.mean())
263 + ls = list(out.values()); mu, t = I.nw_t(ls); return {"mean_pct": round(mu * 100, 2) if ls else None, "t": round(t, 2) if t == t else None, "t_plain": round(I.plain_t(ls), 2) if len(ls) > 1 else None, "ci95_pct": [round(x * 100, 2) for x in I.ci95(ls)] if len(ls) > 1 else [None, None], "n_cohorts": len(ls), "by_cohort_pct": {k: round(v * 100, 1) for k, v in out.items()}}
264 +Q = {"reversal": {"high_language_with_vs_without_projects": {"month_1_12_w": contrast_by_cohort(hl, "proj", "fwd_w"), "month_13_24_w": contrast_by_cohort(hl, "proj", "r13_24_w"), "month_1_24_raw": contrast_by_cohort(hl, "proj", "fwd_ret_24m"), "n_high_language_with_projects": int(hl.proj.sum()), "n_high_language": int(len(hl))},
265 + "Q5_minus_Q1_month_13_24": {"raw": ss(rv, "z_AI_COMPOS", range(2010, 2024), "r13_24"), "w": ss(rv, "z_AI_COMPOS", range(2010, 2024), "r13_24_w"), "pre_w": ss(rv, "z_AI_COMPOS", range(2010, 2022), "r13_24_w"), "post_w": ss(rv, "z_AI_COMPOS", range(2022, 2024), "r13_24_w")},
266 + "Q5_minus_Q1_month_1_24": {"raw": ss(rv, "z_AI_COMPOS", range(2010, 2024), "fwd_ret_24m"), "w": ss(rv, "z_AI_COMPOS", range(2010, 2024), "fwd24_w")}}}
267 +df["cyber_terc"] = df.groupby("fiscal_year")["z_CYBER"].transform(lambda s: pd.qcut(s.rank(method="first"), 3, labels=False)); df["z_AI_COMPOS_cyberneu"] = df.groupby(["fiscal_year", "cyber_terc"])["z_AI_COMPOS"].transform(zc)
268 +Q["cyber"] = {"corr_cyber_composite": round(float(df[["z_CYBER", "z_AI_COMPOS"]].corr().iloc[0, 1]), 3), "corr_cyber_adopt": round(float(df[["z_CYBER", "z_AI_ADOPT"]].corr().iloc[0, 1]), 3), "cyberneutral_sort": {"full": ss(df, "z_AI_COMPOS_cyberneu", ALL), "full_w": ss(df, "z_AI_COMPOS_cyberneu", ALL, "fwd_w"), "pre_w": ss(df, "z_AI_COMPOS_cyberneu", PRE, "fwd_w"), "post_w": ss(df, "z_AI_COMPOS_cyberneu", POST, "fwd_w")},
269 + "cyber_sort": {"full_w": ss(df, "z_CYBER", ALL, "fwd_w"), "pre_w": ss(df, "z_CYBER", PRE, "fwd_w"), "post_w": ss(df, "z_CYBER", POST, "fwd_w")}, "fm_compos_with_cyber_w": RES["fama_macbeth"]["compos_cyber_w"]}
270 +R2["cluster_k"] = Q; print("Q reversal 13-24 hi-lang proj vs no", Q["reversal"]["high_language_with_vs_without_projects"]["month_13_24_w"], "| cyber-neutral full w", Q["cyber"]["cyberneutral_sort"]["full_w"]["mean"], "corr", Q["cyber"]["corr_cyber_composite"], flush=True)
271 +
272 +# ---------------- R. multiple-testing families [N] ----------------
273 +rng = np.random.default_rng(I.RNG_SEED); B = 9999
274 +def family_cohort(series_dict):
275 + """series_dict: name -> list of cohort values. Shared cohort-resampling bootstrap (same draws for every member), recentred t's."""
276 + names = list(series_dict); n = min(len(v) for v in series_dict.values()); M = np.array([np.asarray(series_dict[k][:n], float) for k in names]).T # cohorts x members (aligned on the first n cohorts)
277 + t_obs = M.mean(0) / (M.std(0, ddof=1) / np.sqrt(n)); idx = rng.integers(0, n, size=(B, n)); Mb = M[idx] # B x n x m
278 + mb = Mb.mean(1); sb = Mb.std(1, ddof=1) / np.sqrt(n); t_boot = (mb - M.mean(0)) / sb
279 + p_plain = 2 * (1 - __import__("scipy").stats.t.cdf(np.abs(t_obs), n - 1)); return {"members": names, "t": [round(float(x), 2) for x in t_obs], "p": [round(float(x), 4) for x in p_plain], "p_holm": [round(float(x), 4) for x in I.holm(p_plain)], "p_bh": [round(float(x), 4) for x in I.bh(p_plain)], "p_romano_wolf": [round(float(x), 4) for x in I.romano_wolf(t_obs, t_boot)], "family_size": len(names), "n_cohorts": int(n), "hlz_t_gt_3": [bool(abs(x) > 3) for x in t_obs]}
280 +MHT = {}
281 +# (a) Table 5 / 9 family: 5 signals x {EW, VW} x {full, pre, post}, cohort-level spreads on winsorized returns (full = all cohorts; pre/post subsets need aligned lengths -> separate sub-families)
282 +for per, yrs in [("full", ALL), ("pre", PRE), ("post", POST)]:
283 + fam = {}
284 + for k in ["AI_ADOPT", "AI_INFRA", "AI_RISK", "AI_COMPOS", "SPID_TECH"]:
285 + for wt in ["EW", "VW"]:
286 + s = RES["sorts"][f"{k}_{wt}_w"]["by_cohort"]; fam[f"{k}_{wt}"] = [s[str(y)] if str(y) in s else s[y] for y in sorted(int(x) for x in s) if int(y) in yrs] if False else [v for y, v in sorted((int(y), v) for y, v in s.items()) if y in yrs]
287 + MHT[f"sorts_{per}_w"] = family_cohort(fam)
288 +# (b) Table 14: 8 placebo interactions on shared wild-bootstrap draws (identical seed and G -> identical weight matrix)
289 +names = list(PL); t_obs = np.array([PL[n]["interaction"]["t_year"] for n in names]); tb = np.column_stack([tboot_pl[n]["webb"] for n in names])
290 +p_year = 2 * (1 - __import__("scipy").stats.t.cdf(np.abs(t_obs), 14)); MHT["placebo_interactions_yearcl"] = {"members": names, "t_year": [round(float(x), 2) for x in t_obs], "t_firm": [PL[n]["interaction"]["t_firm"] for n in names], "p_wild_webb": [PL[n]["interaction"]["p_wild_webb"] for n in names], "p_holm_on_wild": [round(float(x), 4) for x in I.holm([PL[n]["interaction"]["p_wild_webb"] for n in names])], "p_bh_on_wild": [round(float(x), 4) for x in I.bh([PL[n]["interaction"]["p_wild_webb"] for n in names])], "p_romano_wolf": [round(float(x), 4) for x in I.romano_wolf(t_obs, tb)], "family_size": len(names)}
291 +# (c) ten break dates: post - pre mean difference, cohort bootstrap shared
292 +cs_w = cohort_spreads(df, "z_AI_COMPOS", ALL, "fwd_w"); yrs_ = sorted(cs_w); v = np.array([cs_w[y] for y in yrs_]); stats_obs = []; names = []
293 +for b in range(2013, 2023):
294 + m = np.array([y >= b for y in yrs_]); names.append(str(b)); stats_obs.append(v[m].mean() - v[~m].mean())
295 +idx = rng.integers(0, len(v), size=(B, len(v))); boot = []
296 +for bi in range(B):
297 + vb = v[idx[bi]]; yb = np.array(yrs_)[idx[bi]]; boot.append([(vb[yb >= b].mean() - vb[yb < b].mean()) if (yb >= b).any() and (yb < b).any() else 0.0 for b in range(2013, 2023)])
298 +boot = np.array(boot) - np.array(stats_obs); sd = boot.std(0, ddof=1); t_obs = np.array(stats_obs) / sd; t_boot = boot / sd
299 +MHT["break_dates_w"] = {"members": names, "diff_pct": [round(x * 100, 1) for x in stats_obs], "t_boot_se": [round(float(x), 2) for x in t_obs], "p_boot": [round(float((np.abs(t_boot[:, i]) >= abs(t_obs[i])).mean()), 4) for i in range(len(names))], "p_romano_wolf": [round(float(x), 4) for x in I.romano_wolf(t_obs, t_boot)], "family_size": len(names)}
300 +# (e) three forms: FM slope series (winsorized), cohort bootstrap
301 +MHT["forms_fm_w"] = family_cohort({k: v for k, v in zip(["10k", "10q", "8k"], [[], [], []])} if False else {k: [] for k in []}) if False else None
302 +fmf = {};
303 +for yr, g in df.groupby("fiscal_year"):
304 + g = g.dropna(subset=["fwd_w", "z_AI_10k", "z_AI_10q", "z_AI_8k", "z_size", "z_mom", "z_vol"])
305 + if len(g) < 40: continue
306 + for k in ["z_AI_10k", "z_AI_10q", "z_AI_8k"]: fmf.setdefault(k, []).append(sm.OLS(g["fwd_w"], sm.add_constant(g[[k, "z_size", "z_mom", "z_vol"]])).fit().params[k])
307 +MHT["forms_fm_w"] = family_cohort(fmf)
308 +# (f) robustness rows (Table 19): 12 specifications x {pre, post} on winsorized cohort spreads -> two families (pre, post)
309 +rob = {"baseline": "z_AI_COMPOS", "sector_neutral": "z_AI_COMPOS_sn", "ex_mag7": None, "ex_top500": None, "signal_mean": "z_AI_mean", "signal_top5": "z_AI_top5pct", "iwv_adj": None, "h6": None, "h24": None, "pca": "z_AI_PC1", "seasoned": None, "listed36": None}
310 +for per, yrs in [("pre", PRE), ("post", POST)]:
311 + fam = {}
312 + for nm in rob:
313 + if nm == "ex_mag7": s = cohort_spreads(df[~df.ticker.isin(MAG7)], "z_AI_COMPOS", yrs, "fwd_w")
314 + elif nm == "ex_top500": s = cohort_spreads(df[df.cap_rank > 500], "z_AI_COMPOS", yrs, "fwd_w")
315 + elif nm == "iwv_adj": s = cohort_spreads(df, "z_AI_COMPOS", yrs, "fwd_madj_w")
316 + elif nm == "h6": s = cohort_spreads(df, "z_AI_COMPOS", yrs, "fwd6_w")
317 + elif nm == "h24": s = cohort_spreads(df, "z_AI_COMPOS", range(yrs[0], min(yrs[-1] + 1, 2024)), "fwd24_w")
318 + elif nm == "seasoned": s = cohort_spreads(seasoned, "z_AI_COMPOS", yrs, "fwd_w")
319 + elif nm == "listed36": s = cohort_spreads(l36, "z_AI_COMPOS", yrs, "fwd_w")
320 + else: s = cohort_spreads(df, rob[nm], yrs, "fwd_w")
321 + fam[nm] = [s[y] for y in sorted(s)]
322 + MHT[f"robustness_{per}_w"] = family_cohort(fam)
323 +# deflated Sharpe ratio of the post-launch and full-sample composite LS, trials = every LS series tried in r3k_15 (all signals, samples, groups)
324 +la = RES["ls_alpha"]["AI_COMPOS"]; trials = RES["sr_trials_ann"]
325 +MHT["deflated_sharpe"] = {"full": I.deflated_sharpe(la["full"]["sharpe"], la["full"]["n_months"], la["full"]["skew"], la["full"]["kurt"], trials), "post": I.deflated_sharpe(la["post"]["sharpe"], la["post"]["n_months"], la["post"]["skew"], la["post"]["kurt"], trials), "pre": I.deflated_sharpe(la["pre"]["sharpe"], la["pre"]["n_months"], la["pre"]["skew"], la["pre"]["kurt"], trials), "n_trials": len(trials)}
326 +# count of t > 3 (HLZ) across every statistic printed in the headline families
327 +MHT["hlz_summary"] = {k: int(sum(v["hlz_t_gt_3"])) for k, v in MHT.items() if isinstance(v, dict) and "hlz_t_gt_3" in v}
328 +R2["mht"] = MHT; print("R MHT sorts_full_w RW p:", dict(zip(MHT["sorts_full_w"]["members"], MHT["sorts_full_w"]["p_romano_wolf"])), "| placebo RW", dict(zip(MHT["placebo_interactions_yearcl"]["members"], MHT["placebo_interactions_yearcl"]["p_romano_wolf"])), "| DSR", MHT["deflated_sharpe"], flush=True)
329 +
330 +# ---------------- S. specification curve [N] ----------------
331 +# signals by form x aggregation from filing_scores
332 +agg_cols = {"wmean": "wmean", "mean": "mean", "top5": "top5pct"}; spec_sig = {}
333 +for form in ["all", "10-K", "10-Q", "8-K"]:
334 + sub = fs if form == "all" else fs[fs.form_type == form]
335 + for ag, col in agg_cols.items():
336 + a = sub.groupby(["ticker", "fiscal_year"]).agg(ad=(f"ai_adoption_{col}", "mean"), inf=(f"ai_infrastructure_{col}", "mean")).reset_index(); nm = f"sig_{form}_{ag}"
337 + t_ = df[["ticker", "fiscal_year"]].merge(a, on=["ticker", "fiscal_year"], how="left"); zz = (t_.groupby(t_.fiscal_year)["ad"].transform(zc) + t_.groupby(t_.fiscal_year)["inf"].transform(zc)) / 2; df[nm] = zz.values; spec_sig[(form, ag)] = nm
338 +df["fwd_w5"] = df.groupby("fiscal_year")["fwd_ret"].transform(lambda s: s.clip(s.quantile(.05), s.quantile(.95))); df["fwd6_w5"] = df.groupby("fiscal_year")["fwd_ret_6m"].transform(lambda s: s.clip(s.quantile(.05), s.quantile(.95))); df["fwd24_w5"] = df.groupby("fiscal_year")["fwd_ret_24m"].transform(lambda s: s.clip(s.quantile(.05), s.quantile(.95)))
339 +OCS = {("raw", 6): ("fwd_ret_6m", "mean"), ("raw", 12): ("fwd_ret", "mean"), ("raw", 24): ("fwd_ret_24m", "mean"), ("w199", 6): ("fwd6_w", "mean"), ("w199", 12): ("fwd_w", "mean"), ("w199", 24): ("fwd24_w", "mean"), ("w595", 6): ("fwd6_w5", "mean"), ("w595", 12): ("fwd_w5", "mean"), ("w595", 24): ("fwd24_w5", "mean"), ("median", 6): ("fwd_ret_6m", "median"), ("median", 12): ("fwd_ret", "median"), ("median", 24): ("fwd_ret_24m", "median")}
340 +FILT = {"all": df, "ex_top500": df[df.cap_rank > 500], "large": df[df.size_group == "large"], "small": df[df.size_group == "small"]}
341 +def sort_spread(d, z, weight, oc, stat, yrs):
342 + ls = []
343 + for yr, g in d[d.fiscal_year.isin(list(yrs))].groupby("fiscal_year"):
344 + g = g.dropna(subset=[z, oc]).copy()
345 + if len(g) < 40: continue
346 + g["q"] = pd.qcut(g[z].rank(method="first"), 5, labels=False) + 1; vals = {}
347 + for q in (1, 5):
348 + s = g[g.q == q]
349 + if weight is None: vals[q] = s[oc].median() if stat == "median" else s[oc].mean()
350 + else: w = s[weight].clip(lower=0).fillna(0); vals[q] = np.average(s[oc], weights=w) if w.sum() > 0 else s[oc].mean()
351 + ls.append(vals[5] - vals[1])
352 + mu, t = I.nw_t(ls); return mu, t, I.plain_t(ls), len(ls)
353 +rows = []; t1 = time.time()
354 +for (wz, hz), (oc, st) in OCS.items():
355 + yrs = ALL if hz < 24 else range(2010, 2024)
356 + for wt in ["EW", "VW"]:
357 + for (form, ag), z in spec_sig.items():
358 + for fname, d in FILT.items():
359 + mu, t, tp, n = sort_spread(d, z, None if wt == "EW" else "mktcap", oc, st, yrs); rows.append({"winsor": wz, "horizon": hz, "weight": wt, "form": form, "aggregation": ag, "size_filter": fname, "spread_pct": round(mu * 100, 2), "t_nw": round(t, 2), "t_plain": round(tp, 2), "n_years": n})
360 +sc = pd.DataFrame(rows); sc.to_csv(f"{T}/spec_curve.csv", index=False)
361 +R2["spec_curve"] = {"n_specs": int(len(sc)), "share_positive": round(float((sc.spread_pct > 0).mean()), 3), "share_positive_sig_nw": round(float(((sc.spread_pct > 0) & (sc.t_nw.abs() > 1.96)).mean()), 3), "share_positive_sig_plain": round(float(((sc.spread_pct > 0) & (sc.t_plain.abs() > 1.96)).mean()), 3), "share_negative_sig_nw": round(float(((sc.spread_pct < 0) & (sc.t_nw.abs() > 1.96)).mean()), 3), "median_spread_pct": round(float(sc.spread_pct.median()), 2), "iqr_spread_pct": [round(float(sc.spread_pct.quantile(.25)), 2), round(float(sc.spread_pct.quantile(.75)), 2)], "by_size_filter_median": sc.groupby("size_filter").spread_pct.median().round(2).to_dict(), "by_winsor_median": sc.groupby("winsor").spread_pct.median().round(2).to_dict(), "by_weight_median": sc.groupby("weight").spread_pct.median().round(2).to_dict(), "by_form_median": sc.groupby("form").spread_pct.median().round(2).to_dict(), "seconds": round(time.time() - t1)}
362 +print("S spec curve", R2["spec_curve"], flush=True)
363 +
364 +# ---------------- coverage ----------------
365 +cov = df.groupby("fiscal_year").agg(n_firms=("ticker", "nunique"), median_chunks=("total_chunks", "median"), median_filings=("filing_count", "median"), n_large=("size_group", lambda s: (s == "large").sum()), n_small=("size_group", lambda s: (s == "small").sum()), n_nocap=("mktcap", lambda s: s.isna().sum()), n_large_pct42=("size_group_pct42", lambda s: (s == "large").sum()), n_sp500=("in_sp500", "sum"))
366 +R2["coverage"] = {c: cov[c].to_dict() for c in cov.columns} | {"firms_per_cohort": cov.n_firms.to_dict()}
367 +df["rk"] = df.groupby("fiscal_year")["z_AI_COMPOS"].rank(pct=True); wide = df.pivot_table(index="ticker", columns="fiscal_year", values="rk")
368 +def lag_corr(Lg): return round(float(np.nanmean([wide[y].corr(wide[y + Lg], method="spearman") for y in wide.columns if y + Lg in wide.columns])), 3)
369 +R2["coverage"]["rank_autocorr_lag1"] = lag_corr(1); R2["coverage"]["rank_autocorr_lag5"] = lag_corr(5); R2["runtime_s"] = round(time.time() - T0)
370 +json.dump(R2, open(f"{T}/results2.json", "w"), indent=2, default=str); print("saved results/r3k_v2/results2.json in", R2["runtime_s"], "s", flush=True)
added scripts/r3k_inference.py +146 −0
@@ -0,0 +1,146 @@
1 +"""r3k_inference.py — Stage D (FIRC-0002, package L4, clusters E / M / N) inference helpers, imported by the r3k_15+ scripts.
2 +Everything here is self-contained (numpy / scipy / statsmodels): no remembered critical values — the fixed-b and sup-Wald
3 +reference distributions are simulated with a fixed seed and the simulation parameters are returned with the statistic.
4 + cgm_cov(res, g1, g2) two-way clustered covariance (Cameron-Gelbach-Miller 2011) for a fitted statsmodels OLS
5 + ols_cluster(y, X, groups) OLS with one-way cluster-robust SE (small-sample G/(G-1) x (N-1)/(N-k) adjustment)
6 + wild_cluster_boot(y, X, g, j, ...) restricted wild-cluster bootstrap p-value and bootstrap t draws (Rademacher / Webb)
7 + nw_t(x, lags) Newey-West (Bartlett, `lags`) t of a mean; plain_t(x) its i.i.d. counterpart; ci95
8 + fixed_b_cv(T, bw, ...) Kiefer-Vogelsang (2005) fixed-b critical value for the Bartlett HAC t, simulated
9 + block_boot_ci(x, B, seed) percentile CI of a mean by resampling the series elements (cohort block bootstrap)
10 + mde(se, power=0.8) minimum detectable effect at 5 % two-sided
11 + holm(p), bh(p), romano_wolf(t_obs, t_boot) multiple-testing adjustments (RW: step-down max-t on shared bootstrap draws)
12 + deflated_sharpe(sr, T, skew, kurt, sr_trials) Bailey-Lopez de Prado (2014) DSR probability
13 + permutation_three_cohorts(spreads, post_years) exact permutation p over all C(n,3) subsets and over contiguous blocks
14 + sup_wald(x, trim, nsim, seed) Andrews (1993) sup-Wald mean-shift test with a simulated finite-sample p-value
15 +"""
16 +import itertools, numpy as np, pandas as pd
17 +from scipy import stats
18 +from statsmodels.stats.sandwich_covariance import cov_hac
19 +import statsmodels.api as sm
20 +RNG_SEED = 20261105
21 +
22 +def cgm_cov(res, g1, g2):
23 + v1 = res.get_robustcov_results(cov_type="cluster", groups=np.asarray(g1)).cov_params()
24 + v2 = res.get_robustcov_results(cov_type="cluster", groups=np.asarray(g2)).cov_params()
25 + g12 = pd.factorize(pd.Series(list(zip(np.asarray(g1), np.asarray(g2)))))[0]
26 + v12 = res.get_robustcov_results(cov_type="cluster", groups=g12).cov_params()
27 + return np.asarray(v1) + np.asarray(v2) - np.asarray(v12)
28 +
29 +def ols_cluster(y, X, groups):
30 + y = np.asarray(y, float); X = np.asarray(X, float); g = pd.factorize(np.asarray(groups))[0]; n, k = X.shape; G = g.max() + 1
31 + XtX_inv = np.linalg.pinv(X.T @ X); b = XtX_inv @ X.T @ y; e = y - X @ b
32 + meat = np.zeros((k, k))
33 + for gi in range(G):
34 + m = g == gi; s = X[m].T @ e[m]; meat += np.outer(s, s)
35 + adj = (G / (G - 1)) * ((n - 1) / (n - k)); V = adj * XtX_inv @ meat @ XtX_inv
36 + return b, np.sqrt(np.diag(V)), e, XtX_inv, g
37 +
38 +def wild_cluster_boot(y, X, groups, j, B=9999, weights="rademacher", seed=RNG_SEED):
39 + """Restricted (null-imposed) wild-cluster bootstrap for H0: beta_j = 0 (Cameron-Gelbach-Miller 2008; Webb 2023 six-point
40 + weights for few clusters). Returns dict(p, t_obs, t_boot (B,), G)."""
41 + y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape; g = pd.factorize(np.asarray(groups))[0]; G = g.max() + 1
42 + b, se, e, XtX_inv, _ = ols_cluster(y, X, g); t_obs = b[j] / se[j]
43 + Xr = np.delete(X, j, axis=1); br = np.linalg.pinv(Xr.T @ Xr) @ Xr.T @ y; er = y - Xr @ br; yr = Xr @ br
44 + rng = np.random.default_rng(seed); adj = (G / (G - 1)) * ((n - 1) / (n - k))
45 + blocks = [np.where(g == gi)[0] for gi in range(G)]; Xg = [X[ix] for ix in blocks]
46 + if weights == "webb":
47 + pts = np.array([-np.sqrt(1.5), -1, -np.sqrt(0.5), np.sqrt(0.5), 1, np.sqrt(1.5)]); W = rng.choice(pts, size=(B, G))
48 + else:
49 + W = rng.choice([-1.0, 1.0], size=(B, G))
50 + XtX_inv_Xt = XtX_inv @ X.T; t_boot = np.empty(B)
51 + for bi in range(B):
52 + ys = yr + er * W[bi][g]; bs = XtX_inv_Xt @ ys; es = ys - X @ bs; meat = np.zeros((k, k))
53 + for gi in range(G):
54 + s = Xg[gi].T @ es[blocks[gi]]; meat += np.outer(s, s)
55 + V = adj * XtX_inv @ meat @ XtX_inv; t_boot[bi] = bs[j] / np.sqrt(V[j, j])
56 + p = float((np.abs(t_boot) >= abs(t_obs)).mean()); return {"p": round(p, 4), "t_obs": round(float(t_obs), 2), "t_boot": t_boot, "G": int(G), "B": B, "weights": weights}
57 +
58 +def nw_t(x, lags=3):
59 + x = pd.Series(x).dropna().values
60 + if len(x) < 3: return (float(x.mean()) if len(x) else np.nan), np.nan
61 + r = sm.OLS(x, np.ones((len(x), 1))).fit(); return float(x.mean()), float(x.mean() / np.sqrt(cov_hac(r, nlags=lags)[0, 0]))
62 +def plain_t(x):
63 + x = pd.Series(x).dropna().values
64 + if len(x) < 2: return np.nan
65 + return float(x.mean() / (x.std(ddof=1) / np.sqrt(len(x))))
66 +def ci95(x):
67 + x = pd.Series(x).dropna().values
68 + if len(x) < 2: return [None, None]
69 + se = x.std(ddof=1) / np.sqrt(len(x)); q = stats.t.ppf(0.975, len(x) - 1); return [float(x.mean() - q * se), float(x.mean() + q * se)]
70 +def nw_ci95(x, lags=3):
71 + x = pd.Series(x).dropna().values
72 + if len(x) < 3: return [None, None]
73 + r = sm.OLS(x, np.ones((len(x), 1))).fit(); se = np.sqrt(cov_hac(r, nlags=lags)[0, 0]); return [float(x.mean() - 1.96 * se), float(x.mean() + 1.96 * se)]
74 +
75 +_FB_CACHE = {}
76 +def fixed_b_cv(T, bw, alpha=0.05, nsim=20000, steps=1000, seed=RNG_SEED):
77 + """Kiefer-Vogelsang (2005) fixed-b limit of the Bartlett-kernel HAC t-statistic of a mean, b = (bw+1)/T as in statsmodels
78 + cov_hac(nlags=bw) (Bartlett weights 1 - l/(bw+1)). Limit: W(1)/sqrt(Q_b), Q_b = (2/b) int_0^1 B(r)^2 dr - (2/b) int_0^{1-b} B(r)B(r+b) dr,
79 + B = Brownian bridge. Two-sided critical value at `alpha`, simulated once per (T, bw) with a fixed seed."""
80 + b = (bw + 1) / T; key = (round(b, 4), alpha)
81 + if key in _FB_CACHE: return _FB_CACHE[key]
82 + rng = np.random.default_rng(seed); dt = 1.0 / steps; cvs = []
83 + for chunk in range(nsim // 2000):
84 + dW = rng.normal(0, np.sqrt(dt), size=(2000, steps)); W = np.cumsum(dW, axis=1); r = np.arange(1, steps + 1) / steps
85 + Bb = W - r * W[:, -1:]; lag = int(round(b * steps)); Q = (2 / b) * (Bb ** 2).sum(1) * dt - (2 / b) * (Bb[:, :steps - lag] * Bb[:, lag:]).sum(1) * dt
86 + cvs.append(np.abs(W[:, -1]) / np.sqrt(np.maximum(Q, 1e-12)))
87 + cv = float(np.quantile(np.concatenate(cvs), 1 - alpha)); _FB_CACHE[key] = {"cv": round(cv, 3), "b": round(b, 4), "nsim": nsim, "steps": steps}; return _FB_CACHE[key]
88 +
89 +def block_boot_ci(x, B=9999, seed=RNG_SEED):
90 + x = pd.Series(x).dropna().values; rng = np.random.default_rng(seed)
91 + if len(x) < 2: return [None, None]
92 + m = rng.choice(x, size=(B, len(x)), replace=True).mean(1); return [float(np.quantile(m, 0.025)), float(np.quantile(m, 0.975))]
93 +def mde(se, power=0.8, alpha=0.05): return float((stats.norm.ppf(1 - alpha / 2) + stats.norm.ppf(power)) * se)
94 +
95 +def holm(p):
96 + p = np.asarray(p, float); m = len(p); o = np.argsort(p); adj = np.empty(m); run = 0.0
97 + for i, idx in enumerate(o): run = max(run, (m - i) * p[idx]); adj[idx] = min(1.0, run)
98 + return adj
99 +def bh(p):
100 + p = np.asarray(p, float); m = len(p); o = np.argsort(p)[::-1]; adj = np.empty(m); run = 1.0
101 + for i, idx in enumerate(o): rank = m - i; run = min(run, m / rank * p[idx]); adj[idx] = min(1.0, run)
102 + return adj
103 +def romano_wolf(t_obs, t_boot):
104 + """Step-down max-|t| adjustment (Romano-Wolf 2005) from a (B, m) matrix of bootstrap t's computed under the null with the SAME
105 + draws for all m hypotheses (so cross-sectional dependence is preserved). Returns adjusted p-values."""
106 + t_obs = np.abs(np.asarray(t_obs, float)); tb = np.abs(np.asarray(t_boot, float)); m = len(t_obs); o = np.argsort(-t_obs); adj = np.empty(m); prev = 0.0
107 + for i, idx in enumerate(o):
108 + rest = o[i:]; maxt = tb[:, rest].max(1); p = float((maxt >= t_obs[idx]).mean()); prev = max(prev, p); adj[idx] = prev
109 + return adj
110 +
111 +def deflated_sharpe(sr_ann, T, skew, kurt, sr_trials_ann):
112 + """Bailey & Lopez de Prado (2014). Inputs: annualised SR of the candidate, number of monthly observations T, skewness and
113 + (non-excess) kurtosis of its monthly returns, the annualised SR of every strategy tried (list). Monthly SR = SR_ann / sqrt(12)."""
114 + sr = sr_ann / np.sqrt(12); trials = np.asarray(sr_trials_ann, float) / np.sqrt(12); N = len(trials); gamma = 0.5772156649
115 + v = float(trials.var(ddof=1)) if N > 1 else 0.0
116 + sr0 = np.sqrt(v) * ((1 - gamma) * stats.norm.ppf(1 - 1 / N) + gamma * stats.norm.ppf(1 - 1 / (N * np.e))) if N > 1 else 0.0
117 + denom = np.sqrt(max(1 - skew * sr + (kurt - 1) / 4 * sr ** 2, 1e-9)); z = (sr - sr0) * np.sqrt(T - 1) / denom
118 + return {"dsr_prob": round(float(stats.norm.cdf(z)), 4), "sr0_ann": round(float(sr0 * np.sqrt(12)), 3), "n_trials": int(N), "psr_z": round(float(z), 3)}
119 +
120 +def permutation_three_cohorts(spreads, post_years):
121 + """spreads: dict year -> spread; post_years: the observed post block. Exact permutation over all C(n, k) subsets of k cohorts
122 + (statistic: mean(post) - mean(rest)) and over the n-k+1 contiguous blocks of length k."""
123 + yrs = sorted(spreads); v = np.array([spreads[y] for y in yrs]); k = len(post_years); obs = v[[yrs.index(y) for y in post_years]].mean() - v[[i for i, y in enumerate(yrs) if y not in post_years]].mean()
124 + diffs = []
125 + for comb in itertools.combinations(range(len(yrs)), k):
126 + m = np.zeros(len(yrs), bool); m[list(comb)] = True; diffs.append(v[m].mean() - v[~m].mean())
127 + diffs = np.array(diffs); contig = []
128 + for s in range(len(yrs) - k + 1):
129 + m = np.zeros(len(yrs), bool); m[s:s + k] = True; contig.append((yrs[s], float(v[m].mean() - v[~m].mean())))
130 + cvals = np.array([c[1] for c in contig])
131 + return {"observed_diff": float(obs), "n_subsets": int(len(diffs)), "p_two_sided": float((np.abs(diffs) >= abs(obs) - 1e-12).mean()), "p_one_sided": float((diffs >= obs - 1e-12).mean()), "rank_among_subsets": int((diffs > obs + 1e-12).sum() + 1),
132 + "contiguous_blocks": {str(y): round(d, 4) for y, d in contig}, "p_contiguous_one_sided": float((cvals >= obs - 1e-12).mean()), "rank_among_contiguous": int((cvals > obs + 1e-12).sum() + 1)}
133 +
134 +def sup_wald(x, trim=0.15, nsim=20000, seed=RNG_SEED):
135 + """Andrews (1993) sup-Wald test for a single mean shift at an unknown date in a short series, trimming `trim` at each end;
136 + the Wald statistic uses the full-sample residual variance (homoskedastic, no HAC: T is 15). The finite-sample null distribution
137 + is simulated (i.i.d. normal series of the same length, same trimming), so the p-value needs no tabulated critical value."""
138 + x = np.asarray(pd.Series(x).dropna(), float); T = len(x); lo, hi = max(2, int(np.floor(trim * T))), T - max(2, int(np.floor(trim * T)))
139 + def sw(v):
140 + best, bk = -1.0, None
141 + for k in range(lo, hi + 1):
142 + a, b = v[:k], v[k:]; s2 = (((a - a.mean()) ** 2).sum() + ((b - b.mean()) ** 2).sum()) / (T - 2); w = (a.mean() - b.mean()) ** 2 / (s2 * (1 / k + 1 / (T - k)))
143 + if w > best: best, bk = w, k
144 + return best, bk
145 + obs, k = sw(x); rng = np.random.default_rng(seed); sims = np.array([sw(rng.normal(size=T))[0] for _ in range(nsim)])
146 + return {"sup_wald": round(float(obs), 3), "break_after_index": int(k), "p_simulated": round(float((sims >= obs).mean()), 4), "T": T, "trim": trim, "nsim": nsim, "cv05_simulated": round(float(np.quantile(sims, 0.95)), 3), "cv10_simulated": round(float(np.quantile(sims, 0.90)), 3)}
147
Menu

Theme
↑↓ navigate ↵ open esc close / or ⌘K anywhere