From ab3acfe4cb3386c3d4f17e7391c8fc8a4db59a11 Mon Sep 17 00:00:00 2001 From: Gahow Wang Date: Wed, 30 Sep 2026 20:25:08 +0800 Subject: [PATCH] feat(research): supply-chain links, theme propagation, cycles, attention, small caps Tests of the "rank and buy top N misses the chain / cycle / sentiment" question (report: https://claude.ai/artifact/HUnPtcy7D7Zp1yNQsDLRHY): - review_live_us: Apr-Sep 2026 replay; book return vs AI-chain weight corr 0.71, overnight_recovery_top10 excess almost all allocation (memory stage +28pp). - fetch_supply_links: 10-K mentions of S&P 1500 companies by role, PIT by filing date (12,117 S&P 500 10-Ks; --universe mid for 14,382 S&P 400/600 10-Ks). - lead_lag: corr-peer and 10-K link momentum worked 2006-16 (FM t 3.3-4.3) and vanished 2017-26 (t -2.1..-0.6); Granger lead-lag network and earnings spillover drift null after 2016. - theme_propagation: linked laggards of each era's hot industry underperform 2017-26 (-2% to -12%/yr). - lead_lag_small: the same signals on S&P 400/600, 2019-26, all |t| <= 1.8. - sector_cycle: 21 cycle/rotation variants, none with holdout |t| > 2. - fetch_wiki_attention + attention: Wikipedia pageview attention, |IC| <= 0.013. - concentration: industry caps on momentum top 10 cut alpha, not drawdown. Co-Authored-By: Claude Opus 5.5 --- research/lab/attention.py | 217 ++++++++ research/lab/concentration.py | 88 +++ research/lab/fetch_supply_links.py | 350 ++++++++++++ research/lab/fetch_wiki_attention.py | 285 ++++++++++ research/lab/lead_lag.py | 409 ++++++++++++++ research/lab/lead_lag_small.py | 360 ++++++++++++ research/lab/review_live_us.py | 303 +++++++++++ research/lab/sector_cycle.py | 781 +++++++++++++++++++++++++++ research/lab/theme_propagation.py | 268 +++++++++ 9 files changed, 3061 insertions(+) create mode 100644 research/lab/attention.py create mode 100644 research/lab/concentration.py create mode 100644 research/lab/fetch_supply_links.py create mode 100644 research/lab/fetch_wiki_attention.py create mode 100644 research/lab/lead_lag.py create mode 100644 research/lab/lead_lag_small.py create mode 100644 research/lab/review_live_us.py create mode 100644 research/lab/sector_cycle.py create mode 100644 research/lab/theme_propagation.py diff --git a/research/lab/attention.py b/research/lab/attention.py new file mode 100644 index 0000000..78412e9 --- /dev/null +++ b/research/lab/attention.py @@ -0,0 +1,217 @@ +""" +Retail attention from Wikipedia pageviews: does it predict S&P 500 returns, +and does it improve the momentum sleeve? + + uv run python -m research.lab.fetch_wiki_attention # data/research/wiki_pageviews.pkl + PYTHONPATH=. uv run --with numba python -m research.lab.attention [--views views_main] + +Why Wikipedia: it is the only free, daily, point-in-time attention series that +covers every S&P 500 name back to 2015-07. X full-archive search needs +enterprise/pay-per-use access; Reddit's own API has no history (Pushshift is +moderator-only; Arctic Shift's full-text aggregates time out on +r/wallstreetbets; the academictorrents dumps are the remaining route); +StockTwits closed developer sign-up; Google Trends caps daily data at +~9-month windows with per-query scaling. + +Timing: a pageview day is a UTC day, which ends at 20:00 ET, after the close. +The signal on trading day t uses views through UTC day t-1 only; trades fill +at the open of t+1 (research engine convention). + +Signals (x = log(1 + daily views), article + its redirects): + abn_att mean x over the last 5 days minus the median x of the 252 days before + (Da, Engelberg & Gao 2011 abnormal search volume, on Wikipedia) + att_shock x on the last day minus mean x over the last 21 days + att_trend mean x over the last 63 days minus mean x over the 189 days before + att_level mean x over the last 252 days (popularity; mostly a size proxy) + +Tests, fixed before looking at results: + (a) daily cross-sectional rank IC with the next-open to next-open 5- and + 21-day return; t from non-overlapping samples. + (b) top-20 sleeves, each signal high and low (8 trials): abn_att/att_shock + every 5 days, att_trend/att_level every 21 days, buffer 1. Alpha vs QQQ + (Newey-West), design 2016-07..2020-12, holdout 2021..2026-09. + (c) the momentum sleeve (overnight_mom + rec_126, top 10, every 10 days, + buffer 2 = live overnight_recovery_top10) with 5 attention rules: + veto_hi_abn drop names in the top 10% of abn_att + veto_hi_shock drop names in the top 10% of att_shock + resort_low_abn top 30 by momentum, keep the 10 with the lowest abn_att + resort_high_abn top 30 by momentum, keep the 10 with the highest abn_att + resort_trend top 30 by momentum, keep the 10 with the highest att_trend + A name without pageviews is never vetoed and is neutral in a re-sort. +""" + +from __future__ import annotations + +import dataclasses +import json +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.alpha import _nw_se, alpha_stats, bench_returns +from research.lab.factor_scan import RESULTS, load +from research.lab.report import stats +from research.lab.rotation import factor_ranks +from research.lab.us_sweep import _blend, log_trials, run + +warnings.filterwarnings("ignore") +VIEWS = os.path.join("data", "research", "wiki_pageviews.pkl") +SPAN = ("2016-07-01", "2026-09-24") # first year of views is warm-up +DESIGN = ("2016-07-01", "2020-12-31") +HOLDOUT = ("2021-01-01", "2026-09-24") + + +def load_views(p, which: str = "views") -> pd.DataFrame: + """which = "views" (article + redirects) or "views_main" (article only).""" + with open(VIEWS, "rb") as fh: + d = pickle.load(fh) + print(f"{which}; redirects fetched: {d.get('complete', True)}") + return d[which].reindex(columns=p.close.columns) + + +def signals(p, views: pd.DataFrame) -> dict[str, pd.DataFrame]: + """Calendar-day signals, lagged one day, sampled at each trading close.""" + x = np.log1p(views) + base = x.shift(5).rolling(252, min_periods=200).median() + sig = { + "abn_att": x.rolling(5, min_periods=4).mean() - base, + "att_shock": x - x.rolling(21, min_periods=15).mean(), + "att_trend": x.rolling(63, min_periods=50).mean() - x.shift(63).rolling(189, min_periods=150).mean(), + "att_level": x.rolling(252, min_periods=200).mean(), + } + out = {} + for k, s in sig.items(): + s = s.shift(1) # through UTC day t-1 + s = s.reindex(s.index.union(p.dates)).ffill(limit=3).reindex(p.dates) + out[k] = s.where(p.member) + return out + + +def forward(p, h: int) -> pd.DataFrame: + o = p.open.where(p.open.notna(), p.close.ffill()) + return o.shift(-1 - h) / o.shift(-1) - 1 + + +def ic_table(p, S: dict[str, pd.DataFrame]) -> pd.DataFrame: + rows = [] + for h in (5, 21): + f = forward(p, h).where(p.member) + for k, s in S.items(): + s2, f2 = s.where(f.notna()), f.where(s.notna()) + ic = s2.rank(axis=1).corrwith(f2.rank(axis=1), axis=1, method="pearson") + ic = ic.loc[SPAN[0]:SPAN[1]].dropna() + for label, (a, z) in (("full", SPAN), ("design", DESIGN), ("holdout", HOLDOUT)): + x = ic.loc[a:z] + nov = x.iloc[::h] + rows.append({"signal": k, "h": h, "period": label, "ic": x.mean(), + "t": nov.mean() / nov.std() * np.sqrt(len(nov)), "n": len(nov), + "hit": (x > 0).mean()}) + return pd.DataFrame(rows) + + +def paired_t(d: pd.Series) -> float: + d = d.dropna() + se = _nw_se(np.ones((len(d), 1)), (d - d.mean()).to_numpy())[0] + return float(d.mean() / se) + + +def line(r: pd.Series, b: pd.Series) -> tuple[str, dict]: + rec, cells = {}, [] + for label, (a, z) in (("full", SPAN), ("design", DESIGN), ("holdout", HOLDOUT)): + x = r.loc[a:z] + s = stats((1 + x).cumprod()) + al = alpha_stats(x, b) + rec[label] = {"cagr": s["cagr"], "sharpe": s["sharpe"], "mdd": s["mdd"], + "alpha": al["alpha"], "t": al["t"], "beta": al["beta"]} + cells.append(f"{s['cagr']*100:5.1f}% sh {s['sharpe']:.2f} a {al['alpha']*100:5.1f}% t {al['t']:5.2f}") + return " | ".join(cells), rec + + +def main() -> None: + import argparse + ap = argparse.ArgumentParser() + ap.add_argument("--views", default="views", choices=["views", "views_main"]) + args = ap.parse_args() + tag = "" if args.views == "views" else "_main" + p, _ = load("us") + views = load_views(p, args.views) + b = bench_returns(p.dates) + S = signals(p, views) + cov = S["abn_att"].notna().sum(axis=1) / p.member.sum(axis=1) + print("coverage of PIT members with abn_att:", + {y: round(float(cov.loc[str(y)].mean()), 3) for y in range(2016, 2027)}) + + out = {"coverage": {str(y): float(cov.loc[str(y)].mean()) for y in range(2016, 2027)}} + ic = ic_table(p, S) + print("\n(a) rank IC (next-open fwd return), t from non-overlapping samples") + print(ic.pivot_table(index=["signal", "h"], columns="period", values=["ic", "t"]).round(4).to_string()) + out["ic"] = ic.to_dict("records") + + print(f"\n{'':26s}{'2016-07..2026-09':>40s} | {'design 16-20':>40s} | {'holdout 21-26':>40s}") + R, rows = {}, [] + every = {"abn_att": 5, "att_shock": 5, "att_trend": 21, "att_level": 21} + for k, s in S.items(): + for side, sc in (("high", s), ("low", -s)): + st = run(p, sc.where(p.member), [SPAN], top_n=20, buffer=1.0, every=every[k]) + name = f"{k} {side} top20" + R[name] = st["_rets"] + rows.append({"score": f"attention:{name}", **{kk: v for kk, v in st.items() if not kk.startswith("_")}}) + txt, rec = line(st["_rets"], b) + out[name] = rec + print(f"{name:26s}{txt}", flush=True) + ew = p.close.pct_change(fill_method=None).where(p.member.shift(1, fill_value=False)).mean(axis=1) + txt, out["EW S&P 500"] = line(ew, b) + print(f"{'EW S&P 500 (reference)':26s}{txt}") + + # (c) momentum interaction + pr = factor_ranks(p, "us", ["overnight_mom", "rec_126"]) + mom = _blend(pr["overnight_mom"], pr["rec_126"]).where(p.member) + m = p.member + pct = {k: S[k].rank(axis=1, pct=True) for k in ("abn_att", "att_shock", "att_trend")} + mom_pct = mom.rank(axis=1, pct=True) + rc = {k: float(v.loc[SPAN[0]:].iloc[::21].corrwith(mom_pct.loc[SPAN[0]:].iloc[::21], axis=1).mean()) + for k, v in pct.items()} + print("\nmean cross-sectional rank corr with the momentum score:", {k: round(v, 3) for k, v in rc.items()}) + out["corr_with_momentum"] = rc + top30 = mom.rank(axis=1, ascending=False) <= 30 + ok = lambda veto: m & ~veto.fillna(False) # noqa: E731 + V = { + "base": (mom, m, 2.0), + "veto_hi_abn": (mom, ok(pct["abn_att"] > 0.9), 2.0), + "veto_hi_shock": (mom, ok(pct["att_shock"] > 0.9), 2.0), + "resort_low_abn": ((1 - pct["abn_att"]).fillna(0.5).where(top30), m, 1.0), + "resort_high_abn": (pct["abn_att"].fillna(0.5).where(top30), m, 1.0), + "resort_trend": (pct["att_trend"].fillna(0.5).where(top30), m, 1.0), + } + print(f"\n(c) momentum top10 with attention rules\n{'':26s}{'2016-07..2026-09':>40s} | {'design 16-20':>40s} | {'holdout 21-26':>40s} | vs base") + M = {} + for name, (score, elig, buf) in V.items(): + pv = dataclasses.replace(p, member=elig) + st = run(pv, score.where(elig), [SPAN], top_n=10, buffer=buf, every=10) + M[name] = st["_rets"] + if name != "base": + rows.append({"score": f"attention_mom:{name}", **{kk: v for kk, v in st.items() if not kk.startswith("_")}}) + txt, rec = line(st["_rets"], b) + extra = "" + if name != "base": + d, dd = (M[name] - M["base"]).loc[HOLDOUT[0]:], (M[name] - M["base"]).loc[DESIGN[0]:DESIGN[1]] + extra = (f"holdout {d.mean()*252*100:+5.1f}%/yr t {paired_t(d):+.2f} " + f"(design {dd.mean()*252*100:+5.1f}% t {paired_t(dd):+.2f})") + rec["diff"] = {"holdout": float(d.mean() * 252), "holdout_t": paired_t(d), + "design": float(dd.mean() * 252), "design_t": paired_t(dd)} + out["mom:" + name] = rec + print(f"{name:26s}{txt} | {extra}", flush=True) + if not tag: # the article-only run is a robustness check + log_trials(rows, "attention") + print(f"\ntrials logged: {len(rows)}") + + with open(os.path.join(RESULTS, f"attention{tag}.json"), "w") as fh: + json.dump(out, fh, indent=1, default=float) + pd.DataFrame({**R, **{"mom:" + k: v for k, v in M.items()}}).to_pickle(os.path.join(RESULTS, f"attention{tag}.pkl")) + + +if __name__ == "__main__": + main() diff --git a/research/lab/concentration.py b/research/lab/concentration.py new file mode 100644 index 0000000..261665c --- /dev/null +++ b/research/lab/concentration.py @@ -0,0 +1,88 @@ +""" +Theme concentration in the momentum sleeve. In the 2026-04..09 replay the +live overnight_recovery_top10 held about 92% in AI-hardware names; this +checks whether capping names per industry changes its return or drawdown +over 2006-2026. + + PYTHONPATH=. uv run --with numba python -m research.lab.concentration + +A cap of k per group keeps a stock eligible only while it ranks in the top k +of its group by the same score; the usual top-10 / buffer-2 / every-10-days +selection then runs on the eligible names. Groups: GICS sector (current +classification, 11 groups) and 3-digit SIC (sec_sic.csv, ~150 groups). +Reported: alpha vs QQQ per half, MaxDD, worst 3-month return, and the mean +largest group weight. +""" + +from __future__ import annotations + +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab.alpha import alpha_stats, bench_returns +from research.lab.factor_scan import RESULTS, load +from research.lab.rotation import factor_ranks +from research.lab.us_sweep import _blend, run + +warnings.filterwarnings("ignore") +LONG = [("2006-09-27", "2026-09-24")] + + +def capped(score: pd.DataFrame, groups: pd.Series, k: int) -> pd.DataFrame: + g = groups.reindex(score.columns) + out = score.copy() + for name in g.dropna().unique(): + cols = list(g.index[g == name]) + rk = score[cols].rank(axis=1, ascending=False) + out[cols] = score[cols].where(rk <= k) + return out + + +def top_group_weight(score: pd.DataFrame, groups: pd.Series, n: int = 10) -> float: + rk = score.rank(axis=1, ascending=False) <= n + g = groups.reindex(score.columns) + share = [rk.loc[:, g == x].sum(axis=1) for x in g.dropna().unique()] + m = pd.concat(share, axis=1).max(axis=1) / n + return float(m[rk.sum(axis=1) >= n].mean()) + + +def main() -> None: + p, _ = load("us") + pr = factor_ranks(p, "us", ["overnight_mom", "rec_126"]) + mom = _blend(pr["overnight_mom"], pr["rec_126"]).where(p.member) + sic = pd.read_csv(os.path.join("data", "research", "sec_sic.csv")).drop_duplicates("ticker").set_index("ticker")["sic"] + sic3 = (pd.to_numeric(sic, errors="coerce") // 10).dropna().astype(int).astype(str) + groups = {"gics": p.extra["sector"], "sic3": sic3} + b = bench_returns(p.dates) + variants = [("uncapped", None, None), ("gics <=4", "gics", 4), ("gics <=3", "gics", 3), ("gics <=2", "gics", 2), + ("sic3 <=3", "sic3", 3), ("sic3 <=2", "sic3", 2), ("sic3 <=1", "sic3", 1)] + rows, R = [], {} + for name, gk, k in variants: + s = mom if gk is None else capped(mom, groups[gk], k) + r = run(p, s, LONG, top_n=10, buffer=2.0, every=10)["_rets"] + e = (1 + r).cumprod() + q = (1 + r).rolling(63).apply(np.prod, raw=True) - 1 + row = {"variant": name, "mdd": float((e / e.cummax() - 1).min()), "worst_3m": float(q.min()), + "vol": float(r.std() * np.sqrt(252)), + "max_gics_w": top_group_weight(s.loc["2006-09-27":], groups["gics"]), + "max_sic3_w": top_group_weight(s.loc["2006-09-27":], groups["sic3"])} + for h, x in (("all", r), ("design", r.loc[:"2016"]), ("holdout", r.loc["2017":])): + a = alpha_stats(x, b) + row.update({f"{h}_alpha": a["alpha"], f"{h}_t": a["t"], f"{h}_cagr": a["cagr"]}) + row["sharpe"] = float(r.mean() / r.std() * np.sqrt(252)) + rows.append(row) + R[name] = r + print(f"{name:10s} alpha {row['all_alpha']*100:5.1f}% (t {row['all_t']:.2f}) | 06-16 {row['design_alpha']*100:5.1f}% " + f"t {row['design_t']:.2f} | 17-26 {row['holdout_alpha']*100:5.1f}% t {row['holdout_t']:.2f} | " + f"cagr {row['all_cagr']*100:5.1f}% sharpe {row['sharpe']:.2f} vol {row['vol']*100:4.1f}% " + f"mdd {row['mdd']*100:4.0f}% worst3m {row['worst_3m']*100:4.0f}% | max sector w {row['max_gics_w']*100:3.0f}% " + f"max sic3 w {row['max_sic3_w']*100:3.0f}%", flush=True) + pd.DataFrame(rows).to_csv(os.path.join(RESULTS, "concentration.csv"), index=False) + pd.DataFrame(R).to_pickle(os.path.join(RESULTS, "concentration_returns.pkl")) + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_supply_links.py b/research/lab/fetch_supply_links.py new file mode 100644 index 0000000..7a2cc23 --- /dev/null +++ b/research/lab/fetch_supply_links.py @@ -0,0 +1,350 @@ +""" +Economic links between S&P 500 ever-members from their 10-K filings: which +other member companies each 10-K names, and in what role (customer, supplier, +competitor). The links become known on the filing date, so a signal built +from them is point-in-time. + + PYTHONPATH=. uv run python -m research.lab.fetch_supply_links --workers 4 + PYTHONPATH=. uv run python -m research.lab.fetch_supply_links --universe mid # S&P 400/600 10-Ks from 2015 + +Inputs: data/research/lazy_prices_index.csv (10-K accessions and primary +documents, from research/lab/fetch_lazy_prices.py) and sec_sic.csv (ticker, +CIK, current SEC name). Names: current and former SEC names (submissions JSON, +cached in data/research/sec_names.json), the Wikipedia security name of +current members, and a short manual list of brands (Google, Facebook, AWS). +Matching is case-sensitive on the proper-case or upper-case form of each +alias, with word boundaries. A one-word alias that is an English dictionary +word ("Target", "Ball") is dropped unless listed in ALLOW, and the filer's +own names are skipped. + +A mention is counted in every role whose keywords appear in the same +sentence: customer ("customer", "sales to", "revenue(s) from", "accounted +for ... % of"), supplier ("supplier", "vendor", "supply", "purchase", +"sole source", "foundry", "manufactur-"), competitor ("compet-", "rival"). +In filer i's 10-K, "sales to J accounted for 12% of revenue" marks J as a +customer of i. + +Output: data/research/sec_10k_links.csv, one row per (10-K, named company): +filer, filed, accession, target, n_customer, n_supplier, n_competitor, n_all, +example (first customer sentence, else the first sentence, 300 chars). +data/research/sec_10k_links_done.txt lists processed accessions (resumable). +""" + +from __future__ import annotations + +import argparse +import csv +import json +import multiprocessing as mp +import os +import re +import time + +import pandas as pd + +from research.lab import fetch_lazy_prices as lp + +DIR = os.path.join("data", "research") +OUT = os.path.join(DIR, "sec_10k_links.csv") +DONE = os.path.join(DIR, "sec_10k_links_done.txt") +NAMES = os.path.join(DIR, "sec_names.json") +COLS = ["filer", "filed", "accession", "target", "n_customer", "n_supplier", "n_competitor", "n_all", "example"] + +EXTRA = {"GOOGL": ["Google", "Alphabet"], "GOOG": ["Google", "Alphabet"], "META": ["Facebook", "Meta Platforms"], + "AMZN": ["Amazon", "Amazon.com", "Amazon Web Services", "AWS"], "MSFT": ["Microsoft", "Azure"], + "AAPL": ["Apple"], "HPQ": ["Hewlett-Packard", "HP Inc"], "HPE": ["Hewlett Packard Enterprise"], + "WMT": ["Wal-Mart", "Walmart"], "IBM": ["IBM"], "AMD": ["AMD", "Advanced Micro Devices"], + "NVDA": ["NVIDIA", "Nvidia"], "INTC": ["Intel"], "ORCL": ["Oracle"], "CSCO": ["Cisco"], + "DELL": ["Dell", "Dell Technologies"], "GE": ["General Electric"], "TGT": ["Target Corporation"], + "BA": ["Boeing"], "F": ["Ford Motor"], "GM": ["General Motors"], "TSLA": ["Tesla"], + "CVS": ["CVS"], "UNH": ["UnitedHealth"], "MCK": ["McKesson"], "CAH": ["Cardinal Health"], + "COR": ["AmerisourceBergen", "Cencora"], "HD": ["Home Depot"], "LOW": ["Lowe's", "Lowe’s"], + "COST": ["Costco"], "KR": ["Kroger"], "T": ["AT&T"], "VZ": ["Verizon"], "TMUS": ["T-Mobile"], + "CMCSA": ["Comcast"], "DIS": ["Disney"], "QCOM": ["Qualcomm", "QUALCOMM"], "AVGO": ["Broadcom"], + "TXN": ["Texas Instruments"], "MU": ["Micron"], "WDC": ["Western Digital"], "STX": ["Seagate"], + "LRCX": ["Lam Research"], "AMAT": ["Applied Materials"], "KLAC": ["KLA-Tencor", "KLA Corporation"], + "SNDK": ["SanDisk", "Sandisk"], "GLW": ["Corning Incorporated"], "LMT": ["Lockheed Martin"], + "RTX": ["Raytheon", "United Technologies"], "NOC": ["Northrop Grumman"], "GD": ["General Dynamics"], + "UPS": ["UPS", "United Parcel Service"], "FDX": ["FedEx"], "BBY": ["Best Buy"], "EBAY": ["eBay"], + "MMM": ["3M"], "COHR": ["Coherent Corp", "II-VI"], "BALL": ["Ball Corporation"], + "PGR": ["Progressive Corporation"], "CRM": ["Salesforce", "salesforce.com"], + "NDAQ": ["Nasdaq, Inc", "NASDAQ OMX"], "PFG": ["Principal Financial Group"], "DOW": ["Dow Chemical", "Dow Inc"]} +# aliases that are mostly something else in a 10-K (the exchange, the officer title, the stock index) +DROP = {"NDAQ": {"Nasdaq", "NASDAQ", "Nasdaq Stock Market", "NASDAQ STOCK MARKET"}, + "PFG": {"Principal Financial", "PRINCIPAL FINANCIAL"}, "DOW": {"Dow", "DOW"}} +# sentences about people and financing documents, not business relationships +NOISE = re.compile(r"served as\b|\bserves as\b|\bpositions? (?:with|at)\b|\broles? (?:with|at)\b|\bworked (?:at|for)\b|" + r"\bboard of directors\b|\bprior to joining\b|\bChief \w+ Officer\b|\bVice President\b|\bIndenture\b|" + r"\bCredit Agreement\b|\bRegistration Rights\b|\bunderwrit\w*\b|\btrustee\b|\bas (?:administrative )?agent\b|" + r"\bRule 13a\b|\bCertification\b|\bExhibit\b|\bpar value\b|\bStock Market\b|\bSelect Market\b|" + r"\bDow Jones\b|\bS&P\b", re.I) +ALLOW = {"Apple", "Oracle", "Amazon", "Intel", "Micron", "Boeing", "Tesla", "Google", "Facebook", "Seagate", + "Dell", "Cisco", "Comcast", "Disney", "Verizon", "Qualcomm", "Broadcom", "Microsoft", "Walmart", + "Kroger", "Costco", "Nvidia", "McKesson", "Lumentum", "Corning", "Chevron", "Exxon", "Pfizer", + "Merck", "Amgen", "Nike", "Starbucks", "Visa", "Mastercard", "Netflix", "Adobe", "Salesforce", + "Alphabet"} +SUFFIX = re.compile(r"(?:[\s,]+(?:inc|incorporated|corp|corporation|co|company|companies|ltd|limited|plc|llc|lp|l\.p|" + r"holdings?|group|the|n\.v|s\.a|ag|se|de|del|new|\&))+\.?$", re.I) +SLASH = re.compile(r"\s*/[A-Z]{2,}/?\s*$") +KW = {"customer": re.compile(r"\bcustomers?\b|\bsales to\b|\brevenues? from\b|\baccounted for\b|\bnet sales to\b", re.I), + "supplier": re.compile(r"\bsuppl(?:y|ier|iers|ied)\b|\bvendors?\b|\bpurchas\w*\b|\bsole[- ]source\b|\bfoundr(?:y|ies)\b|" + r"\bmanufactur\w*\b", re.I), + "competitor": re.compile(r"\bcompet\w*\b|\brivals?\b", re.I)} +_SENT = re.compile(r"(?<=[.;!?])\s+(?=[A-Z(\"“])") + + +def _clean_name(n: str) -> str: + n = SLASH.sub("", n.strip()) + prev = None + while prev != n: + prev = n + n = SUFFIX.sub("", n).strip(" ,.") + return n + + +def _variants(n: str) -> set[str]: + out = {n} + if n.isupper() and len(n) > 4: + out.add(" ".join(w.capitalize() if len(w) > 3 or i == 0 else w.lower() if w in {"OF", "AND", "THE"} else w + for i, w in enumerate(n.split()))) + if not n.isupper(): + out.add(n.upper()) + return {v for v in out if len(v) >= 3} + + +def former_names(sic: pd.DataFrame) -> dict[str, list[str]]: + if os.path.exists(NAMES): + return json.load(open(NAMES)) + out = {} + for cik in sic["cik"].dropna().astype(int).unique(): + r = lp._get(f"https://data.sec.gov/submissions/CIK{cik:010d}.json") + if r is None: + continue + d = r.json() + out[str(cik)] = [d.get("name", "")] + [f.get("name", "") for f in d.get("formerNames", [])] + json.dump(out, open(NAMES, "w")) + return out + + +def aliases(sic: pd.DataFrame | None = None, fn: dict | None = None) -> dict[str, set[str]]: + """Name aliases per ticker. Defaults: the S&P 500 ever-members (sec_sic.csv, + sec_names.json); the mid/small run passes the S&P 1500 frames.""" + if sic is None: + sic = pd.read_csv(os.path.join(DIR, "sec_sic.csv")) + if fn is None: + fn = former_names(sic) + wiki = pd.read_csv(os.path.join("data", "us_sectors.csv")).set_index("ticker")["Security"] + words = {w.strip() for w in open("/usr/share/dict/words") if w.strip().islower()} + out: dict[str, set[str]] = {} + for r in sic.itertuples(): + names = list(fn.get(str(int(r.cik)), [])) + [r.name] + if r.ticker in wiki.index: + names.append(wiki[r.ticker]) + al = set() + for n in names: + if isinstance(n, str) and n.strip(): + al |= _variants(_clean_name(n)) + al |= set(EXTRA.get(r.ticker, [])) + al -= DROP.get(r.ticker, set()) + keep = set() + for a in al: + ws = a.split() + if len(a) < 3: + continue + if len(ws) == 1 and (a.lower() in words or a.capitalize().lower() in words) and a not in ALLOW \ + and a.capitalize() not in ALLOW: + continue + if len(ws) == 1 and a.isupper() and len(a) <= 3 and a not in {"IBM", "AMD", "UPS", "CVS", "AWS", "HP"}: + continue + keep.add(a) + out[r.ticker] = keep + # an alias shared by two tickers (share classes, reused names) stays only if they are one company + owner: dict[str, set[str]] = {} + for t, al in out.items(): + for a in al: + owner.setdefault(a, set()).add(t) + same_co = [{"GOOGL", "GOOG"}, {"FOXA", "FOX"}, {"NWSA", "NWS"}, {"BRK-B", "BRK.B"}] + for a, ts in owner.items(): + if len(ts) > 1 and not any(ts <= s for s in same_co): + for t in ts: + out[t].discard(a) + return out + + +_RX = None +_ALIAS2T: dict[str, str] = {} + + +def _compile(al: dict[str, set[str]]): + global _RX, _ALIAS2T + _ALIAS2T = {a: t for t, s in al.items() for a in s} + alts = sorted(_ALIAS2T, key=len, reverse=True) + _RX = re.compile(r"(? dict[str, dict]: + found: dict[str, dict] = {} + for line in text.split("\n"): + if len(line) < 20 or not _RX.search(line): + continue + for s in _SENT.split(line): + hits = {_ALIAS2T[m.group(0)] for m in _RX.finditer(s)} - own + if not hits or NOISE.search(s): + continue + roles = {k for k, rx in KW.items() if rx.search(s)} + if "competitor" in roles: # "our competitors include X ... customers" is not a customer link + roles = {"competitor"} + for t in hits: + d = found.setdefault(t, {"n_customer": 0, "n_supplier": 0, "n_competitor": 0, "n_all": 0, "example": ""}) + d["n_all"] += 1 + for k in roles: + d[f"n_{k}"] += 1 + if (not d["example"]) or ("customer" in roles and d["n_customer"] == 1): + d["example"] = s[:300] + return found + + +def work(job) -> tuple[str, list[dict]] | None: + filer, own, cik, filed, acc, doc = job + a = acc.replace("-", "") + name = doc if isinstance(doc, str) and doc.strip() else f"{acc}.txt" + r = lp._get(f"https://www.sec.gov/Archives/edgar/data/{cik}/{a}/{name}") + if r is None: + return acc, [] + if name.lower().endswith(".pdf"): + return acc, [] + try: + text = lp.clean_text(r.content) + except Exception: # noqa: BLE001 + return acc, [] + rows = [{"filer": filer, "filed": filed, "accession": acc, "target": t, **d} + for t, d in links_in(text, set(own)).items()] + return acc, rows + + +def _init(lock, nxt, al): + lp._init_worker(lock, nxt) + _compile(al) + + +# --------------------------------------------------------------------------- +# S&P 400 + 600 (research/lab/panel_mid.py): 10-Ks filed from 2015 on, names +# matched against the S&P 1500 (the S&P 500 ever-members plus these) +# --------------------------------------------------------------------------- + +MID_INDEX = os.path.join(DIR, "mid_10k_index.csv") +MID_SIC = os.path.join(DIR, "mid_sec_sic.csv") +MID_NAMES = os.path.join(DIR, "mid_sec_names.json") +MID_OUT = os.path.join(DIR, "mid_10k_links.csv") +MID_DONE = os.path.join(DIR, "mid_10k_links_done.txt") +MID_START = "2015-01-01" + + +def build_mid_index() -> None: + """CIKs of the mid/small panel's tickers (company_tickers.json, i.e. names + still listed today, plus sec_sic.csv), their SEC names and 10-K filings.""" + from research.lab.fetch_us_fundamentals import cik_map + from research.lab.panel_mid import load_us_mid + tickers = [t for t in load_us_mid("midsmall").close.columns if t not in lp.US_ETFS] + cm = cik_map() + old = pd.read_csv(os.path.join(DIR, "sec_sic.csv")).set_index("ticker")["cik"].to_dict() + tc = {t: int(old.get(t, cm.get(t, 0))) for t in tickers} + tc = {t: c for t, c in tc.items() if c} + print(f"{len(tc)}/{len(tickers)} mid/small tickers map to a CIK", flush=True) + by_cik: dict[int, list[str]] = {} + for t, c in tc.items(): + by_cik.setdefault(c, []).append(t) + sic_rows, names, idx = [], {}, [] + t0 = time.time() + for k, (cik, ts) in enumerate(by_cik.items(), 1): + r = lp._get(f"https://data.sec.gov/submissions/CIK{cik:010d}.json") + if r is None: + continue + js = r.json() + names[str(cik)] = [js.get("name", "")] + [x.get("name", "") for x in js.get("formerNames", [])] + parts = [lp._columns(js["filings"]["recent"])] + for fl in js["filings"].get("files", []): + if fl.get("filingTo", "9999") < MID_START: + continue + r2 = lp._get(f"https://data.sec.gov/submissions/{fl['name']}") + if r2 is not None: + parts.append(lp._columns(r2.json())) + df = pd.concat(parts, ignore_index=True).drop_duplicates("accessionNumber") + for t in sorted(ts): + sic_rows.append({"ticker": t, "cik": cik, "sic": js.get("sic", ""), + "sicDescription": js.get("sicDescription", ""), "name": js.get("name", "")}) + fk = df[df["form"].isin(["10-K", "10-K405"]) & (df["filingDate"] >= MID_START)] + idx.append(pd.DataFrame({"cik": cik, "tickers": "|".join(sorted(ts)), "form": fk["form"], + "filed": fk["filingDate"], "period": fk["reportDate"], + "accession": fk["accessionNumber"], "doc": fk["primaryDocument"]})) + if k % 100 == 0: + print(f" submissions {k}/{len(by_cik)} [{time.time() - t0:.0f}s]", flush=True) + pd.DataFrame(sic_rows).sort_values("ticker").to_csv(MID_SIC, index=False) + json.dump(names, open(MID_NAMES, "w")) + ix = pd.concat(idx, ignore_index=True).sort_values(["cik", "filed"]) + ix.to_csv(MID_INDEX, index=False) + print(f"mid index: {len(ix)} 10-Ks for {ix['cik'].nunique()} CIKs [{time.time() - t0:.0f}s]", flush=True) + + +def sp1500_frames() -> tuple[pd.DataFrame, dict]: + big = pd.read_csv(os.path.join(DIR, "sec_sic.csv")) + mid = pd.read_csv(MID_SIC) + sic = pd.concat([big, mid[~mid["ticker"].isin(big["ticker"])]], ignore_index=True) + fn = {**json.load(open(MID_NAMES)), **json.load(open(NAMES))} + return sic, fn + + +def main() -> None: + global OUT, DONE + ap = argparse.ArgumentParser() + ap.add_argument("--workers", type=int, default=4) + ap.add_argument("--limit", type=int, default=0) + ap.add_argument("--universe", default="sp500", choices=["sp500", "mid"]) + args = ap.parse_args() + lock, nxt = mp.Lock(), mp.Value("d", 0.0) + lp._init_worker(lock, nxt) + if args.universe == "mid": + if not os.path.exists(MID_INDEX): + build_mid_index() + OUT, DONE = MID_OUT, MID_DONE + sic_all, fn = sp1500_frames() + al = aliases(sic_all, fn) + ix = pd.read_csv(MID_INDEX, dtype={"doc": str}) + sic = pd.read_csv(MID_SIC) + else: + al = aliases() + ix = pd.read_csv(os.path.join(DIR, "lazy_prices_index.csv"), dtype={"doc": str}) + sic = pd.read_csv(os.path.join(DIR, "sec_sic.csv")) + print(f"{len(al)} tickers, {sum(map(len, al.values()))} aliases", flush=True) + ix = ix[ix["form"].isin(["10-K", "10-K405"])] + cik2t = sic.groupby("cik")["ticker"].apply(list).to_dict() + done = set(open(DONE).read().split()) if os.path.exists(DONE) else set() + jobs = [] + for r in ix.itertuples(): + if r.accession in done: + continue + ts = cik2t.get(r.cik, str(r.tickers).split("|")) + jobs.append((ts[0], [t for t in ts], int(r.cik), r.filed, r.accession, r.doc)) + if args.limit: + jobs = jobs[:args.limit] + print(f"{len(jobs)} 10-Ks to scan ({len(done)} done)", flush=True) + if not os.path.exists(OUT): + with open(OUT, "w", newline="") as fh: + csv.DictWriter(fh, fieldnames=COLS).writeheader() + t0 = time.time() + with mp.get_context("fork").Pool(args.workers, initializer=_init, initargs=(lock, nxt, al)) as pool, \ + open(OUT, "a", newline="") as out, open(DONE, "a") as dn: + w = csv.DictWriter(out, fieldnames=COLS) + for k, (acc, rows) in enumerate(pool.imap_unordered(work, jobs, chunksize=4), 1): + for row in rows: + w.writerow(row) + dn.write(acc + "\n") + if k % 200 == 0 or k == len(jobs): + out.flush(); dn.flush() + el = time.time() - t0 + print(f" {k}/{len(jobs)} [{el/60:.0f} min, ETA {(len(jobs)-k)*el/k/60:.0f} min]", flush=True) + print("done", flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_wiki_attention.py b/research/lab/fetch_wiki_attention.py new file mode 100644 index 0000000..193da1e --- /dev/null +++ b/research/lab/fetch_wiki_attention.py @@ -0,0 +1,285 @@ +""" +Daily English Wikipedia pageviews for every US research-panel ticker, as a +point-in-time attention proxy (the free source with the longest daily history; +X/Reddit history is not available for free, see research/lab/attention.py). + + uv run python -m research.lab.fetch_wiki_attention # mapping + pageviews + uv run python -m research.lab.fetch_wiki_attention --map-only + +Ticker -> article, first hit wins: + 0. OVERRIDES below (hand-checked) + 1. "List of S&P 500 companies" constituents table (current members) + 2. "Historical components of the S&P 500" change log (added/removed names, + most recent entry first, so a reused ticker gets its latest S&P company) + 3. Wikidata: SEC CIK (P5531) of data/research/sec_sic.csv + 4. Wikidata: NYSE/Nasdaq ticker qualifier (P414 / P249) +Titles are resolved to the canonical article; views are summed over the +article and every redirect to it, so a rename (Facebook -> Meta Platforms) +does not break the series. Agent "user" (spiders and automated excluded). + +Writes data/research/wiki_pageviews.pkl: + {"views": DataFrame dates x tickers (canonical + redirects), "complete": bool, + "views_main": same, canonical title only, + "map": DataFrame ticker/title/source, "fetched": timestamp} +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import os +import pickle +import re +import threading +import time +import urllib.parse +from concurrent.futures import ThreadPoolExecutor + +import pandas as pd +import requests +from bs4 import BeautifulSoup + +UA = "quant-research-attention/0.1 (academic factor research; gahow@example.com)" +OUT = os.path.join("data", "research", "wiki_pageviews.pkl") +CACHE = os.path.join("data", "research", "wiki_cache") +START, END = "20150701", "20260929" +SPARQL = "https://query.wikidata.org/sparql" +API = "https://en.wikipedia.org/w/api.php" +PV = ("https://wikimedia.org/api/rest_v1/metrics/pageviews/per-article/en.wikipedia/" + "all-access/user/{t}/daily/{a}/{b}") +# change-log links that point at the wrong article, or names missing from both tables +OVERRIDES = {"ITT": "ITT Inc.", "CNX": "CNX Resources", "PODD": "Insulet", "MBC": "MasterBrand"} +TICKER_RE = re.compile(r"^[A-Z][A-Z0-9.\-]{0,6}$") + +_session = threading.local() + + +def _s() -> requests.Session: + if not hasattr(_session, "s"): + _session.s = requests.Session() + _session.s.headers["User-Agent"] = UA + return _session.s + + +def _get(url, **kw): + for i in range(20): + try: + r = _s().get(url, timeout=60, **kw) + except requests.RequestException: + time.sleep(2 ** i) + continue + if r.status_code == 429 or r.status_code >= 500: + ra = r.headers.get("Retry-After") + time.sleep(float(ra) + 1 if ra else 2 ** i) + continue + return r + raise RuntimeError(f"failed: {url}") + + +def _norm(t: str) -> str: + return t.strip().upper().replace(".", "-") + + +def _title(href: str) -> str | None: + if not href or not href.startswith("./"): + return None + t = urllib.parse.unquote(href[2:]).split("#")[0] + return t.replace("_", " ") + + +def _page(title: str) -> BeautifulSoup: + r = _get("https://en.wikipedia.org/api/rest_v1/page/html/" + urllib.parse.quote(title, safe="")) + return BeautifulSoup(r.text, "lxml") + + +def _from_wikipedia() -> tuple[dict, dict]: + cur = {} + s = _page("List_of_S&P_500_companies") + for tr in s.find("table", id="constituents").find_all("tr")[1:]: + td = tr.find_all("td") + if len(td) >= 2: + a = td[1].find("a") + t = _title(a.get("href")) if a else None + if t: + cur[_norm(td[0].get_text(strip=True))] = t + hist = {} + s = _page("Historical_components_of_the_S&P_500") + for tr in s.find("table", id="changes").find_all("tr")[2:]: # newest first + cells = tr.find_all("td") + for c, nxt in zip(cells, cells[1:]): + tk = c.get_text(strip=True) + a = nxt.find("a") + if TICKER_RE.match(tk) and a is not None: + t = _title(a.get("href")) + if t: + hist.setdefault(_norm(tk), t) + return cur, hist + + +def _sparql(q: str) -> pd.DataFrame: + r = _get(SPARQL, params={"query": q}, headers={"Accept": "text/csv"}) + from io import StringIO + return pd.read_csv(StringIO(r.text), dtype=str) + + +def _from_wikidata() -> tuple[dict, dict]: + c = _sparql("""SELECT ?cik ?article WHERE { ?item wdt:P5531 ?cik . + ?article schema:about ?item ; schema:isPartOf . }""") + by_cik = {} + for cik, art in zip(c["cik"].str.lstrip("0"), c["article"]): + by_cik.setdefault(cik, urllib.parse.unquote(art.rsplit("/", 1)[1]).replace("_", " ")) + t = _sparql("""SELECT ?ticker ?article WHERE { VALUES ?exch { wd:Q13677 wd:Q82059 } + ?item p:P414 ?st . ?st ps:P414 ?exch ; pq:P249 ?ticker . + ?article schema:about ?item ; schema:isPartOf . }""") + by_tk = {} + for tk, art in zip(t["ticker"], t["article"]): + by_tk.setdefault(_norm(tk), urllib.parse.unquote(art.rsplit("/", 1)[1]).replace("_", " ")) + return by_cik, by_tk + + +def _canonical(titles: list[str]) -> dict[str, str]: + out = {} + for i in range(0, len(titles), 50): + batch = titles[i:i + 50] + r = _get(API, params={"action": "query", "titles": "|".join(batch), "redirects": 1, + "format": "json", "formatversion": 2}).json()["query"] + norm = {n["from"]: n["to"] for n in r.get("normalized", [])} + red = {n["from"]: n["to"] for n in r.get("redirects", [])} + missing = {p["title"] for p in r["pages"] if p.get("missing")} + for t in batch: + x = norm.get(t, t) + x = red.get(x, x) + if x not in missing: + out[t] = x + return out + + +def _redirects(titles: list[str]) -> dict[str, list[str]]: + """Redirects into each title, 50 titles per request, sequential (the action + API throttles parallel clients); cached in CACHE/redirects.json.""" + path = os.path.join(CACHE, "redirects.json") + os.makedirs(CACHE, exist_ok=True) + out = json.load(open(path)) if os.path.exists(path) else {} + todo = [t for t in titles if t not in out] + for i in range(0, len(todo), 50): + batch, cont = todo[i:i + 50], {} + got = {t: [] for t in batch} + while True: + r = _get(API, params={"action": "query", "titles": "|".join(batch), "prop": "redirects", + "rdnamespace": 0, "rdlimit": "max", "format": "json", + "formatversion": 2, **cont}).json() + for p in r["query"]["pages"]: + got.setdefault(p["title"], []).extend(x["title"] for x in p.get("redirects", [])) + if "continue" not in r: + break + cont = r["continue"] + out.update({t: got.get(t, []) for t in batch}) + with open(path, "w") as fh: + json.dump(out, fh) + return {t: out[t] for t in titles} + + +def build_map(tickers: list[str]) -> pd.DataFrame: + cur, hist = _from_wikipedia() + by_cik, by_tk = _from_wikidata() + sic = pd.read_csv(os.path.join("data", "research", "sec_sic.csv"), dtype={"cik": str}) + cik = dict(zip(sic["ticker"], sic["cik"].str.lstrip("0"))) + rows = [] + for t in tickers: + n = _norm(t) + for src, title in (("override", OVERRIDES.get(t)), ("sp500_list", cur.get(n)), ("sp500_changes", hist.get(n)), + ("wikidata_cik", by_cik.get(cik.get(t, ""))), + ("wikidata_ticker", by_tk.get(n))): + if title: + rows.append({"ticker": t, "title": title, "source": src}) + break + else: + rows.append({"ticker": t, "title": None, "source": None}) + m = pd.DataFrame(rows) + canon = _canonical(sorted(m["title"].dropna().unique())) + m["canonical"] = m["title"].map(canon) + return m + + +def _views(title: str) -> pd.Series: + os.makedirs(CACHE, exist_ok=True) + key = os.path.join(CACHE, re.sub(r"[^A-Za-z0-9]+", "_", title)[:120] + "_" + hashlib.md5(title.encode()).hexdigest()[:8] + ".json") + if os.path.exists(key): + with open(key) as fh: + d = json.load(fh) + else: + url = PV.format(t=urllib.parse.quote(title.replace(" ", "_"), safe=""), a=START, b=END) + r = _get(url) + d = r.json().get("items", []) if r.status_code == 200 else [] + with open(key, "w") as fh: + json.dump(d, fh) + time.sleep(0.05) + if not d: + return pd.Series(dtype=float) + return pd.Series({pd.Timestamp(x["timestamp"][:8]): float(x["views"]) for x in d}) + + +def _assemble(m, titles, reds, V, with_redirects: bool) -> tuple[pd.DataFrame, pd.DataFrame]: + idx = pd.date_range(pd.Timestamp(START), pd.Timestamp(END), freq="D") + main, tot = {}, {} + for t in titles: + s0 = V[t].reindex(idx) + parts = [s0.fillna(0.0)] + ([V[r].reindex(idx).fillna(0.0) for r in reds[t]] if with_redirects else []) + s = pd.concat(parts, axis=1).sum(axis=1) + main[t], tot[t] = s0, s.where(s0.notna() | (s > 0)) + tick = m.dropna(subset=["canonical"]) + views = pd.DataFrame({r.ticker: tot[r.canonical] for r in tick.itertuples()}) + views_main = pd.DataFrame({r.ticker: main[r.canonical] for r in tick.itertuples()}) + return views, views_main + + +def _save(m, views, views_main, complete: bool) -> None: + with open(OUT, "wb") as fh: + pickle.dump({"views": views, "views_main": views_main, "map": m, "complete": complete, + "fetched": pd.Timestamp.now()}, fh) + print(f"wrote {OUT} (redirects {'included' if complete else 'NOT included'}): {views.shape}", flush=True) + + +def fetch_views(m: pd.DataFrame, workers: int = 4) -> None: + """Canonical articles first (saved without redirects), then every redirect.""" + titles = sorted(m["canonical"].dropna().unique()) + reds = _redirects(titles) + rest = sorted({r for v in reds.values() for r in v} - set(titles)) + print(f"{len(titles)} articles, {len(rest)} redirects", flush=True) + done = [0] + + def one(t): + s = _views(t) + done[0] += 1 + if done[0] % 500 == 0: + print(f" {done[0]}/{len(titles) + len(rest)}", flush=True) + return t, s + + with ThreadPoolExecutor(workers) as ex: + V = dict(ex.map(one, titles)) + _save(m, *_assemble(m, titles, reds, V, False), complete=False) + with ThreadPoolExecutor(workers) as ex: + V.update(dict(ex.map(one, rest))) + _save(m, *_assemble(m, titles, reds, V, True), complete=True) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--map-only", action="store_true") + ap.add_argument("--workers", type=int, default=4) + args = ap.parse_args() + with open(os.path.join("data", "research", "us.panel.pkl"), "rb") as fh: + tickers = list(pickle.load(fh).close.columns) + m = build_map(tickers) + print(m["source"].value_counts(dropna=False).to_string()) + print(f"mapped {m['canonical'].notna().sum()}/{len(m)}; shared articles:", + m.dropna(subset=["canonical"]).groupby("canonical")["ticker"].agg(list).loc[lambda s: s.str.len() > 1].to_dict()) + m.to_csv(os.path.join("data", "research", "wiki_map.csv"), index=False) + if args.map_only: + return + fetch_views(m, args.workers) + + +if __name__ == "__main__": + main() diff --git a/research/lab/lead_lag.py b/research/lab/lead_lag.py new file mode 100644 index 0000000..15a34eb --- /dev/null +++ b/research/lab/lead_lag.py @@ -0,0 +1,409 @@ +""" +Does a linked company's move predict a stock that has not moved yet? Tests +of the "X already rose, its supplier Y has not priced it in" idea on PIT +S&P 500 members, with links of four kinds: + + corr_peer the 10 members whose daily residual returns (market beta + removed) correlated most with the stock over the last 252 + days; re-estimated each month end. A proxy for economic + similarity (same industry, same customers, same inputs). + lead_peer the 10 members whose weekly residual return in week w-1 + best predicted the stock's return in week w over the last + 104 weeks (lagged cross-correlation); re-estimated each + month end. A data-driven lead-lag network. + 10k_customer members named in the stock's latest 10-K in a customer + sentence (research/lab/fetch_supply_links.py), valid from + the filing day for 15 months. Cohen & Frazzini (2008) + customer momentum. + 10k_supplier members named in a supplier sentence (upstream links). + 10k_any any member named in the 10-K, either direction (a filer's + links plus the filers that name it). + +Signals, known at the close of day t: + _mom mean return of the linked names over the last 21 trading + days (5 for lead_peer). + _gap that minus the stock's own return over the same days: the + "linked names moved, this one has not" catch-up signal. + +Measures (design 2006-2016, holdout 2017-2026, 10 bp, next-open fills): + * Fama-MacBeth: next-21-day (next-5-day for lead_peer) return rank on the + signal rank, controlling for the stock's own 1-month return, 12-1 + momentum and overnight momentum ranks; non-overlapping periods. + * a top-20 sleeve on the signal alone: alpha against QQQ. + * momentum top 10 (the live overnight_recovery_top10) with the top 30 by + momentum re-sorted by the signal, vs momentum alone. +Plus an event test: after a member's earnings day (8-K item 2.02), do its +corr_peers drift in the direction of its announcement return over the next +1-21 days (under-reaction, the idea) or reverse (over-reaction, Thomas & +Zhang 2008)? + + PYTHONPATH=. uv run --with numba python -m research.lab.lead_lag [--part signals|events|chain] +""" + +from __future__ import annotations + +import argparse +import json +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factors import _ret, resid_returns +from research.lab.factor_scan import RESULTS, load +from research.lab.rotation import factor_ranks + +warnings.filterwarnings("ignore") +LONG = [("2006-09-27", "2026-09-24")] +DESIGN_END = "2016-12-31" +K = 10 +CACHE = os.path.join("data", "research", "lead_lag_links.pkl") + + +# --------------------------------------------------------------------------- +# link matrices: {month_end: (tickers_idx, W)} with W[i, j] = weight of j in i's links +# --------------------------------------------------------------------------- + +def month_ends(dates: pd.DatetimeIndex) -> list[pd.Timestamp]: + s = pd.Series(dates, index=dates) + return list(s.groupby(dates.to_period("M")).max()) + + +def corr_peer_links(p, win: int = 252, k: int = K) -> dict: + e = resid_returns(p).where(p.member) + out = {} + X = e.to_numpy() + for me in month_ends(p.dates): + i1 = p.dates.get_loc(me) + 1 + if i1 < win: + continue + blk = X[i1 - win:i1] + mem = p.member.iloc[i1 - 1].to_numpy() & (np.isfinite(blk).sum(axis=0) >= int(win * 0.8)) + idx = np.where(mem)[0] + Z = blk[:, idx] + Z = (Z - np.nanmean(Z, axis=0)) / np.nanstd(Z, axis=0) + Z = np.nan_to_num(Z) + C = Z.T @ Z / win + np.fill_diagonal(C, -np.inf) + top = np.argsort(-C, axis=1)[:, :k] + out[me] = (idx, idx[top]) + return out + + +def lead_peer_links(p, weeks: int = 104, k: int = K) -> dict: + e = resid_returns(p).where(p.member).fillna(0.0) + cum = e.cumsum().to_numpy() + out = {} + for me in month_ends(p.dates): + i1 = p.dates.get_loc(me) + if i1 < 5 * (weeks + 1) + 252: + continue + pts = np.arange(i1, i1 - 5 * (weeks + 1), -5)[::-1] + W = np.diff(cum[pts], axis=0) # (weeks, N) weekly residual returns + mem = p.member.iloc[i1].to_numpy() & (np.abs(W).sum(axis=0) > 0) + idx = np.where(mem)[0] + A, B = W[:-1, idx], W[1:, idx] + A = (A - A.mean(0)) / (A.std(0) + 1e-12) + B = (B - B.mean(0)) / (B.std(0) + 1e-12) + L = A.T @ B / len(A) # L[j, i]: j's week w-1 vs i's week w + np.fill_diagonal(L, -np.inf) + top = np.argsort(-L.T, axis=1)[:, :k] # for each i, its k best leaders + out[me] = (idx, idx[top]) + return out + + +def tenk_links(p, role: str) -> pd.DataFrame: + """Rows of (filer, target, filed) for the role ("customer", "supplier", + "any"); "any" is made symmetric.""" + d = pd.read_csv(os.path.join("data", "research", "sec_10k_links.csv"), parse_dates=["filed"]) + if role != "any": + d = d[d[f"n_{role}"] > 0] + d = d[["filer", "target", "filed"]] + if role == "any": + d = pd.concat([d, d.rename(columns={"filer": "target", "target": "filer"})]) + return d + + +def tenk_link_sets(p, role: str, valid_days: int = 456) -> dict: + """Month-end link sets from 10-Ks filed within the last 15 months (by + filing date, so point-in-time).""" + d = tenk_links(p, role) + col = {t: i for i, t in enumerate(p.close.columns)} + d = d[d["filer"].isin(col) & d["target"].isin(col) & (d["filer"] != d["target"])] + out = {} + for me in month_ends(p.dates): + live = d[(d["filed"] <= me) & (d["filed"] > me - pd.Timedelta(days=valid_days))] + if live.empty: + continue + mem = p.member.loc[me] + g = {} + for f, t in zip(live["filer"], live["target"]): + if mem.get(f, False) and mem.get(t, False): + g.setdefault(col[f], set()).add(col[t]) + if g: + out[me] = {i: np.array(sorted(s)) for i, s in g.items()} + return out + + +def link_signal(p, links: dict, horizon: int, kind: str) -> pd.DataFrame: + """Mean linked-name log return over the last ``horizon`` days (kind="mom") + or that minus the stock's own (kind="gap"), using each month's links for + the days until the next month end.""" + lr = np.log1p(_ret(p).fillna(0.0)) + R = lr.rolling(horizon, min_periods=horizon).sum().to_numpy() + out = np.full(R.shape, np.nan) + mes = sorted(links) + for a, b in zip(mes, mes[1:] + [p.dates[-1] + pd.Timedelta(days=1)]): + rows = np.where((p.dates > a) & (p.dates <= b))[0] + if not len(rows): + continue + L = links[a] + if isinstance(L, tuple): + idx, nb = L + val = np.nanmean(R[np.ix_(rows, nb.ravel())].reshape(len(rows), *nb.shape), axis=2) + out[np.ix_(rows, idx)] = val + else: + for i, nb in L.items(): + out[rows, i] = np.nanmean(R[np.ix_(rows, nb)], axis=1) + S = pd.DataFrame(out, index=p.dates, columns=p.close.columns) + if kind == "gap": + S = S - pd.DataFrame(R, index=p.dates, columns=p.close.columns) + return S.where(p.member) + + +# --------------------------------------------------------------------------- +# measures +# --------------------------------------------------------------------------- + +def fama_macbeth(p, sig: pd.DataFrame, controls: dict[str, pd.DataFrame], horizon: int) -> dict: + """Cross-sectional regressions every ``horizon`` days of the forward + return rank on the signal rank and control ranks (all in [-0.5, 0.5]).""" + o = p.open.where(p.open.notna(), p.close).ffill() + fwd = (o.shift(-(horizon + 1)) / o.shift(-1) - 1).where(p.member) # next open to the open h days later + rows = [] + for i in range(260, len(p.dates) - horizon, horizon): + d = p.dates[i] + y = fwd.iloc[i] + cols = {"sig": sig.iloc[i], **{k: v.iloc[i] for k, v in controls.items()}} + df = pd.DataFrame({"y": y, **cols}).dropna() + if len(df) < 50: + continue + df = df.rank(pct=True) - 0.5 + X = np.column_stack([np.ones(len(df)), df.drop(columns="y").to_numpy()]) + b, *_ = np.linalg.lstsq(X, df["y"].to_numpy(), rcond=None) + ic = df["y"].corr(df["sig"]) + rows.append({"date": d, "b_sig": b[1], "ic": ic, "n": len(df)}) + fm = pd.DataFrame(rows).set_index("date") + + def tstat(x): + x = x.dropna() + return float(x.mean() / x.std() * np.sqrt(len(x))) if len(x) > 2 else np.nan + out = {} + for name, sub in (("all", fm), ("design", fm.loc[:DESIGN_END]), ("holdout", fm.loc["2017":])): + out[name] = {"b": float(sub["b_sig"].mean()), "t": tstat(sub["b_sig"]), + "ic": float(sub["ic"].mean()), "ic_t": tstat(sub["ic"]), "n_obs": float(sub["n"].mean()), + "periods": len(sub)} + return out + + +def sleeve_stats(p, score: pd.DataFrame, top_n: int, every: int, bench) -> dict: + from research.lab.alpha import alpha_stats + from research.lab.us_sweep import run + r = run(p, score.where(p.member), LONG, top_n=top_n, buffer=2.0, every=every)["_rets"] + e = (1 + r).cumprod() + out = {"mdd": float((e / e.cummax() - 1).min())} + for name, x in (("all", r), ("design", r.loc[:DESIGN_END]), ("holdout", r.loc["2017":])): + a = alpha_stats(x, bench) + out[name] = {"alpha": a["alpha"], "t": a["t"], "cagr": a["cagr"], "bench": a["bench_cagr"], "beta": a["beta"]} + return out, r + + +def resort(mom: pd.DataFrame, sig: pd.DataFrame, pool: int = 30) -> pd.DataFrame: + """Keep the top ``pool`` by momentum, order them by the signal (the rest + NaN). A name without links gets the neutral rank 0.5; momentum breaks ties.""" + rk = mom.rank(axis=1, ascending=False) + s = sig.rank(axis=1, pct=True).fillna(0.5) + 1e-3 * mom.rank(axis=1, pct=True) + return s.where(rk <= pool) + + +def fmt_sleeve(name: str, s: dict) -> str: + d, h = s["design"], s["holdout"] + return (f"{name:34s} alpha {s['all']['alpha']*100:5.1f}% t {s['all']['t']:5.2f} | 06-16 {d['alpha']*100:5.1f}% " + f"t {d['t']:5.2f} cagr {d['cagr']*100:5.1f} | 17-26 {h['alpha']*100:5.1f}% t {h['t']:5.2f} " + f"cagr {h['cagr']*100:5.1f} vs QQQ {h['bench']*100:5.1f} | mdd {s['mdd']*100:4.0f}%") + + +def part_signals(p, only: list[str] | None) -> dict: + from research.lab.alpha import bench_returns + from research.lab.us_sweep import _blend + pr = factor_ranks(p, "us", ["overnight_mom", "rec_126", "mom_12_1"]) + mom = _blend(pr["overnight_mom"], pr["rec_126"]).where(p.member) + c = p.close.ffill() + own21 = (c / c.shift(21) - 1).where(p.member) + own5 = (c / c.shift(5) - 1).where(p.member) + bench = bench_returns(p.dates) + cache = pickle.load(open(CACHE, "rb")) if os.path.exists(CACHE) else {} + builders = {"corr_peer": lambda: corr_peer_links(p), "lead_peer": lambda: lead_peer_links(p)} + if os.path.exists(os.path.join("data", "research", "sec_10k_links.csv")): + for role in ("customer", "supplier", "any"): + builders[f"10k_{role}"] = (lambda role=role: tenk_link_sets(p, role)) + res = {} + base, base_r = sleeve_stats(p, mom, 10, 10, bench) + print(fmt_sleeve("momentum top10 (reference)", base), flush=True) + res["momentum_top10"] = base + R = {"momentum_top10": base_r} + for link, build in builders.items(): + if only and link not in only: + continue + if link not in cache: + cache[link] = build() + pickle.dump(cache, open(CACHE, "wb")) + L = cache[link] + cover = np.mean([len(v[0]) if isinstance(v, tuple) else len(v) for v in L.values()]) + hz = 5 if link == "lead_peer" else 21 + every = 5 if link == "lead_peer" else 21 + own = own5 if hz == 5 else own21 + print(f"\n== {link}: {len(L)} months, {cover:.0f} stocks with links on average", flush=True) + for kind in ("mom", "gap"): + sig = link_signal(p, L, hz, kind) + fm = fama_macbeth(p, sig, {"own": own, "mom_12_1": pr["mom_12_1"], "overnight": pr["overnight_mom"]}, hz) + print(f" {link}_{kind}: FM b {fm['all']['b']:+.4f} t {fm['all']['t']:5.2f} | design t {fm['design']['t']:5.2f} " + f"| holdout t {fm['holdout']['t']:5.2f} | IC {fm['all']['ic']:+.4f} (t {fm['all']['ic_t']:.2f}) " + f"n {fm['all']['n_obs']:.0f}", flush=True) + sl, r = sleeve_stats(p, sig, 20, every, bench) + print(" " + fmt_sleeve(f"{link}_{kind} top20", sl), flush=True) + rs, r2 = sleeve_stats(p, resort(mom, sig), 10, 10, bench) + print(" " + fmt_sleeve(f"momentum top30 -> top10 by {kind}", rs), flush=True) + res[f"{link}_{kind}"] = {"fm": fm, "sleeve": sl, "resort": rs, + "corr_with_momentum": float(r.corr(base_r))} + R[f"{link}_{kind}"] = r + R[f"resort_{link}_{kind}"] = r2 + pd.DataFrame(R).to_pickle(os.path.join(RESULTS, "lead_lag_returns.pkl")) + return res + + +# --------------------------------------------------------------------------- +# earnings spillover +# --------------------------------------------------------------------------- + +def part_events(p) -> dict: + """For each earnings day of member j, the abnormal return of its + corr_peers (that did not report within 2 days) on the event day and over + days 1-5, 1-10 and 1-21 after it, regressed on j's event-day abnormal return.""" + ev = pd.read_csv(os.path.join("data", "research", "sec_8k_earnings.csv")) + ev["t"] = pd.to_datetime(ev["acceptanceDateTime"], utc=True).dt.tz_convert("America/New_York") + # after 16:00 the market reacts on the next session + ev["d"] = ev["t"].dt.normalize().dt.tz_localize(None) + pd.to_timedelta((ev["t"].dt.hour >= 16).astype(int), "D") + dates = p.dates + pos = dates.searchsorted(ev["d"]) + ev = ev[pos < len(dates)].assign(i=pos[pos < len(dates)]) + col = {t: k for k, t in enumerate(p.close.columns)} + ev = ev[ev["ticker"].isin(col)] + r = _ret(p) + mkt = r.where(p.member).mean(axis=1) + ab = r.sub(mkt, axis=0).to_numpy() + cache = pickle.load(open(CACHE, "rb")) if os.path.exists(CACHE) else {} + if "corr_peer" not in cache: + cache["corr_peer"] = corr_peer_links(p) + pickle.dump(cache, open(CACHE, "wb")) + L = cache["corr_peer"] + mes = np.array(sorted(L)) + rep = np.zeros(ab.shape, bool) + for t, i in zip(ev["ticker"], ev["i"]): + rep[max(i - 2, 0):i + 3, col[t]] = True + rows = [] + for t, i in zip(ev["ticker"], ev["i"]): + if i + 22 >= len(dates) or i < 300: + continue + me = mes[mes < dates[i]] + if not len(me): + continue + idx, nb = L[me[-1]] + j = col[t] + w = np.where(idx == j)[0] + if not len(w) or not p.member.iat[i, j]: + continue + a_j = ab[i, j] + if not np.isfinite(a_j): + continue + for q in nb[w[0]]: + if rep[i, q]: + continue + x = ab[i:i + 22, q] + if not np.isfinite(x[0]): + continue + x = np.nan_to_num(x) + rows.append((dates[i], t, p.close.columns[q], a_j, x[0], x[1:6].sum(), x[1:11].sum(), x[1:22].sum())) + E = pd.DataFrame(rows, columns=["date", "announcer", "peer", "ar_j", "ar0", "car1_5", "car1_10", "car1_21"]) + E.to_pickle(os.path.join(RESULTS, "lead_lag_events.pkl")) + out = {} + print(f"\n{len(E)} (event, peer) pairs, {E['announcer'].nunique()} announcers") + for name, sub in (("all", E), ("design", E[E["date"] <= DESIGN_END]), ("holdout", E[E["date"] > DESIGN_END])): + res = {} + for y in ("ar0", "car1_5", "car1_10", "car1_21"): + # cluster by event date: average the slope within each month, t over months + g = sub.groupby(sub["date"].dt.to_period("M")) + sl = g.apply(lambda d: np.polyfit(d["ar_j"], d[y], 1)[0] if len(d) > 20 and d["ar_j"].std() > 0 else np.nan).dropna() + res[y] = {"slope": float(sl.mean()), "t": float(sl.mean() / sl.std() * np.sqrt(len(sl)))} + # the trade: peers of big positive (negative) surprises that barely moved themselves + big = sub[sub["ar_j"].abs() > 0.05] + lag = big[(np.sign(big["ar0"]) != np.sign(big["ar_j"])) | (big["ar0"].abs() < 0.01)] + sgn = np.sign(lag["ar_j"]) + res["laggard_trade_car1_10"] = float((sgn * lag["car1_10"]).mean()) + res["laggard_trade_n"] = len(lag) + res["laggard_trade_t"] = float((sgn * lag["car1_10"]).groupby(lag["date"].dt.to_period("M")).mean().pipe( + lambda s: s.mean() / s.std() * np.sqrt(len(s)))) + out[name] = res + print(f" {name:8s} " + " | ".join(f"{y} slope {v['slope']:+.3f} t {v['t']:5.2f}" for y, v in res.items() if isinstance(v, dict)) + + f" | laggard trade 10d {res['laggard_trade_car1_10']*100:+.2f}% (t {res['laggard_trade_t']:.2f}, n {res['laggard_trade_n']})") + return out + + +# --------------------------------------------------------------------------- +# AI chain case study (hindsight grouping, descriptive only) +# --------------------------------------------------------------------------- + +def part_chain(p) -> dict: + from research.lab.review_live_us import AI_CHAIN + r = _ret(p) + out = {} + stages = {s: [t for t in ts if t in r.columns] for s, ts in AI_CHAIN.items()} + W = pd.DataFrame({s: r[ts].mean(axis=1) for s, ts in stages.items()}) + mkt = r.where(p.member).mean(axis=1) + W = W.sub(mkt, axis=0) # stage return over the EW market + wk = (1 + W).resample("W-FRI").prod() - 1 + for span in (("2016-01-01", "2022-12-31"), ("2023-01-01", "2026-09-24")): + x = wk.loc[span[0]:span[1]] + lead = pd.DataFrame(index=x.columns, columns=x.columns, dtype=float) + for a in x.columns: + for b in x.columns: + lead.loc[a, b] = x[a].shift(1).corr(x[b]) + n = len(x) + print(f"\nstage excess returns, weekly, {span[0][:4]}-{span[1][:4]} (n {n}, |corr| > {2/np.sqrt(n):.2f} is ~2 s.e.):" + "\nrow = stage in week w-1, column = stage in week w") + print(lead.round(2).to_string()) + out[f"{span[0][:4]}_{span[1][:4]}"] = lead.round(3).to_dict() + return out + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--part", default="signals", choices=["signals", "events", "chain", "all"]) + ap.add_argument("--only", nargs="*") + args = ap.parse_args() + p, _ = load("us") + out_path = os.path.join(RESULTS, "lead_lag.json") + out = json.load(open(out_path)) if os.path.exists(out_path) else {} + if args.part in ("signals", "all"): + out.setdefault("signals", {}).update(part_signals(p, args.only)) + if args.part in ("events", "all"): + out["events"] = part_events(p) + if args.part in ("chain", "all"): + out["chain"] = part_chain(p) + json.dump(out, open(out_path, "w"), indent=1, default=float) + + +if __name__ == "__main__": + main() diff --git a/research/lab/lead_lag_small.py b/research/lab/lead_lag_small.py new file mode 100644 index 0000000..94c8abe --- /dev/null +++ b/research/lab/lead_lag_small.py @@ -0,0 +1,360 @@ +""" +Lead-lag and supply-chain signals on S&P 400 + 600 stocks, where the +literature places these effects (low attention, large customers leading small +suppliers). Companion to research/lab/lead_lag.py (S&P 500). + + PYTHONPATH=. uv run python -m research.lab.fetch_supply_links --universe mid # 10-K links, ~1 h + PYTHONPATH=. uv run --with numba python -m research.lab.lead_lag_small + +Universe: PIT S&P 400 + 600 members (research/lab/panel_mid.py, iShares +holdings; covers 66% of members in 2019 and 83% in 2023; the missing names are +mostly delisted ones, so long-only returns carry survivorship bias) as the +stocks traded, and the S&P 1500 (those plus PIT S&P 500 members) as the pool +of linked names. Main sample 2019-01..2026-09, halves 2019-22 / 2023-26, +2017-26 for reference. Every signal is fixed before the results are seen and +tested at both horizons h = 21 and 63 trading days (signal window = holding +period = h), because the AI-chain stages lagged each other by 2-6 months. + +Signals for small stock i, known at the close of day t: + peer_mom / peer_gap the 10 S&P 1500 members whose daily residual returns + (S&P 1500 equal-weight beta removed) correlated most + with i over the last 252 days, re-estimated each + month end: their mean log return over the last h days + (mom), and that minus i's own (gap). + bigpeer_mom the same with the 5 most correlated S&P 500 members + only (large firms leading small ones, Hou 2007). + cust_mom / cust_gap S&P 1500 members named as customers in i's 10-Ks + (data/research/mid_10k_links.csv plus S&P 500-run + links whose filer is a mid/small member), valid for + 15 months from the filing day. +Event test: after an S&P 500 member's earnings day (8-K item 2.02), the +abnormal returns of the small stocks that have it among their 5 big peers, on +the day and over days 1-21 and 1-63, regressed on its day-0 abnormal return. + +Measures per signal and horizon (samples every 21 days; Newey-West t with +h/21 - 1 lags): + * Fama-MacBeth slope of the next-h open-to-open return rank on the signal + rank, controlling for the stock's 1-month return and 12-1 momentum ranks. + * top-minus-bottom quintile return, and the top quintile minus the + equal-weight small universe after 20 bp a side on its measured turnover. + * a top-30 sleeve (engine.simulate, next-open fills, 20 bp, buffer 2, + rebalanced every h days, phase-averaged): alpha against the IJH/IJR blend + and against QQQ. +The S&P 500 members get the same peer_mom / peer_gap / cust (S&P 500 10-K +links) FM rows over the same years for comparison. +""" + +from __future__ import annotations + +import json +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factors import _rolling_beta + +warnings.filterwarnings("ignore") +DIR = os.path.join("data", "research") +RESULTS = os.path.join(DIR, "results") +END = "2026-09-24" +SPANS = {"2019-26": ("2019-01-01", END), "2019-22": ("2019-01-01", "2022-12-31"), + "2023-26": ("2023-01-01", END), "2017-26": ("2017-01-01", END)} +HS = (21, 63) + + +# --------------------------------------------------------------------------- +# S&P 1500 frame +# --------------------------------------------------------------------------- + +def frame(): + from research.lab.factor_scan import load + from research.lab.panel_mid import load_us_mid + pb, _ = load("us") + ps = load_us_mid("midsmall") + dates = pb.dates.intersection(ps.dates) + dates = dates[dates <= pd.Timestamp(END)] + cols = sorted(set(pb.close.columns) | set(ps.close.columns)) + + def comb(a, b): + return a.reindex(index=dates, columns=cols).combine_first(b.reindex(index=dates, columns=cols)) + close = comb(pb.close, ps.close) + opn = comb(pb.open, ps.open) + mb = pb.member.reindex(index=dates, columns=cols).fillna(False).astype(bool) + ms = ps.member.reindex(index=dates, columns=cols).fillna(False).astype(bool) & ~mb + return {"pb": pb, "ps": ps, "dates": dates, "cols": cols, "close": close, "open": opn, + "mb": mb, "ms": ms, "mall": mb | ms} + + +def residuals(F) -> pd.DataFrame: + r = F["close"].pct_change(fill_method=None) + m = r.where(F["mall"]).mean(axis=1).fillna(0.0) + beta = _rolling_beta(r, m, 252, 126).shift(1) + F["ret"], F["mkt"] = r, m + return r - beta.mul(m, axis=0) + + +def month_ends(dates) -> list[pd.Timestamp]: + s = pd.Series(dates, index=dates) + return list(s.groupby(dates.to_period("M")).max()) + + +def peer_links(F, e: pd.DataFrame, focal: pd.DataFrame, pool: pd.DataFrame, k: int, win: int = 252, + start: str = "2015-06-01") -> dict: + """{month_end: {focal col index: peer col indices}} by residual correlation.""" + X = e.where(F["mall"]).to_numpy() + out = {} + fo, po = focal.to_numpy(), pool.to_numpy() + for me in month_ends(F["dates"]): + if me < pd.Timestamp(start): + continue + i1 = F["dates"].get_loc(me) + 1 + blk = X[i1 - win:i1] + ok = np.isfinite(blk).sum(axis=0) >= int(win * 0.8) + fi = np.where(fo[i1 - 1] & ok)[0] + pi = np.where(po[i1 - 1] & ok)[0] + if not len(fi) or not len(pi): + continue + Z = blk[:, np.union1d(fi, pi)] + allidx = np.union1d(fi, pi) + Z = np.nan_to_num((Z - np.nanmean(Z, 0)) / np.nanstd(Z, 0)) + pos = {c: j for j, c in enumerate(allidx)} + A = Z[:, [pos[c] for c in fi]] + B = Z[:, [pos[c] for c in pi]] + C = A.T @ B / win + C[fi[:, None] == pi[None, :]] = -np.inf + top = np.argsort(-C, axis=1)[:, :k] + out[me] = {int(f): pi[t] for f, t in zip(fi, top)} + return out + + +def customer_links(F, focal: pd.DataFrame, files: list[str], valid_days: int = 456) -> dict: + col = {t: i for i, t in enumerate(F["cols"])} + d = [] + for f in files: + if os.path.exists(f): + x = pd.read_csv(f, parse_dates=["filed"]) + d.append(x[x["n_customer"] > 0][["filer", "target", "filed"]]) + d = pd.concat(d).drop_duplicates() + d = d[d["filer"].isin(col) & d["target"].isin(col) & (d["filer"] != d["target"])] + out = {} + fo, ma = focal, F["mall"] + for me in month_ends(F["dates"]): + live = d[(d["filed"] <= me) & (d["filed"] > me - pd.Timedelta(days=valid_days))] + g = {} + for f, t in zip(live["filer"], live["target"]): + if fo.at[me, f] and ma.at[me, t]: + g.setdefault(col[f], set()).add(col[t]) + if g: + out[me] = {i: np.array(sorted(s)) for i, s in g.items()} + return out + + +def link_signal(F, links: dict, h: int, kind: str) -> pd.DataFrame: + lr = np.log1p(F["ret"].fillna(0.0)) + R = lr.rolling(h, min_periods=h).sum().to_numpy() + out = np.full(R.shape, np.nan) + mes = sorted(links) + for a, b in zip(mes, mes[1:] + [F["dates"][-1] + pd.Timedelta(days=1)]): + rows = np.where((F["dates"] > a) & (F["dates"] <= b))[0] + for i, nb in links[a].items(): + out[rows, i] = np.nanmean(R[np.ix_(rows, nb)], axis=1) + S = pd.DataFrame(out, index=F["dates"], columns=F["cols"]) + if kind == "gap": + S = S - pd.DataFrame(R, index=F["dates"], columns=F["cols"]) + return S + + +# --------------------------------------------------------------------------- +# measures +# --------------------------------------------------------------------------- + +def _nw_t(x: pd.Series, lags: int) -> float: + x = x.dropna() + n = len(x) + if n < 6: + return np.nan + u = x - x.mean() + s = (u ** 2).sum() / n + for l in range(1, lags + 1): + s += 2 * (1 - l / (lags + 1)) * (u[l:].to_numpy() * u[:-l].to_numpy()).sum() / n + return float(x.mean() / np.sqrt(s / n)) + + +def forward(F, h: int) -> pd.DataFrame: + o = F["open"].where(F["open"].notna(), F["close"]).ffill(limit=3) + return o.shift(-(h + 1)) / o.shift(-1) - 1 + + +def cross_section_tests(F, sig: pd.DataFrame, focal: pd.DataFrame, h: int, fwd: pd.DataFrame, + own21: pd.DataFrame, mom: pd.DataFrame, cost: float = 0.002) -> dict: + """FM slope, quintile spread and top-quintile excess, sampled every 21 days.""" + rows = [] + idx = np.arange(260, len(F["dates"]) - h - 1, 21) + for i in idx: + d = F["dates"][i] + f = focal.iloc[i] + df = pd.DataFrame({"y": fwd.iloc[i], "s": sig.iloc[i], "own": own21.iloc[i], "m": mom.iloc[i]})[f.values] + uni = df["y"].dropna() + df = df.dropna() + if len(df) < 60: + continue + rk = df.rank(pct=True) - 0.5 + X = np.column_stack([np.ones(len(rk)), rk[["s", "own", "m"]].to_numpy()]) + b, *_ = np.linalg.lstsq(X, rk["y"].to_numpy(), rcond=None) + q = pd.qcut(df["s"].rank(method="first"), 5, labels=False) + top, bot = df["y"][q == 4], df["y"][q == 0] + rows.append({"date": d, "b": b[1], "ls": top.mean() - bot.mean(), + "top_gross": top.mean() - uni.mean(), "names": frozenset(top.index), "n": len(df)}) + T = pd.DataFrame(rows).set_index("date") + out = {} + lags = max(h // 21 - 1, 0) + per_year = 252 / h + for name, (a, z) in SPANS.items(): + x = T.loc[a:z] + # non-overlapping portfolio periods for the h-day returns: every (h/21)-th sample + step = max(h // 21, 1) + xp = x.iloc[::step].copy() + prev = [None] + list(xp["names"].iloc[:-1]) + turn = [1.0 if p_ is None else len(n_ - p_) / max(len(n_), 1) for n_, p_ in zip(xp["names"], prev)] + xp["top_ex"] = xp["top_gross"] - 2 * cost * np.array(turn) + out[name] = {"fm_b": float(x["b"].mean()), "fm_t": _nw_t(x["b"], lags), + "ls_ann": float(xp["ls"].mean() * per_year), "ls_t": _nw_t(xp["ls"], 0), + "top_ex_ann": float(xp["top_ex"].mean() * per_year), "top_ex_t": _nw_t(xp["top_ex"], 0), + "n": float(x["n"].mean()), "periods": len(x)} + return out + + +def sleeve(F, sig: pd.DataFrame, h: int) -> dict: + from research.lab.alpha import alpha_stats, bench_returns + from research.lab.us_sweep import run + ps = F["ps"] + score = sig.reindex(index=ps.dates, columns=ps.close.columns).where(ps.member) + r = run(ps, score, [("2017-01-03", END)], top_n=30, buffer=2.0, every=h)["_rets"] + bq = bench_returns(r.index) + bm = ps.benchmark.reindex(r.index).pct_change(fill_method=None) if ps.benchmark is not None else None + out = {} + for name, (a, z) in SPANS.items(): + x = r.loc[a:z] + aq = alpha_stats(x, bq) + am = alpha_stats(x, bm) if bm is not None else {"alpha": np.nan, "t": np.nan} + out[name] = {"alpha_mid": am["alpha"], "t_mid": am["t"], "alpha_qqq": aq["alpha"], "t_qqq": aq["t"], + "cagr": aq["cagr"], "bench_mid_cagr": am.get("bench_cagr", np.nan)} + return out, r + + +def events(F, big_links: dict, e_abn: np.ndarray) -> dict: + ev = pd.read_csv(os.path.join(DIR, "sec_8k_earnings.csv")) + ev["t"] = pd.to_datetime(ev["acceptanceDateTime"], utc=True).dt.tz_convert("America/New_York") + ev["d"] = ev["t"].dt.normalize().dt.tz_localize(None) + pd.to_timedelta((ev["t"].dt.hour >= 16).astype(int), "D") + dates = F["dates"] + col = {t: i for i, t in enumerate(F["cols"])} + pos = dates.searchsorted(ev["d"]) + ev = ev[(pos < len(dates))].assign(i=pos[pos < len(dates)]) + ev = ev[ev["ticker"].isin(col) & (ev["d"] >= pd.Timestamp("2016-06-01"))] + mes = np.array(sorted(big_links)) + # invert: big j -> small followers, per month end + inv = {me: {} for me in mes} + for me in mes: + for i, nb in big_links[me].items(): + for j in nb: + inv[me].setdefault(int(j), []).append(i) + mb = F["mb"].to_numpy() + rows = [] + for t, i in zip(ev["ticker"], ev["i"]): + if i + 64 >= len(dates): + continue + j = col[t] + if not mb[i, j] or not np.isfinite(e_abn[i, j]): + continue + me = mes[mes < dates[i]] + if not len(me): + continue + for q in inv[me[-1]].get(j, []): + x = e_abn[i:i + 64, q] + if not np.isfinite(x[0]): + continue + x = np.nan_to_num(x) + rows.append((dates[i], t, F["cols"][q], e_abn[i, j], x[0], x[1:22].sum(), x[1:64].sum())) + E = pd.DataFrame(rows, columns=["date", "announcer", "small", "ar_j", "ar0", "car1_21", "car1_63"]) + E.to_pickle(os.path.join(RESULTS, "lead_lag_small_events.pkl")) + out = {} + for name, (a, z) in SPANS.items(): + sub = E[(E["date"] >= a) & (E["date"] <= z)] + res = {"pairs": len(sub)} + for y in ("ar0", "car1_21", "car1_63"): + g = sub.groupby(sub["date"].dt.to_period("M")) + sl = g.apply(lambda d: np.polyfit(d["ar_j"], d[y], 1)[0] if len(d) > 20 and d["ar_j"].std() > 0 else np.nan).dropna() + lags = 0 if y != "car1_63" else 2 + res[y] = {"slope": float(sl.mean()), "t": _nw_t(sl, lags)} + out[name] = res + return out + + +def main() -> None: + os.makedirs(RESULTS, exist_ok=True) + F = frame() + e = residuals(F) + c = F["close"].ffill() + own21 = c / c.shift(21) - 1 + mom = c.shift(21) / c.shift(252) - 1 + fwd = {h: forward(F, h) for h in HS} + ms, mb, mall = F["ms"], F["mb"], F["mall"] + cov = pd.DataFrame({"small_members": ms.sum(axis=1), "big_members": mb.sum(axis=1)}).resample("YE").mean() + print("members per day (yearly mean):\n", cov.loc["2016":].round(0).to_string(), flush=True) + + L = {"peer": peer_links(F, e, ms, mall, 10), "bigpeer": peer_links(F, e, ms, mb, 5), + "peer_sp500": peer_links(F, e, mb, mall, 10)} + have_mid = os.path.exists(os.path.join(DIR, "mid_10k_links.csv")) + if have_mid: + L["cust"] = customer_links(F, ms, [os.path.join(DIR, "mid_10k_links.csv"), os.path.join(DIR, "sec_10k_links.csv")]) + L["cust_sp500"] = customer_links(F, mb, [os.path.join(DIR, "sec_10k_links.csv")]) + for k, v in L.items(): + n = np.mean([len(x) for x in list(v.values()) if True][-90:]) if v else 0 + print(f"links {k}: {len(v)} months, {n:.0f} focal stocks with links (last 90 months avg)", flush=True) + + variants = [] + for h in HS: + variants += [("peer_mom", "peer", "mom", h, "small"), ("peer_gap", "peer", "gap", h, "small"), + ("bigpeer_mom", "bigpeer", "mom", h, "small")] + if have_mid: + variants += [("cust_mom", "cust", "mom", h, "small"), ("cust_gap", "cust", "gap", h, "small")] + variants += [("sp500_peer_mom", "peer_sp500", "mom", h, "big"), ("sp500_peer_gap", "peer_sp500", "gap", h, "big"), + ("sp500_cust_mom", "cust_sp500", "mom", h, "big")] + res, rets, trials = {}, {}, [] + for name, lk, kind, h, uni in variants: + sig = link_signal(F, L[lk], h, kind) + focal = ms if uni == "small" else mb + cs = cross_section_tests(F, sig, focal, h, fwd[h], own21, mom) + row = {"cs": cs} + if uni == "small": + sl, r = sleeve(F, sig, h) + row["sleeve"] = sl + rets[f"{name}_{h}"] = r + res[f"{name}_{h}"] = row + m = cs["2019-26"] + line = (f"{name:16s} h{h:<3d} FM t 19-26 {m['fm_t']:5.2f} | 19-22 {cs['2019-22']['fm_t']:5.2f} | 23-26 " + f"{cs['2023-26']['fm_t']:5.2f} | 17-26 {cs['2017-26']['fm_t']:5.2f} | L-S {m['ls_ann']*100:5.1f}%/yr " + f"(t {m['ls_t']:.2f}) | top-EW net {m['top_ex_ann']*100:5.1f}%/yr (t {m['top_ex_t']:.2f}) n {m['n']:.0f}") + if "sleeve" in row: + s = row["sleeve"]["2019-26"] + line += f" | top30 a_mid {s['alpha_mid']*100:5.1f}% (t {s['t_mid']:.2f}) a_qqq {s['alpha_qqq']*100:5.1f}% (t {s['t_qqq']:.2f})" + print(line, flush=True) + trials.append({"name": name, "h": h, "universe": uni, **{f"fm_t_{k}": v["fm_t"] for k, v in cs.items()}, + "ls_ann_1926": m["ls_ann"], "top_ex_ann_1926": m["top_ex_ann"]}) + # event test: S&P 500 earnings -> small followers + abn = (F["ret"].sub(F["mkt"], axis=0)).to_numpy() + res["events"] = events(F, L["bigpeer"], abn) + for k, v in res["events"].items(): + print(f"events {k}: pairs {v['pairs']} | " + " | ".join(f"{y} slope {v[y]['slope']:+.3f} t {v[y]['t']:.2f}" + for y in ("ar0", "car1_21", "car1_63")), flush=True) + trials += [{"name": f"event_{y}", "h": 0, "universe": "small"} for y in ("ar0", "car1_21", "car1_63")] + json.dump(res, open(os.path.join(RESULTS, "lead_lag_small.json"), "w"), indent=1, default=float) + pd.DataFrame(rets).to_pickle(os.path.join(RESULTS, "lead_lag_small_returns.pkl")) + pd.DataFrame(trials).to_csv(os.path.join(RESULTS, "lead_lag_small_trials.csv"), index=False) + print(f"{len(trials)} trials", flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/review_live_us.py b/research/lab/review_live_us.py new file mode 100644 index 0000000..b91bb37 --- /dev/null +++ b/research/lab/review_live_us.py @@ -0,0 +1,303 @@ +""" +Review of the US replay books over 2026-04-08 .. 2026-09-28: what drove the +market, where each strategy's return came from, and how good its entries and +exits were. + + PYTHONPATH=. uv run python -m research.lab.review_live_us + +Inputs are the live data files (data/us*.csv, data/etfs.csv, the original +book's data/trader_us_*.json). The books are point-in-time replays, not live +trading (see project memory); every decision in them used only data up to +its date. Books trade at the close, so the weights recorded on day d earn +the close-to-close return of day d+1. + +AI_CHAIN groups S&P 500 members into stages of the AI build-out. The grouping +was written in 2026 with hindsight and is used here to describe the period, +never as a trading signal. +""" + +from __future__ import annotations + +import glob +import json +import os +import warnings + +import numpy as np +import pandas as pd + +from universe_history import membership_mask + +warnings.filterwarnings("ignore") +START, END = "2026-04-08", "2026-09-28" +RESULTS = os.path.join("data", "research", "results") + +AI_CHAIN = { + "hyperscaler": ["MSFT", "GOOGL", "AMZN", "META", "ORCL"], + "ai_software": ["PLTR", "NOW", "CRM", "ADBE", "APP", "INTU", "DDOG", "CRWD", "SNOW", "WDAY"], + "servers_network": ["DELL", "SMCI", "HPE", "ANET", "CSCO", "JBL", "FLEX"], + "compute_chips": ["NVDA", "AMD", "AVGO", "MRVL", "INTC", "QCOM"], + "memory_storage": ["MU", "SNDK", "WDC", "STX", "NTAP", "PSTG"], + "optical_components": ["LITE", "COHR", "CIEN", "GLW", "APH", "TEL"], + "semi_equipment": ["AMAT", "LRCX", "KLAC", "TER", "Q"], + "power_datacenter": ["VST", "CEG", "NRG", "GEV", "ETN", "VRT", "BE", "PWR", "FIX", "EME", "EQIX", "DLR"], + "analog_semis": ["TXN", "ADI", "MCHP", "ON", "NXPI", "MPWR"], +} +FOCUS = ["overnight_recovery_top10", "us_momentum_qqq_trend", "us_three_sleeve", "lowrisk_top20", + "fc_up_cap_mom_gap_weekly", "fc_rec_mfilt_deep_upvol_monthly", "recovery_mom_top10", + "momentum", "trend_rider_v7_vt24"] + + +def load_prices(): + c = pd.read_csv("data/us.csv", index_col=0, parse_dates=True) + e = pd.read_csv("data/etfs.csv", index_col=0, parse_dates=True) + c = c.join(e[[x for x in e.columns if x not in c.columns]], how="left") + c = c.sort_index() + etf = [x for x in e.columns if x in c.columns] + c[etf] = c[etf].ffill(limit=3) + return c + + +def load_books() -> dict[str, dict]: + out = {} + for f in sorted(glob.glob("data/trader_us_*.json")): + d = json.load(open(f)) + if d.get("daily_log"): + out[d["strategy"]] = d + return out + + +def book_weights(book: dict, cols) -> tuple[pd.DataFrame, pd.Series]: + """Weights (value / equity) after each day's close, and daily equity.""" + rows, eq = {}, {} + for r in book["daily_log"]: + e = r["equity"] + eq[r["date"]] = e + rows[r["date"]] = {t: h["value"] / e for t, h in r["holdings"].items()} + w = pd.DataFrame.from_dict(rows, orient="index").fillna(0.0) + w.index = pd.to_datetime(w.index) + eq = pd.Series(eq); eq.index = pd.to_datetime(eq.index) + w, eq = w.sort_index(), eq.sort_index() + return w[~w.index.duplicated(keep="last")], eq[~eq.index.duplicated(keep="last")] + + +def stage_of() -> dict[str, str]: + return {t: s for s, ts in AI_CHAIN.items() for t in ts} + + +def market_view(C: pd.DataFrame, members: list[str], sectors: pd.Series) -> dict: + r = C.loc[START:END, members].ffill() + tot = (r.iloc[-1] / r.iloc[0] - 1).dropna() + st = stage_of() + df = pd.DataFrame({"ret": tot, "sector": sectors.reindex(tot.index), + "stage": [st.get(t, "") for t in tot.index]}) + out = {"n": len(df), "ew": float(df["ret"].mean()), "median": float(df["ret"].median()), + "spy": float(C.loc[END, "SPY"] / C.loc[START, "SPY"] - 1), + "qqq": float(C.loc[END, "QQQ"] / C.loc[START, "QQQ"] - 1), + "smh": float(C.loc[END, "SMH"] / C.loc[START, "SMH"] - 1)} + out["by_sector"] = df.groupby("sector")["ret"].agg(["mean", "median", "count"]).sort_values("mean", ascending=False) + out["by_stage"] = df[df["stage"] != ""].groupby("stage")["ret"].agg(["mean", "median", "count"]).sort_values("mean", ascending=False) + top = df.sort_values("ret", ascending=False).head(40) + out["top40"] = top + out["top40_ai_share"] = float((top["stage"].isin([s for s in AI_CHAIN if s != "analog_semis"])).mean()) + return out + + +def attribution(w: pd.DataFrame, R: pd.DataFrame, bench_members: pd.DataFrame, groups: pd.Series) -> dict: + """Daily Brinson against the equal-weight PIT S&P 500: group allocation and + selection within groups, summed arithmetically over the period.""" + dates = w.index.intersection(R.index) + alloc = sel = act = 0.0 + per_group = {} + for d0, d1 in zip(dates[:-1], dates[1:]): + wp = w.loc[d0] + wp = wp[wp > 0] + r1 = R.loc[d1] + mem = bench_members.loc[d0] + bm = [t for t in mem.index[mem] if np.isfinite(r1.get(t, np.nan))] + g_b = groups.reindex(bm).fillna("Other") + rb = r1[bm] + Rb = rb.mean() + g_p = groups.reindex(wp.index).fillna("Other") + rp = r1.reindex(wp.index).fillna(0.0) + for g in set(g_b) | set(g_p): + wb_g = float((g_b == g).mean()) + rb_g = float(rb[g_b == g].mean()) if (g_b == g).any() else Rb + wp_g = float(wp[g_p == g].sum()) + rp_g = float((wp[g_p == g] * rp[g_p == g]).sum() / wp_g) if wp_g > 0 else rb_g + a = (wp_g - wb_g) * (rb_g - Rb) + s = wp_g * (rp_g - rb_g) + alloc += a; sel += s + pg = per_group.setdefault(g, [0.0, 0.0]) + pg[0] += a; pg[1] += s + cash = 1 - wp.sum() + act += float((wp * rp).sum()) - Rb # cash earns 0 in the books + pg = pd.DataFrame(per_group, index=["allocation", "selection"]).T + pg["total"] = pg.sum(axis=1) + return {"allocation": alloc, "selection": sel, "active_sum": act, "by_group": pg.sort_values("total")} + + +def trade_quality(book: dict, C: pd.DataFrame, ew: pd.Series, horizons=(10, 21)) -> pd.DataFrame: + """Forward return after each full exit / entry, minus the equal-weight index.""" + rows = [] + idx = C.index + for t in book["trade_log"]: + d = pd.Timestamp(t["date"]) + if t["ticker"] not in C.columns or d not in idx: + continue + i = idx.get_loc(d) + row = {"date": d, "ticker": t["ticker"], "action": t["action"]} + for h in horizons: + if i + h < len(idx): + p0, p1 = C[t["ticker"]].iloc[i], C[t["ticker"]].iloc[i + h] + b = ew.iloc[i + 1:i + h + 1] + row[f"fwd{h}"] = p1 / p0 - 1 - ((1 + b).prod() - 1) + rows.append(row) + return pd.DataFrame(rows) + + +def exits_entries(w: pd.DataFrame) -> tuple[list, list]: + held = w > 1e-6 + prev = held.shift(1, fill_value=False) + entries = [(d, t) for d, t in zip(*np.where(held & ~prev)) if d > 0] + exits = [(d, t) for d, t in zip(*np.where(~held & prev))] + return ([(w.index[d], w.columns[t]) for d, t in entries], + [(w.index[d], w.columns[t]) for d, t in exits]) + + +def winner_capture(w: pd.DataFrame, C: pd.DataFrame, winners: list[str]) -> pd.DataFrame: + """For each period winner: days held, first entry, and the log-return earned + while held as a share of its total log-return over the period.""" + R = np.log(C.loc[START:END].ffill()).diff() + rows = [] + for t in winners: + tot = float(R[t].sum()) + if t in w.columns: + h = (w[t] > 1e-6).reindex(R.index).fillna(False).shift(1, fill_value=False) + got = float(R[t][h].sum()) + first = w.index[w[t] > 1e-6][0] if (w[t] > 1e-6).any() else None + days = int(h.sum()) + else: + got, first, days = 0.0, None, 0 + rows.append({"ticker": t, "total_logret": tot, "captured_logret": got, + "capture": got / tot if tot else np.nan, "days_held": days, "first_entry": first}) + return pd.DataFrame(rows) + + +def stage_timing(C: pd.DataFrame) -> pd.DataFrame: + """Equal-weight stage indices from 2025-01; the date each one first reached + 25/50/75% of its log gain up to END.""" + rows = {} + px = C.loc["2025-01-01":END].ffill() + for s, ts in AI_CHAIN.items(): + ts = [t for t in ts if t in px.columns] + r = np.log(px[ts]).diff().mean(axis=1).fillna(0.0) + cum = r.cumsum() + tot = cum.iloc[-1] + rows[s] = {"logret_2025_01_to_end": float(tot)} + for q in (0.25, 0.5, 0.75): + hit = cum.index[cum >= q * tot] + rows[s][f"t{int(q*100)}"] = hit[0].date() if len(hit) and tot > 0 else None + return pd.DataFrame(rows).T + + +def main() -> None: + os.makedirs(RESULTS, exist_ok=True) + C = load_prices() + sec = pd.read_csv("data/us_sectors.csv").set_index("ticker") + from research.lab.factor_scan import load + p, _ = load("us") + sectors = p.extra["sector"].combine_first(sec["GICS Sector"]) + dates = C.loc[START:END].index + M = membership_mask(dates, tickers=[c for c in C.columns if c in sectors.index]) + members0 = list(M.columns[M.iloc[0]]) + R = C.pct_change(fill_method=None) + ew = R[members0].loc[dates].mean(axis=1) + + mv = market_view(C, members0, sectors) + print(f"period {START}..{END}: SPY {mv['spy']*100:.1f}% QQQ {mv['qqq']*100:.1f}% SMH {mv['smh']*100:.1f}% " + f"EW members {mv['ew']*100:.1f}% (median {mv['median']*100:.1f}%, n {mv['n']})") + print("\nby GICS sector (equal-weight member return):") + print((mv["by_sector"].assign(mean=lambda x: x["mean"] * 100, median=lambda x: x["median"] * 100)).round(1).to_string()) + print("\nby AI-chain stage (hindsight grouping):") + print((mv["by_stage"].assign(mean=lambda x: x["mean"] * 100, median=lambda x: x["median"] * 100)).round(1).to_string()) + print(f"\ntop 40 winners: {mv['top40_ai_share']*100:.0f}% are AI-chain names") + print((mv["top40"].assign(ret=lambda x: x["ret"] * 100)).round(1).to_string()) + print("\nstage timing (EW index from 2025-01; date first reaching 25/50/75% of the gain):") + timing = stage_timing(C) + print(timing.to_string()) + + books = load_books() + st = stage_of() + grp_stage = pd.Series({t: st.get(t, "not_ai") for t in C.columns}) + winners = list(mv["top40"].index) + summary, attrib, caps, tq_all = [], {}, {}, {} + for name, b in books.items(): + w, eq = book_weights(b, C.columns) + w = w.loc[START:END] + ret = float(eq.loc[:END].iloc[-1] / b["initial_capital"] - 1) + sd = str(eq.index[0].date()) + spy = float(C.loc[END, "SPY"] / C.loc[sd:, "SPY"].iloc[0] - 1) if sd else np.nan + stk = w[[c for c in w.columns if c in sectors.index]] + ai = stk[[c for c in stk.columns if st.get(c, "") and st.get(c) != "analog_semis"]].sum(axis=1) + tq = trade_quality(b, C, R[members0].mean(axis=1)) + row = {"strategy": name, "start": sd, "ret": ret, "spy_same_span": spy, + "stock_w": float(stk.sum(axis=1).mean()), "ai_w": float(ai.mean()), + "n_trades": len(b["trade_log"])} + if len(tq): + s = tq[tq["action"] == "SELL"]; bu = tq[tq["action"] == "BUY"] + row.update({"sell_fwd21": s["fwd21"].mean() if "fwd21" in s else np.nan, + "buy_fwd21": bu["fwd21"].mean() if "fwd21" in bu else np.nan}) + cap = winner_capture(w, C, winners) + row["winner_capture"] = float(cap["captured_logret"].sum() / cap["total_logret"].sum()) + row["winners_held"] = int((cap["days_held"] > 0).sum()) + summary.append(row) + if name in FOCUS and sd <= START: + attrib[name] = {"sector": attribution(w, R, M, sectors), "ai_stage": attribution(w, R, M, grp_stage)} + caps[name] = cap + tq_all[name] = tq + S = pd.DataFrame(summary).set_index("strategy").sort_values("ret", ascending=False) + print("\nall books (ret since own start; ai_w = mean weight in AI-chain names; winner_capture = share of the top-40 " + "winners' log gain earned; sell/buy_fwd21 = 21-day excess return after the trade):") + print(S.assign(**{c: S[c] * 100 for c in ["ret", "spy_same_span", "stock_w", "ai_w", "sell_fwd21", "buy_fwd21", "winner_capture"]}) + .round(1).to_string()) + full = S[S["start"] <= START] + print(f"\nbooks from {START}: corr(ret, ai_w) {full['ret'].corr(full['ai_w']):.2f}, " + f"corr(ret, winner_capture) {full['ret'].corr(full['winner_capture']):.2f}") + + for name, a in attrib.items(): + print(f"\n== {name}: Brinson vs EW S&P 500 (sum of daily terms, %)") + for k in ("sector", "ai_stage"): + x = a[k] + print(f" by {k}: allocation {x['allocation']*100:6.1f} selection {x['selection']*100:6.1f} " + f"active {x['active_sum']*100:6.1f}") + print(" ", ", ".join(f"{g} {v*100:+.1f}" for g, v in x["by_group"]["total"].items() if abs(v) > 0.01)) + cap = caps[name] + held = cap[cap["days_held"] > 0].sort_values("first_entry") + print(" top-40 winners held (first entry, capture):", + ", ".join(f"{r.ticker} {str(r.first_entry.date())[5:]} {r.capture*100:.0f}%" for r in held.itertuples())) + miss = cap[cap["days_held"] == 0].sort_values("total_logret", ascending=False) + print(" never held:", ", ".join(miss["ticker"].head(20))) + tq = tq_all[name] + if len(tq): + s = tq[tq["action"] == "SELL"] + print(f" exits: {len(s)}, 21d excess after exit mean {s['fwd21'].mean()*100:+.1f}% " + f"median {s['fwd21'].median()*100:+.1f}%, share >+10%: {(s['fwd21'] > .10).mean()*100:.0f}%") + + out = {"market": {k: v for k, v in mv.items() if not isinstance(v, pd.DataFrame)}, + "summary": S.reset_index().to_dict(orient="records"), + "by_sector": mv["by_sector"].reset_index().to_dict(orient="records"), + "by_stage": mv["by_stage"].reset_index().to_dict(orient="records"), + "top40": mv["top40"].reset_index().rename(columns={"index": "ticker"}).to_dict(orient="records"), + "timing": timing.reset_index().astype(str).to_dict(orient="records"), + "attribution": {n: {k: {"allocation": v["allocation"], "selection": v["selection"], + "active": v["active_sum"], "by_group": v["by_group"]["total"].to_dict()} + for k, v in a.items()} for n, a in attrib.items()}, + "capture": {n: c.astype({"first_entry": str}).to_dict(orient="records") for n, c in caps.items()}} + with open(os.path.join(RESULTS, "review_live_us.json"), "w") as fh: + json.dump(out, fh, indent=1, default=str) + pd.to_pickle({"trades": tq_all}, os.path.join(RESULTS, "review_live_us_trades.pkl")) + + +if __name__ == "__main__": + main() diff --git a/research/lab/sector_cycle.py b/research/lab/sector_cycle.py new file mode 100644 index 0000000..20ba863 --- /dev/null +++ b/research/lab/sector_cycle.py @@ -0,0 +1,781 @@ +""" +Can sector cyclicality be modelled well enough to improve the US book? +Business-cycle sector rotation, sector momentum and reversal, cycle-position +signals, and cycle-conditioned versions of the momentum sleeve. + + PYTHONPATH=. uv run python -m research.lab.sector_cycle fetch # FRED series + sector ETFs + PYTHONPATH=. uv run python -m research.lab.sector_cycle # all tests + +Yardstick (as in beat_qqq.py): alpha against QQQ and SPY (Newey-West t), +design 2006-2016 vs holdout 2017-2026. ETF-level tests also run from 1999 +and are split at 2008. Every variant is fixed before it is run and logged in +data/research/results/sector_cycle_trials.csv. + +Inputs and their lags + sector ETFs the nine 1998 SPDRs (XLB XLE XLF XLI XLK XLP XLU XLV XLY), + Yahoo adjusted closes. XLC (2018) and XLRE (2015) are too short. + GICS sector p.extra["sector"]: today's classification, not point in time. + SIC 2-digit data/research/sec_sic.csv: today's code, not point in time. + FRED a monthly value for month m (CFNAI, INDPRO, UNRATE, CPIAUCSL) + is used from the last trading day of m+1 (each is released + before then); weekly ICSA 7 days after the week's date; daily + T10Y3M and BAA10Y the next business day. Latest vintage only: + revisions are not modelled (ALFRED vintages were not + reachable), which flatters CFNAI and INDPRO. ISM PMI (NAPM) is + no longer on FRED. + +Tests + (a) Business-cycle rotation. Phase by the OECD four-phase rule + (below/above trend x rising/falling) on CFNAI-MA3 (trend = 0, + direction = 3-month change) and on INDPRO 12-month log growth (trend = + its trailing 10-year mean). Phase -> favoured sectors from Fidelity's + business-cycle approach, not fitted. + (b) Sector momentum: cross-sectional 12-1 top 3, time-series 12-month (a + sector is held only while it beat T-bills over 12 months), long-horizon + reversal (bottom 3 by months 60..13). ETFs and EW PIT GICS sectors. + (c) Cycle position: price position in the 3-year range, rank IC against + the next 3/6/12 months across ETFs, GICS sectors and SIC-2 industries; + EDGAR revenue-growth level and acceleration by sector (2011+); the + semiconductor cycle (^SOX from 1994) with position, relative momentum + and SIC 3674 revenue acceleration. + (d) The momentum sleeve (overnight + rec_126, top 10, every 10 days, + buffer 2 = overnight_recovery_top10) by macro state; a cash switch in + the state picked on 2006-2016 only; its picks restricted to the CFNAI + phase's favoured sectors, or with the phase's unfavoured sectors removed. +""" + +from __future__ import annotations + +import functools +import json +import os +import sys +import warnings + +import numpy as np +import pandas as pd + +from research.lab.alpha import alpha_stats, bench_returns +from research.lab.factor_scan import RESULTS +from research.lab.levered_book import tbill +from research.lab.report import stats +from research.lab.us_final import monthly_blend + +warnings.filterwarnings("ignore") + +EXT = os.path.join("data", "research", "ext") +YAHOO = os.path.join(EXT, "yahoo") +SIC = os.path.join("data", "research", "sec_sic.csv") +ETF9 = ["XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY"] +FRED_SERIES = ["CFNAI", "INDPRO", "UNRATE", "CPIAUCSL", "ICSA", "T10Y3M", "BAA10Y", "USREC"] +MONTHLY = {"CFNAI", "INDPRO", "UNRATE", "CPIAUCSL", "USREC"} +LONG = ("2006-09-27", "2026-09-24") +SPANS = {"full": (None, None), "design": ("2006-01-01", "2016-12-31"), "holdout": ("2017-01-01", None), + "pre08": (None, "2007-12-31"), "post08": ("2008-01-01", None)} +PHASES = ["early", "mid", "late", "recession"] + +# Fidelity, "The Business Cycle Approach to Equity Sector Investing": sectors +# that have tended to beat / trail the market in each phase. +FAVOURED = { + "early": ["Financials", "Real Estate", "Consumer Discretionary", "Industrials", + "Information Technology", "Materials"], + "mid": ["Information Technology", "Communication Services", "Industrials"], + "late": ["Energy", "Materials", "Consumer Staples", "Health Care", "Utilities"], + "recession": ["Consumer Staples", "Utilities", "Health Care"], +} +UNFAVOURED = { + "early": ["Consumer Staples", "Energy", "Health Care", "Communication Services", "Utilities"], + "mid": ["Materials", "Utilities"], + "late": ["Information Technology", "Consumer Discretionary"], + "recession": ["Industrials", "Information Technology", "Real Estate"], +} +ETF_SECTOR = {"XLB": "Materials", "XLE": "Energy", "XLF": "Financials", "XLI": "Industrials", + "XLK": "Information Technology", "XLP": "Consumer Staples", "XLU": "Utilities", + "XLV": "Health Care", "XLY": "Consumer Discretionary"} + +TRIALS: list[dict] = [] + + +# --------------------------------------------------------------------------- +# data +# --------------------------------------------------------------------------- + +def fetch() -> None: + from research.lab.fetch_long_history import fetch_fred, fetch_yahoo + df = fetch_fred(FRED_SERIES) + for s in FRED_SERIES: + df[[s]].dropna().rename_axis("observation_date").to_csv(os.path.join(EXT, f"fred_{s}.csv")) + print(s, df[s].dropna().index[0].date(), df[s].dropna().index[-1].date()) + fetch_yahoo(ETF9 + ["XLC", "XLRE", "^SOX", "SMH", "SOXX"], start="1990-01-01") + + +def fred(s: str) -> pd.Series: + x = pd.read_csv(os.path.join(EXT, f"fred_{s}.csv"), index_col=0, parse_dates=True).iloc[:, 0] + return pd.to_numeric(x, errors="coerce").dropna() + + +def available(x: pd.Series, s: str) -> pd.Series: + """Re-index a FRED series (or a transform of it) by the date it is known.""" + if s in MONTHLY: + idx = x.index + pd.offsets.MonthEnd(2) # obs dated m-01 -> end of m+1 + elif s == "ICSA": + idx = x.index + pd.Timedelta(days=7) + else: + idx = x.index + pd.offsets.BDay(1) + return pd.Series(x.values, index=idx).groupby(level=0).last() + + +def asof(x: pd.Series, dates: pd.DatetimeIndex) -> pd.Series: + return x.reindex(dates.union(x.index)).ffill().reindex(dates) + + +@functools.lru_cache +def yahoo(t: str) -> pd.Series: + path = os.path.join(YAHOO, t.replace("^", "_").replace("=", "_") + ".csv") + return pd.read_csv(path, index_col=0, parse_dates=True)["Close"] + + +def etf_prices(tickers=ETF9) -> pd.DataFrame: + return pd.DataFrame({t: yahoo(t) for t in tickers}).sort_index() + + +def month_ends(dates: pd.DatetimeIndex) -> pd.DatetimeIndex: + s = pd.Series(dates, index=dates) + return pd.DatetimeIndex(s.groupby(dates.to_period("M")).max().values) + + +# --------------------------------------------------------------------------- +# macro phases and states (index = date the value is known) +# --------------------------------------------------------------------------- + +def _four_phase(level_above: pd.Series, rising: pd.Series) -> pd.Series: + ph = np.select([~level_above & rising, level_above & rising, level_above & ~rising, ~level_above & ~rising], + PHASES, default=None) + return pd.Series(ph, index=level_above.index) + + +def phase_cfnai(raw: bool = False) -> pd.Series: + """raw=True: indexed by the observation month (look-ahead; for the oracle + bound and the descriptive table only).""" + m = fred("CFNAI").rolling(3).mean().dropna() + d = m.diff(3) + ok = d.notna() + ph = _four_phase(m[ok] > 0, d[ok] > 0) + return ph if raw else available(ph, "CFNAI") + + +def phase_indpro() -> pd.Series: + g = np.log(fred("INDPRO")).diff(12) + trend = g.rolling(120, min_periods=60).mean() + d = g.diff(3) + ok = g.notna() & trend.notna() & d.notna() + return available(_four_phase((g > trend)[ok], (d > 0)[ok]), "INDPRO") + + +def states(dates: pd.DatetimeIndex) -> pd.DataFrame: + """Macro states known at each date (strings / bools).""" + out = pd.DataFrame(index=dates) + out["cfnai"] = asof(phase_cfnai(), dates) + out["indpro"] = asof(phase_indpro(), dates) + out["curve_inverted"] = asof(available(fred("T10Y3M") < 0, "T10Y3M"), dates) + baa = fred("BAA10Y") + out["credit_widening"] = asof(available(baa - baa.shift(63) > 0, "BAA10Y").astype(float), dates) > 0 + ic = fred("ICSA").rolling(13).mean() + out["claims_rising"] = asof(available(ic > ic.shift(26), "ICSA").astype(float), dates) > 0 + u = fred("UNRATE").rolling(3).mean() + out["sahm"] = asof(available(u - u.rolling(12).min().shift(1) >= 0.5, "UNRATE").astype(float), dates) > 0 + cpi = np.log(fred("CPIAUCSL")).diff(12) + out["inflation_rising"] = asof(available(cpi.diff(3) > 0, "CPIAUCSL").astype(float), dates) > 0 + return out + + +# --------------------------------------------------------------------------- +# books +# --------------------------------------------------------------------------- + +def etf_book(W: pd.DataFrame, px: pd.DataFrame, cost: float = 0.001) -> pd.Series: + """Monthly ETF book. W holds target weights on signal dates (rows sum to at + most 1, the rest earns T-bills); a signal is traded at the next day's close, + and holdings drift between rebalances.""" + px = px[W.columns].ffill() + dates = px.index + ret = px.pct_change(fill_method=None).fillna(0.0).to_numpy() + rf = (tbill(dates).fillna(0.0) / 252).to_numpy() + T = W.reindex(dates).to_numpy() + valid = ~np.isnan(T).all(axis=1) + first = int(np.argmax(valid)) + 1 + h = np.zeros(px.shape[1]); cash = 1.0 + eq = np.full(len(dates), np.nan); eq[first - 1] = 1.0 + for i in range(first, len(dates)): + h = h * (1 + ret[i]); cash *= 1 + rf[i] + V = h.sum() + cash + if valid[i - 1]: + tgt = np.nan_to_num(T[i - 1]) * V + c = np.abs(tgt - h).sum() * cost + V -= c + h = np.nan_to_num(T[i - 1]) * V + cash = V - h.sum() + eq[i] = V + return pd.Series(eq, index=dates).pct_change().dropna() + + +def stock_book(p, W: pd.DataFrame) -> pd.Series: + """Monthly book of stocks from explicit target weights (engine.simulate, + next-open fills, 10 bp, PIT members).""" + from research.lab.engine import simulate + sp = p.slice("2005-01-01", LONG[1]) + r = simulate(sp, target_weights=W.reindex(columns=sp.close.columns), every="M", exec_at="open") + first = W.dropna(how="all").index[0] + return r.returns.loc[first:].iloc[1:] + + +def sector_weights(chosen: dict, members: pd.DataFrame, sector: pd.Series) -> pd.DataFrame: + """chosen: {signal date: {sector: weight}} -> equal weight within each sector.""" + rows = {} + for d, sw in chosen.items(): + m = members.loc[d] + w = pd.Series(0.0, index=members.columns) + for s, x in sw.items(): + cols = m.index[m & (sector == s)] + if len(cols) and x > 0: + w[cols] = x / len(cols) + rows[d] = w + return pd.DataFrame(rows).T + + +# --------------------------------------------------------------------------- +# evaluation +# --------------------------------------------------------------------------- + +@functools.lru_cache +def momentum_sleeve() -> pd.Series: + return pd.read_pickle(os.path.join(RESULTS, "beat_qqq_sleeves.pkl"))["momentum"] + + +def _mdd(r: pd.Series) -> float: + e = (1 + r).cumprod() + return float((e / e.cummax() - 1).min()) + + +def evaluate(name: str, r: pd.Series, test: str, ew: pd.Series | None = None, log: bool = True) -> dict: + r = r.dropna() + out = {"name": name, "test": test, "start": str(r.index[0].date())} + for k, (a, b) in SPANS.items(): + x = r.loc[a:b] + if len(x) < 250: + continue + q = alpha_stats(x, bench_returns(x.index, "QQQ")) + s = alpha_stats(x, bench_returns(x.index, "SPY")) + st = stats((1 + x).cumprod()) + row = {"alpha_qqq": q["alpha"], "t_qqq": q["t"], "beta_qqq": q["beta"], + "alpha_spy": s["alpha"], "t_spy": s["t"], "cagr": st["cagr"], "sharpe": st["sharpe"], + "mdd": _mdd(x)} + if ew is not None: + e = alpha_stats(x, ew.reindex(x.index)) + row.update(alpha_ew=e["alpha"], t_ew=e["t"]) + out[k] = row + mom = momentum_sleeve() + j = pd.concat({"x": r, "m": mom}, axis=1).dropna() + if len(j) > 500: + mm = (1 + j).resample("ME").prod() - 1 + out["corr_mom"] = float(mm.corr().iloc[0, 1]) + bl = monthly_blend(j, {"x": 0.5, "m": 0.5}) + inc = {} + for k in ("full", "design", "holdout"): + a, b = SPANS[k] + xb, xm = bl.loc[a:b], j["m"].loc[a:b] + qb = alpha_stats(xb, bench_returns(xb.index, "QQQ")) + qm = alpha_stats(xm, bench_returns(xm.index, "QQQ")) + inc[k] = {"blend_sharpe": stats((1 + xb).cumprod())["sharpe"], "mom_sharpe": stats((1 + xm).cumprod())["sharpe"], + "blend_alpha_qqq": qb["alpha"], "blend_t_qqq": qb["t"], "mom_alpha_qqq": qm["alpha"], + "blend_mdd": _mdd(xb), "mom_mdd": _mdd(xm)} + out["blend50"] = inc + if log: + TRIALS.append(out) + print(line(out), flush=True) + return out + + +def line(o: dict) -> str: + def cell(k, key="alpha_qqq", t="t_qqq"): + if k not in o: + return " -- " + return f"{o[k][key]*100:6.1f}% t{o[k][t]:5.2f}" + f = o.get("full", {}) + s = (f"{o['name']:44s} QQQ {cell('full')} | des {cell('design')} | hold {cell('holdout')} | " + f"SPY {cell('full', 'alpha_spy', 't_spy')} pre08 {cell('pre08', 'alpha_spy', 't_spy')} " + f"post08 {cell('post08', 'alpha_spy', 't_spy')} | sh {f.get('sharpe', np.nan):4.2f} " + f"mdd {f.get('mdd', np.nan)*100:4.0f}% cagr {f.get('cagr', np.nan)*100:5.1f}%") + if "alpha_ew" in f: + s += f" | vsEW {f['alpha_ew']*100:5.1f}% t{f['t_ew']:5.2f}" + if "corr_mom" in o: + b = o["blend50"] + s += (f" | corr {o['corr_mom']:4.2f} 50/50 sh {b['full']['blend_sharpe']:.2f} vs {b['full']['mom_sharpe']:.2f}" + f" hold {b['holdout']['blend_sharpe']:.2f} vs {b['holdout']['mom_sharpe']:.2f}") + return s + + +def ic_table(sig: pd.DataFrame, level: pd.DataFrame, horizons=(3, 6, 12), label: str = "", + min_groups: int = 5) -> list[dict]: + """Monthly Spearman rank IC of ``sig`` (month-end rows) against the next + h-month return of ``level``; Newey-West t with h lags.""" + rows = [] + me = sig.index + lv = level.reindex(me) + for h in horizons: + fwd = lv.shift(-h) / lv - 1 + ics = [] + for d in me: + x, y = sig.loc[d], fwd.loc[d] + ok = x.notna() & y.notna() + if ok.sum() >= min_groups: + ics.append((d, x[ok].rank().corr(y[ok].rank()))) + s = pd.Series(dict(ics)).dropna() + row = {"signal": label, "h": h, "n": len(s)} + for k, (a, b) in {"full": (None, None), "design": (None, "2016-12-31"), "holdout": ("2017-01-01", None)}.items(): + x = s.loc[a:b] + if len(x) < 24: + continue + e = (x - x.mean()).to_numpy() + X = np.ones((len(e), 1)) + from research.lab.alpha import _nw_se + se = _nw_se(X, e, lags=h)[0] + row[f"ic_{k}"] = float(x.mean()); row[f"t_{k}"] = float(x.mean() / se) + rows.append(row) + print(f" IC {label:42s} h{h:2d} n{len(s):4d} full {row.get('ic_full', np.nan):6.3f} t{row.get('t_full', np.nan):5.2f}" + f" | des {row.get('ic_design', np.nan):6.3f} t{row.get('t_design', np.nan):5.2f}" + f" | hold {row.get('ic_holdout', np.nan):6.3f} t{row.get('t_holdout', np.nan):5.2f}", flush=True) + return rows + + +def pos3y(px: pd.DataFrame, win: int = 756) -> pd.DataFrame: + lo = px.rolling(win, min_periods=win // 2).min() + hi = px.rolling(win, min_periods=win // 2).max() + return (px - lo) / (hi - lo) + + +# --------------------------------------------------------------------------- +# (a) + (b) + (c) at the ETF level +# --------------------------------------------------------------------------- + +def phase_sector_table(px: pd.DataFrame, me: pd.DatetimeIndex, raw: pd.Series) -> dict: + """Does the phase -> sector pattern repeat? Mean monthly sector return in + excess of EW9 by CFNAI phase, per period; the favoured-minus-unfavoured + spread of the Fidelity map; and the correlation of the 4 x 9 pattern + between periods.""" + m = px.ffill().reindex(me) + fwd = (m.shift(-1) / m - 1).iloc[:-1] + ex = fwd.sub(fwd.mean(axis=1), axis=0) + phases = {"real-time": asof(phase_cfnai(), me).iloc[:-1], + "held month (look-ahead)": pd.Series([raw.get((d + pd.offsets.MonthBegin(1)).normalize()) for d in me[:-1]], index=me[:-1])} + periods = {"1999-2007": ("1998-12-01", "2007-12-31"), "2008-2016": ("2008-01-01", "2016-12-31"), + "2017-2026": ("2017-01-01", None)} + out = {} + for pk, ph in phases.items(): + mats, spread = {}, {} + for per, (a, b) in periods.items(): + e, p_ = ex.loc[a:b], ph.loc[a:b] + mat = e.groupby(p_).mean().reindex(PHASES) + mats[per] = mat + sp = {} + for phase in PHASES: + fav = [x for x in ETF9 if ETF_SECTOR[x] in FAVOURED[phase]] + unf = [x for x in ETF9 if ETF_SECTOR[x] in UNFAVOURED[phase]] + rows = e[p_ == phase] + if len(rows) >= 6 and fav and unf: + diff = rows[fav].mean(axis=1) - rows[unf].mean(axis=1) + sp[phase] = {"months": len(rows), "fav_minus_unf_ann": float(diff.mean() * 12), + "t": float(diff.mean() / diff.std() * np.sqrt(len(diff)))} + spread[per] = sp + cors = {} + keys = list(mats) + for i in range(len(keys)): + for j in range(i + 1, len(keys)): + a_, b_ = mats[keys[i]].stack(), mats[keys[j]].stack() + k = a_.index.intersection(b_.index) + cors[f"{keys[i]} vs {keys[j]}"] = float(a_[k].corr(b_[k])) + out[pk] = {"spread": spread, "pattern_corr": cors, + "matrix": {k: (v * 12).round(4).to_dict() for k, v in mats.items()}} + print(f" phase -> sector pattern ({pk}): correlation between periods {cors}") + for per, sp in spread.items(): + print(" " + per + " " + " ".join(f"{ph} {v['fav_minus_unf_ann']*100:5.1f}%/yr t{v['t']:5.2f} ({v['months']}m)" + for ph, v in sp.items())) + return out + + +def etf_tests(res: dict) -> dict: + px = etf_prices().loc["1998-12-22":] + me = month_ends(px.index) + me = me[me < px.index[-1]] + W_ew = pd.DataFrame(1 / 9, index=me, columns=ETF9) + ew = etf_book(W_ew, px) + res["etf_ew9"] = evaluate("ETF EW 9 sectors (baseline)", ew, "baseline") + rets = {"ETF EW9": ew} + + ph = {"cfnai": asof(phase_cfnai(), me), "indpro": asof(phase_indpro(), me)} + for k, s in ph.items(): + W = pd.DataFrame(0.0, index=me, columns=ETF9) + for d in me: + fav = [e for e in ETF9 if ETF_SECTOR[e] in FAVOURED.get(s.loc[d], [])] or ETF9 + W.loc[d, fav] = 1 / len(fav) + name = f"(a) ETF cycle rotation, {k} phase" + rets[name] = etf_book(W, px) + res[name] = evaluate(name, rets[name], "a", ew=ew) + cnt = s.value_counts() + print(" phase months:", cnt.to_dict()) + + # Upper bound, not a trial: the final-vintage CFNAI phase of the month + # about to be held, i.e. perfect knowledge of the current phase. + raw = phase_cfnai(raw=True) + nxt = pd.Series([raw.get((d + pd.offsets.MonthBegin(1)).normalize()) for d in me], index=me) + W = pd.DataFrame(0.0, index=me, columns=ETF9) + for d in me: + fav = [e for e in ETF9 if ETF_SECTOR[e] in FAVOURED.get(nxt.loc[d], [])] or ETF9 + W.loc[d, fav] = 1 / len(fav) + name = "bound: ETF rotation with the held month's CFNAI phase (look-ahead)" + rets[name] = etf_book(W, px) + res[name] = evaluate(name, rets[name], "bound", ew=ew, log=False) + res["phase_sector_table"] = phase_sector_table(px, me, raw) + + c = px.ffill() + rf12 = (1 + tbill(px.index).fillna(0) / 252).cumprod() + m12_1 = (c.shift(21) / c.shift(252) - 1).reindex(me) + r12 = (c / c.shift(252)).reindex(me) + rf_12 = (rf12 / rf12.shift(252)).reindex(me) + lt = (c.shift(252) / c.shift(1260) - 1).reindex(me) + pos = pos3y(c).reindex(me) + + def topk(sig, k, ascending=False): + W = pd.DataFrame(np.nan, index=me, columns=ETF9) + for d in me: + x = sig.loc[d].dropna() + if len(x) < 9: + continue + pick = x.sort_values(ascending=ascending).index[:k] + W.loc[d] = 0.0 + W.loc[d, pick] = 1 / k + return W + + variants = { + "(b) ETF CS momentum 12-1 top3": topk(m12_1, 3), + "(b) ETF TS momentum 12m vs T-bill": (r12.gt(rf_12, axis=0).astype(float) / 9).where(r12.notna()), + "(b) ETF long-term reversal 60..13 bottom3": topk(lt, 3, ascending=True), + "(c) ETF 3y-range position low3 (contrarian)": topk(pos, 3, ascending=True), + "(c) ETF 3y-range position high3": topk(pos, 3), + } + for name, W in variants.items(): + rets[name] = etf_book(W.dropna(how="all"), px) + res[name] = evaluate(name, rets[name], name[1], ew=ew) + + lvl = c.reindex(me) + res["ic_etf"] = (ic_table(pos, lvl, label="ETF 3y position") + ic_table(m12_1, lvl, label="ETF mom 12-1") + + ic_table(lt, lvl, label="ETF return 60..13 (IC<0 = reversal)")) + return rets + + +# --------------------------------------------------------------------------- +# constituent level: EW PIT GICS sectors and SIC-2 industries +# --------------------------------------------------------------------------- + +def group_index(p, groups: pd.Series, min_members: int = 3) -> pd.DataFrame: + """Daily equal-weight index of each group's PIT members.""" + r = p.close.pct_change(fill_method=None).where(p.member.shift(1, fill_value=False)) + out = {} + for g in groups.dropna().unique(): + cols = [c for c in groups.index[groups == g] if c in r.columns] + n = r[cols].notna().sum(axis=1) + gr = r[cols].mean(axis=1).where(n >= min_members) + out[g] = (1 + gr.fillna(0.0)).cumprod().where(n >= min_members) + return pd.DataFrame(out) + + +def constituent_tests(p, res: dict) -> dict: + sec = p.extra["sector"].reindex(p.close.columns) + idx = group_index(p, sec) + dates = p.dates[(p.dates >= "2005-01-01")] + me = month_ends(dates) + me = me[me < dates[-1]] + members = p.member + sectors = [s for s in idx.columns] + + def book(chosen): + return stock_book(p, sector_weights(chosen, members, sec)) + + ew = book({d: {s: 1 / len(sectors) for s in sectors} for d in me}) + res["gics_ew"] = evaluate("GICS EW sectors, EW members (baseline)", ew, "baseline") + rets = {"GICS EW": ew} + ph = asof(phase_cfnai(), me) + ch = {} + for d in me: + fav = FAVOURED[ph.loc[d]] if isinstance(ph.loc[d], str) else sectors + ch[d] = {s: 1 / len(fav) for s in fav} + name = "(a) GICS cycle rotation, cfnai phase" + rets[name] = book(ch) + res[name] = evaluate(name, rets[name], "a", ew=ew) + + c = idx.ffill() + rf12 = (1 + tbill(c.index).fillna(0) / 252).cumprod() + m12_1 = (c.shift(21) / c.shift(252) - 1).reindex(me) + r12 = (c / c.shift(252)).reindex(me) + rf_12 = (rf12 / rf12.shift(252)).reindex(me) + lt = (c.shift(252) / c.shift(1260) - 1).reindex(me) + pos = pos3y(c).reindex(me) + + def topk(sig, k, ascending=False): + out = {} + for d in me: + x = sig.loc[d].dropna() + if len(x) < 9: + continue + out[d] = {s: 1 / k for s in x.sort_values(ascending=ascending).index[:k]} + return out + + ts = {} + for d in me: + x, f = r12.loc[d], rf_12.loc[d] + if x.notna().sum() >= 9: + ts[d] = {s: 1 / len(sectors) for s in x.index[x > f]} + if not ts[d]: + ts[d] = {"__cash__": 0.0} + for name, ch in {"(b) GICS CS momentum 12-1 top3": topk(m12_1, 3), + "(b) GICS TS momentum 12m vs T-bill": ts, + "(b) GICS long-term reversal 60..13 bottom3": topk(lt, 3, ascending=True), + "(c) GICS 3y-range position low3": topk(pos, 3, ascending=True)}.items(): + rets[name] = book(ch) + res[name] = evaluate(name, rets[name], name[1], ew=ew) + lvl = c.reindex(me) + res["ic_gics"] = (ic_table(pos, lvl, label="GICS 3y position") + ic_table(m12_1, lvl, label="GICS mom 12-1") + + ic_table(lt, lvl, label="GICS return 60..13 (IC<0 = reversal)")) + + # SIC 2-digit industries: more groups for the position signal + sic = pd.read_csv(SIC).drop_duplicates("ticker").set_index("ticker")["sic"] + sic2 = (pd.to_numeric(sic, errors="coerce") // 100).reindex(p.close.columns) + ind = group_index(p, sic2, min_members=4).ffill() + ipos = pos3y(ind).reindex(me) + ilvl = ind.reindex(me) + im = (ind.shift(21) / ind.shift(252) - 1).reindex(me) + res["ic_sic2"] = (ic_table(ipos, ilvl, label="SIC2 3y position", min_groups=15) + + ic_table(im, ilvl, label="SIC2 mom 12-1 (reference)", min_groups=15)) + q = {} + for d in me: + x = ipos.loc[d].dropna() + if len(x) < 15: + continue + low = x.sort_values().index[: max(3, len(x) // 5)] + q[d] = {g: 1 / len(low) for g in low} + name = "(c) SIC2 3y-range position bottom quintile" + rets[name] = stock_book(p, sector_weights(q, members, sic2)) + res[name] = evaluate(name, rets[name], "c", ew=ew) + return rets + + +# --------------------------------------------------------------------------- +# (c) EDGAR revenue cycle by sector, and the semiconductor cycle +# --------------------------------------------------------------------------- + +def revenue_cycle(p, res: dict) -> dict: + """Added after the IC table: the sign (acceleration -> lower 12-month + return) was read off 2011-2016 (IC -0.12, t -1.3), so only 2017+ tests + the portfolio rule.""" + from research.lab import fundamentals as fu + F = fu.load(p) + g = F["rev_growth"].where(p.member).clip(-0.9, 3.0) + sec = p.extra["sector"].reindex(p.close.columns) + med = pd.DataFrame({s: g[[c for c in sec.index[sec == s] if c in g.columns]].median(axis=1) + for s in sec.dropna().unique()}) + idx = group_index(p, sec).ffill() + me = month_ends(p.dates[p.dates >= "2011-06-01"]) + me = me[me < p.dates[-1]] + lvl = idx.reindex(me) + level = med.reindex(me) + accel = (med - med.shift(126)).reindex(me) + res["ic_rev_cycle"] = (ic_table(level, lvl, label="GICS median revenue growth", min_groups=9) + + ic_table(accel, lvl, label="GICS revenue growth 6m change", min_groups=9)) + ch = {} + for d in me: + x = accel.loc[d].dropna() + if len(x) >= 9: + ch[d] = {s: 1 / 3 for s in x.sort_values().index[:3]} + name = "(c) GICS 3 sectors with the lowest revenue-growth acceleration" + r = stock_book(p, sector_weights(ch, p.member, sec)) + ew_ch = {d: {s: 1 / len(idx.columns) for s in idx.columns} for d in ch} + ew = stock_book(p, sector_weights(ew_ch, p.member, sec)) + res[name] = evaluate(name, r, "c", ew=ew) + return {name: r} + + +def semis_cycle(p, res: dict) -> dict: + """^SOX vs S&P 500: does the semiconductor cycle position predict the + next 3-12 months? Plus two switch rules SMH <-> SPY.""" + sox, spx = yahoo("^SOX"), yahoo("^GSPC") + j = pd.concat({"sox": sox, "spx": spx}, axis=1).dropna() + me = month_ends(j.index) + me = me[me < j.index[-1]] + rel = j["sox"] / j["spx"] + sigs = {"pos3y": pos3y(j[["sox"]])["sox"].reindex(me), + "rel_mom12": (rel / rel.shift(252) - 1).reindex(me), + "pos3y_rel": pos3y(pd.DataFrame({"r": rel}))["r"].reindex(me)} + # SIC 3674 revenue growth acceleration (EDGAR, PIT members, 2011+) + from research.lab import fundamentals as fu + F = fu.load(p) + sic = pd.read_csv(SIC).drop_duplicates("ticker").set_index("ticker")["sic"] + semis = [t for t in sic.index[pd.to_numeric(sic, errors="coerce") == 3674] if t in p.close.columns] + rg = F["rev_growth"][semis].where(p.member[semis]).clip(-0.9, 3.0).median(axis=1) + sigs["rev_growth"] = asof(rg.dropna(), me) + sigs["rev_accel6m"] = asof((rg - rg.shift(126)).dropna(), me) + lr = np.log(rel.reindex(me)) + rows = [] + from research.lab.alpha import _nw_se + for h in (3, 6, 12): + fwd = lr.shift(-h) - lr + for k, s in sigs.items(): + for span, (a, b) in {"full": (None, None), "to2016": (None, "2016-12-31"), "2017+": ("2017-01-01", None)}.items(): + d = pd.concat({"x": s, "y": fwd}, axis=1).loc[a:b].dropna() + if len(d) < 36: + continue + x = (d["x"] - d["x"].mean()) / d["x"].std() + X = np.column_stack([np.ones(len(d)), x]) + coef, *_ = np.linalg.lstsq(X, d["y"].to_numpy(), rcond=None) + se = _nw_se(X, d["y"].to_numpy() - X @ coef, lags=h) + rows.append({"signal": k, "h": h, "span": span, "n": len(d), "slope_per_sd": coef[1], + "t": coef[1] / se[1], "start": str(d.index[0].date())}) + print(f" SOX {k:12s} h{h:2d} {span:6s} n{len(d):4d} slope/sd {coef[1]*100:6.2f}% t {coef[1]/se[1]:5.2f}", flush=True) + res["semis_regressions"] = rows + + px = pd.DataFrame({"SMH": yahoo("SMH"), "SPY": yahoo("SPY")}).dropna() + pme = month_ends(px.index) + pme = pme[pme < px.index[-1]] + rets = {} + # The last rule's sign comes from the design-period (<= 2016) regression + # above (rel_mom12 -> next 12 months: slope < 0), so only 2017+ tests it. + for name, cond in {"(c) semis switch: SMH when SOX/SPX 12m rel mom > 0": sigs["rel_mom12"] > 0, + "(c) semis switch: SMH when SOX/SPX 12m rel mom < 0 (sign from <=2016)": sigs["rel_mom12"] < 0, + "(c) semis switch: SMH when SOX 3y position < 0.3": sigs["pos3y"] < 0.3, + "(c) semis switch: SMH when SIC3674 rev accel > 0": sigs["rev_accel6m"] > 0}.items(): + cc = cond.reindex(pme, method="ffill") + if name.endswith("rev accel > 0"): + cc = cc.where(asof(sigs["rev_accel6m"], pme).notna()) + W = pd.DataFrame({"SMH": cc.astype(float), "SPY": 1 - cc.astype(float)}).where(cc.notna()).dropna() + rets[name] = etf_book(W, px) + res[name] = evaluate(name, rets[name], "c") + print(" months in SMH:", int(cc.sum()), "of", int(cc.notna().sum())) + rets["SMH buy&hold"] = px["SMH"].pct_change().loc["2000-07-01":] + res["SMH buy&hold"] = evaluate("SMH buy & hold (reference)", rets["SMH buy&hold"], "reference", log=False) + return rets + + +# --------------------------------------------------------------------------- +# (d) the momentum sleeve by macro state, and cycle-conditioned versions +# --------------------------------------------------------------------------- + +def momentum_tests(p, res: dict) -> dict: + from research.lab.rotation import factor_ranks + from research.lab.us_sweep import _blend, run + + mom = momentum_sleeve() + me = month_ends(mom.index) + st = states(me) + st_day = st.reindex(mom.index.union(me)).shift(1).ffill().reindex(mom.index) # known at last month-end + bq = bench_returns(mom.index, "QQQ") + rows = [] + for col in st.columns: + for val in st[col].dropna().unique(): + mask = st_day[col] == val + for span, (a, b) in {"design": ("2006-09-27", "2016-12-31"), "holdout": ("2017-01-01", None)}.items(): + x = mom.loc[a:b][mask.loc[a:b]] + if len(x) < 40: + continue + q = bq.reindex(x.index) + rows.append({"state": col, "value": str(val), "span": span, "days": len(x), + "share": len(x) / len(mom.loc[a:b]), + "ann_ret": float(x.mean() * 252), "sharpe": float(x.mean() / x.std() * np.sqrt(252)), + "qqq_ann_ret": float(q.mean() * 252), + "active_ann": float((x - q).mean() * 252)}) + S = pd.DataFrame(rows) + S.to_csv(os.path.join(RESULTS, "sector_cycle_states.csv"), index=False) + print(S.to_string(float_format=lambda v: f"{v:.2f}"), flush=True) + res["states"] = rows + + # cash switch in the design-worst state (>= 10% of design days) + d = S[(S["span"] == "design") & (S["share"] >= 0.10) & (S["share"] <= 0.9)].sort_values("sharpe") + worst = d.iloc[0] + print(f"design-worst state: {worst['state']}={worst['value']} (Sharpe {worst['sharpe']:.2f}, share {worst['share']:.2f})") + res["cash_switch_state"] = {"state": worst["state"], "value": worst["value"], "design_sharpe": float(worst["sharpe"])} + off = (st_day[worst["state"]].astype(str) == worst["value"]).astype(float) + rf = tbill(mom.index).fillna(0) / 252 + switch_cost = off.diff().abs().fillna(0) * 0.001 # 10 bp to sell all or buy all + rets = {} + name = f"(d) momentum, cash when {worst['state']}={worst['value']}" + rets[name] = mom * (1 - off) + rf * off - switch_cost + res[name] = evaluate(name, rets[name], "d") + res["momentum sleeve"] = evaluate("momentum sleeve (overnight_recovery_top10)", mom, "reference", log=False) + + # picks restricted by the CFNAI phase + pr = factor_ranks(p, "us", ["overnight_mom", "rec_126"]) + score = _blend(pr["overnight_mom"], pr["rec_126"]).where(p.member) + ph = asof(phase_cfnai(), p.dates) + sec = p.extra["sector"].reindex(p.close.columns) + fav = pd.DataFrame(False, index=p.dates, columns=p.close.columns) + unf = pd.DataFrame(True, index=p.dates, columns=p.close.columns) + for phase in PHASES: + rowsel = (ph == phase).to_numpy() + fav.loc[rowsel] = np.repeat(sec.isin(FAVOURED[phase]).to_numpy()[None, :], rowsel.sum(), axis=0) + unf.loc[rowsel] = np.repeat((~sec.isin(UNFAVOURED[phase])).to_numpy()[None, :], rowsel.sum(), axis=0) + base = run(p, score, [LONG], top_n=10, buffer=2.0, every=10)["_rets"] + res["momentum rerun"] = evaluate("momentum sleeve (rerun, same engine)", base, "reference", log=False) + for name, m in {"(d) momentum, favoured sectors of CFNAI phase only": fav, + "(d) momentum, unfavoured sectors of CFNAI phase removed": unf}.items(): + rets[name] = run(p, score.where(m), [LONG], top_n=10, buffer=2.0, every=10)["_rets"] + res[name] = evaluate(name, rets[name], "d") + return rets + + +def main() -> None: + if len(sys.argv) > 1 and sys.argv[1] == "fetch": + fetch() + return + from research.lab.factor_scan import load + only = sys.argv[1:] or ["etf", "constituent", "revenue", "semis", "momentum"] + res, rets = {}, {} + p = None + if set(only) - {"etf"}: + p, _ = load("us") + if "etf" in only: + print("== ETF level (a)(b)(c)", flush=True) + rets.update(etf_tests(res)) + if "constituent" in only: + print("== constituent level (a)(b)(c)", flush=True) + rets.update(constituent_tests(p, res)) + if "revenue" in only: + print("== EDGAR revenue cycle (c)", flush=True) + rets.update(revenue_cycle(p, res)) + if "semis" in only: + print("== semiconductor cycle (c)", flush=True) + rets.update(semis_cycle(p, res)) + if "momentum" in only: + print("== momentum sleeve conditioning (d)", flush=True) + rets.update(momentum_tests(p, res)) + # one results file per kind; a partial run updates only its own keys + jpath, rpath = os.path.join(RESULTS, "sector_cycle.json"), os.path.join(RESULTS, "sector_cycle_returns.pkl") + old = json.load(open(jpath)) if os.path.exists(jpath) else {} + with open(jpath, "w") as fh: + json.dump({**old, **res}, fh, indent=1, default=float) + R = pd.read_pickle(rpath) if os.path.exists(rpath) else pd.DataFrame() + new = pd.DataFrame(rets) + R = pd.concat([R.drop(columns=[c for c in new.columns if c in R.columns]), new], axis=1) + R.to_pickle(rpath) + flat = [] + for t in TRIALS: + row = {"name": t["name"], "test": t["test"], "start": t["start"]} + for k in SPANS: + for m, v in t.get(k, {}).items(): + row[f"{k}_{m}"] = v + row["corr_mom"] = t.get("corr_mom") + flat.append(row) + path = os.path.join(RESULTS, "sector_cycle_trials.csv") + pd.DataFrame(flat).to_csv(path, mode="a", header=not os.path.exists(path), index=False) + print(f"{len(TRIALS)} trials logged to {path}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/theme_propagation.py b/research/lab/theme_propagation.py new file mode 100644 index 0000000..ed5aaea --- /dev/null +++ b/research/lab/theme_propagation.py @@ -0,0 +1,268 @@ +""" +Does a main uptrend propagate along linked industries with a lag of months? +The AI build-out moved stage by stage (hyperscalers and power in H1 2025, +optics / equipment / memory in H2 2025, chips and servers in Q2 2026). This +asks whether that pattern is systematic: when an industry is in a strong +run, do the industries linked to it that have NOT moved yet outperform over +the next 1-6 months? Every era gets its own theme; the rule never names one. + + PYTHONPATH=. uv run python -m research.lab.theme_propagation + +Units: 3-digit SIC industries of PIT S&P 500 members (sec_sic.csv), equal +weight, with at least 3 members that month. Monthly, at each month end t: + + hot industries whose 6-month return over the equal-weight market is + in the top 10% of industries and above +15% + links econ_up supplier industries of the hot one: 10-K customer + sentences filed in the last 3 years (filer = supplier, + named company = customer; research/lab/fetch_supply_links.py) + econ_any any 10-K mention in either direction, last 3 years + corr the 5 industries whose monthly excess returns correlated + most with it over the last 36 months + lagged the 5 industries whose excess return in months m+1..m+3 + correlated most with its return in month m, last 60 months + followers linked industries that are not hot themselves and whose own + 6-month excess return is below the cross-industry median + ("the linked industry has not moved yet"); variant: all linked + control below-median industries not linked to any hot industry + +Portfolio: equal weight over the member stocks of the follower industries, +bought at the first open after t and held 1, 3 or 6 months (overlapping +cohorts averaged, i.e. monthly tranches), against the equal-weight market. +20 bp per unit of turnover is charged at each rebalance. Reported per era +(2006-10, 2011-15, 2016-20, 2021-26) and design (<=2016) / holdout. +A stock-level check follows: the lead_lag link signals over 63- and 126-day +windows predicting the next 63 days (Fama-MacBeth with the usual controls). +""" + +from __future__ import annotations + +import json +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factor_scan import RESULTS, load + +warnings.filterwarnings("ignore") +ERAS = [("2006", "2010"), ("2011", "2015"), ("2016", "2020"), ("2021", "2026")] +HOT_Q, HOT_MIN, LOOK = 0.90, 0.15, 6 + + +def industry_map(p) -> pd.Series: + sic = pd.read_csv(os.path.join("data", "research", "sec_sic.csv")).drop_duplicates("ticker").set_index("ticker")["sic"] + s = (pd.to_numeric(sic, errors="coerce") // 10).dropna().astype(int) + return s.reindex(p.close.columns) + + +def monthly_open_returns(p) -> tuple[pd.DataFrame, pd.DatetimeIndex]: + """Stock return from the first open of month m+1 to the first open of + month m+2, indexed by the signal month end m.""" + o = p.open.where(p.open.notna(), p.close).ffill() + first = pd.Series(p.dates, index=p.dates).groupby(p.dates.to_period("M")).min() + ends = pd.Series(p.dates, index=p.dates).groupby(p.dates.to_period("M")).max() + fo = o.loc[first.values] + r = fo.shift(-1) / fo - 1 # first open of month k -> first open of month k+1 + r.index = ends.values + r = r.shift(-1) # signal at end of m -> return over month m+1 + return r, pd.DatetimeIndex(ends.values) + + +def group_panels(p, ind: pd.Series): + """Monthly member masks and equal-weight monthly close-to-close excess + returns per industry (for the signals).""" + ends = pd.Series(p.dates, index=p.dates).groupby(p.dates.to_period("M")).max().values + c = p.close.ffill() + mc = c.loc[ends] + mr = (mc / mc.shift(1) - 1).where(p.member.loc[ends].shift(1).fillna(False).astype(bool)) + mem = p.member.loc[ends] + mkt = mr.mean(axis=1) + groups = sorted(ind.dropna().unique()) + G = pd.DataFrame(index=mr.index, columns=groups, dtype=float) + N = pd.DataFrame(index=mr.index, columns=groups, dtype=float) + for g in groups: + cols = list(ind.index[ind == g]) + G[g] = mr[cols].mean(axis=1) + N[g] = mem[cols].sum(axis=1) + X = G.sub(mkt, axis=0) + return X.where(N >= 3), N, mem + + +def econ_links(ind: pd.Series, ends, role: str) -> dict: + """{month end: set of (hot industry, linked industry)} from 10-Ks filed + in the last 3 years. role "up": customer industry -> supplier industry.""" + d = pd.read_csv(os.path.join("data", "research", "sec_10k_links.csv"), parse_dates=["filed"]) + if role == "up": + d = d[d["n_customer"] > 0] + d = d.assign(fi=d["filer"].map(ind), ti=d["target"].map(ind)).dropna(subset=["fi", "ti"]) + d = d[d["fi"] != d["ti"]] + out = {} + for t in ends: + live = d[(d["filed"] <= t) & (d["filed"] > t - pd.Timedelta(days=3 * 365))] + if role == "up": # filer supplies the named customer: customer industry leads + pairs = set(zip(live["ti"].astype(int), live["fi"].astype(int))) + else: + pairs = set(zip(live["ti"].astype(int), live["fi"].astype(int))) | set(zip(live["fi"].astype(int), live["ti"].astype(int))) + out[t] = pairs + return out + + +def stat_links(X: pd.DataFrame, kind: str, k: int = 5) -> dict: + out = {} + for i, t in enumerate(X.index): + win = 36 if kind == "corr" else 60 + if i < win: + out[t] = set() + continue + W = X.iloc[i - win + 1:i + 1] + W = W.loc[:, W.notna().sum() >= int(win * 0.8)] + if W.shape[1] < 10: + out[t] = set() + continue + Z = (W - W.mean()) / W.std() + Z = Z.fillna(0.0) + if kind == "corr": + C = Z.T @ Z / len(Z) + else: # C[a, b] = mean over lags 1-3 of corr(a in month m, b in month m+lag) + C = sum(Z.iloc[:-lag].T.to_numpy() @ Z.iloc[lag:].to_numpy() / (len(Z) - lag) for lag in (1, 2, 3)) / 3 + C = pd.DataFrame(C, index=Z.columns, columns=Z.columns) + C = C.to_numpy(copy=True) + np.fill_diagonal(C, -np.inf) + cols = list(Z.columns) + pairs = set() + for a in range(len(cols)): + for b in np.argsort(-C[a])[:k]: + pairs.add((cols[a], cols[b])) + out[t] = pairs + return out + + +def followers(X: pd.DataFrame, links: dict, laggard_only: bool = True) -> tuple[dict, dict, dict]: + """Per month end: hot industries, follower industries, control industries.""" + cum = np.log1p(X.fillna(0.0)).rolling(LOOK, min_periods=LOOK).sum().where(X.notna()) + hot_d, fol_d, ctl_d = {}, {}, {} + for t in X.index: + row = cum.loc[t].dropna() + if len(row) < 20: + continue + thr = max(row.quantile(HOT_Q), np.log1p(HOT_MIN)) + hot = set(row.index[row >= thr]) + med = row.median() + lag = set(row.index[row < med]) + pairs = links.get(t, set()) + linked = {b for a, b in pairs if a in hot} - hot + linked &= set(row.index) + fol = linked & lag if laggard_only else linked + hot_d[t], fol_d[t] = hot, fol + ctl_d[t] = lag - linked - hot + return hot_d, fol_d, ctl_d + + +def portfolio(sel: dict, ind: pd.Series, mem: pd.DataFrame, R: pd.DataFrame, hold: int) -> pd.Series: + """Monthly return of equal-weight member stocks in the selected industries, + held ``hold`` months as overlapping cohorts, net of 20 bp x turnover.""" + months = R.index + W = pd.DataFrame(0.0, index=months, columns=R.columns) + for t, gs in sel.items(): + if t not in W.index or not gs: + continue + cols = [c for c in ind.index[ind.isin(gs)] if mem.at[t, c]] + if cols: + W.loc[t, cols] = 1.0 / len(cols) + Wh = sum(W.shift(k).fillna(0.0) for k in range(hold)) / hold # tranches formed at t, t-1, ... + invested = Wh.sum(axis=1) + gross = (Wh * R.fillna(0.0)).sum(axis=1) / invested.replace(0, np.nan) + turn = (Wh.div(invested.replace(0, np.nan), axis=0).fillna(0.0).diff().abs().sum(axis=1)) / 2 + return (gross - 0.002 * 2 * turn).where(invested > 0) + + +def mstats(r: pd.Series, b: pd.Series) -> dict: + """Excess over the benchmark: mean per month, annualized, Newey-West t (3 lags).""" + x = (r - b).dropna() + if len(x) < 12: + return {"ann": np.nan, "t": np.nan, "n": len(x)} + n, m = len(x), x.mean() + e = x - m + s = (e ** 2).sum() / n + for l in (1, 2, 3): + s += 2 * (1 - l / 4) * (e.iloc[l:].to_numpy() * e.iloc[:-l].to_numpy()).sum() / n + return {"ann": float(m * 12), "t": float(m / np.sqrt(s / n)), "n": n} + + +def stock_level_long(p) -> dict: + """lead_lag link signals over long windows, next-63-day Fama-MacBeth.""" + from research.lab import lead_lag as ll + from research.lab.rotation import factor_ranks + pr = factor_ranks(p, "us", ["overnight_mom", "mom_12_1"]) + L = pickle.load(open(ll.CACHE, "rb")) + c = p.close.ffill() + out = {} + for win in (63, 126): + own = (c / c.shift(win) - 1).where(p.member) + for key in ("corr_peer", "10k_any", "10k_customer"): + for kind in ("mom", "gap"): + sig = ll.link_signal(p, L[key], win, kind) + fm = ll.fama_macbeth(p, sig, {"own": own, "mom_12_1": pr["mom_12_1"], "overnight": pr["overnight_mom"]}, 63) + out[f"{key}_{kind}_{win}d"] = fm + print(f" {key}_{kind} window {win}d -> next 63d: FM t all {fm['all']['t']:5.2f} | design {fm['design']['t']:5.2f} " + f"| holdout {fm['holdout']['t']:5.2f} | IC {fm['all']['ic']:+.3f}", flush=True) + return out + + +def main() -> None: + p, _ = load("us") + ind = industry_map(p) + X, N, mem = group_panels(p, ind) + R, ends = monthly_open_returns(p) # unmasked: a stock bought at t keeps its return after leaving the index + mkt = R.where(mem).mean(axis=1) + from research.lab.alpha import bench_returns + q = bench_returns(p.dates) + qm = (1 + q.fillna(0.0)).groupby(p.dates.to_period("M")).prod() - 1 + qm.index = ends + qm = qm.shift(-1) # align with R (return over month m+1) + print(f"{X.notna().any().sum()} industries with >=3 members at some point; median per month " + f"{int(X.notna().sum(axis=1).median())}", flush=True) + link_sets = {"econ_up": econ_links(ind, X.index, "up"), "econ_any": econ_links(ind, X.index, "any"), + "corr": stat_links(X, "corr"), "lagged": stat_links(X, "lagged")} + res, rets = {}, {} + for name, L in link_sets.items(): + for lag_only in (True, False): + hot, fol, ctl = followers(X, L, lag_only) + nf = np.mean([len(v) for v in fol.values()]) + tag = f"{name}{'' if lag_only else ' (all linked)'}" + for hold in (1, 3, 6): + r = portfolio(fol, ind, mem, R, hold) + rc = portfolio(ctl, ind, mem, R, hold) + row = {"followers_per_month": float(nf)} + for lab, a, b in [("all", "2006-09", "2026-08"), ("design", "2006-09", "2016-12"), ("holdout", "2017-01", "2026-08")] + \ + [(f"{e0}-{e1}", e0, e1) for e0, e1 in ERAS]: + row[lab] = mstats(r.loc[a:b], mkt.loc[a:b]) + row[lab + "_vs_control"] = mstats(r.loc[a:b], rc.loc[a:b]) + row["vs_qqq_holdout"] = mstats(r.loc["2017":"2026-08"], qm.loc["2017":"2026-08"]) + res[f"{tag} hold{hold}"] = row + rets[f"{tag} hold{hold}"] = r + eras = " ".join(f"{k}: {row[k]['ann']*100:+5.1f}% (t {row[k]['t']:+.2f})" for k in [f"{a}-{b}" for a, b in ERAS]) + print(f"{tag:22s} hold {hold}m n {nf:4.1f} | vs mkt design {row['design']['ann']*100:+5.1f}% t {row['design']['t']:+.2f} " + f"holdout {row['holdout']['ann']*100:+5.1f}% t {row['holdout']['t']:+.2f} | vs control holdout " + f"{row['holdout_vs_control']['ann']*100:+5.1f}% t {row['holdout_vs_control']['t']:+.2f} | {eras}", flush=True) + # which themes fired: most frequent hot industries per era + hot, _, _ = followers(X, link_sets["corr"], True) + sic = pd.read_csv(os.path.join("data", "research", "sec_sic.csv")) + desc = sic.assign(g=pd.to_numeric(sic["sic"], errors="coerce") // 10).groupby("g")["sicDescription"].first() + themes = {} + for e0, e1 in ERAS: + cnt = pd.Series([g for t, hs in hot.items() if e0 <= str(t.year) <= e1 for g in hs]).value_counts().head(6) + themes[f"{e0}-{e1}"] = {f"{int(g)} {desc.get(g, '')}": int(n) for g, n in cnt.items()} + print(f"hot industries {e0}-{e1}: " + ", ".join(f"{k} ({v}m)" for k, v in themes[f"{e0}-{e1}"].items()), flush=True) + print("\nstock level, long windows:", flush=True) + res["_stock_level"] = stock_level_long(p) + res["_themes"] = themes + json.dump(res, open(os.path.join(RESULTS, "theme_propagation.json"), "w"), indent=1, default=float) + pd.DataFrame(rets).to_pickle(os.path.join(RESULTS, "theme_propagation_returns.pkl")) + + +if __name__ == "__main__": + main()