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 <noreply@anthropic.com>
This commit is contained in:
2026-09-30 20:25:08 +08:00
parent db1ce69702
commit ab3acfe4cb
9 changed files with 3061 additions and 0 deletions

217
research/lab/attention.py Normal file
View File

@@ -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()

View File

@@ -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()

View File

@@ -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"(?<![\w&-])(?:" + "|".join(re.escape(a) for a in alts) + r")(?![\w-])")
def links_in(text: str, own: set[str]) -> 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()

View File

@@ -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 <https://en.wikipedia.org/> . }""")
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 <https://en.wikipedia.org/> . }""")
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()

409
research/lab/lead_lag.py Normal file
View File

@@ -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:
<link>_mom mean return of the linked names over the last 21 trading
days (5 for lead_peer).
<link>_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()

View File

@@ -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()

View File

@@ -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()

View File

@@ -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()

View File

@@ -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()