From ad569f1cea563d8584715c8f88d3b1395fdc1f7e Mon Sep 17 00:00:00 2001 From: Gahow Wang Date: Tue, 29 Sep 2026 13:40:15 +0800 Subject: [PATCH] feat: us_momentum_qqq_trend; alpha-vs-QQQ research across new data sources Register us_momentum_qqq_trend: half overnight_recovery_top10, half a QQQ trend sleeve (strategies/qqq_trend.py: 200-day switch times a 30% vol target, held as TQQQ + SHY). Mixed stock+ETF, US only (extra tickers QQQ, SHY, TQQQ). Research (research/lab/): alpha against QQQ with a 2006-16 design / 2017-26 holdout split (alpha.py, qqq_scan.py, beat_qqq.py); leverage and drawdown frontier (levered_book.py, us_target50.py); synthetic 3x ETFs from 1990 and TrendRider/QQQ-trend stress tests (fetch_long_history.py, letf_history.py, vol_trend.py); short-horizon sleeves, diversifiers, industry momentum, earnings-announcement returns; new point-in-time data sources: SEC Form 4 insider trades (fetch_insider.py, insider.py), 10-K/10-Q text similarity, 8-K earnings dates and SIC codes (fetch_lazy_prices.py, lazy_prices.py), and an S&P 400/600 universe from iShares holdings (fetch_us_mid.py, panel_mid.py). Co-Authored-By: Claude Opus 5.5 (1M context) --- CLAUDE.md | 19 +- bridge.py | 6 + research/lab/alpha.py | 71 ++++ research/lab/beat_qqq.py | 139 ++++++++ research/lab/diversifiers.py | 140 ++++++++ research/lab/earnings.py | 124 +++++++ research/lab/fetch_insider.py | 372 ++++++++++++++++++++ research/lab/fetch_lazy_prices.py | 527 +++++++++++++++++++++++++++++ research/lab/fetch_long_history.py | 107 ++++++ research/lab/fetch_us_mid.py | 462 +++++++++++++++++++++++++ research/lab/industry.py | 83 +++++ research/lab/insider.py | 432 +++++++++++++++++++++++ research/lab/lazy_prices.py | 336 ++++++++++++++++++ research/lab/letf_history.py | 185 ++++++++++ research/lab/levered_book.py | 172 ++++++++++ research/lab/panel_mid.py | 304 +++++++++++++++++ research/lab/qqq_scan.py | 143 ++++++++ research/lab/short_horizon.py | 106 ++++++ research/lab/us_target50.py | 99 ++++++ research/lab/vol_trend.py | 125 +++++++ strategies/qqq_trend.py | 56 +++ strategies/sleeve_portfolio.py | 27 ++ tests/test_sleeve_portfolio.py | 49 ++- trader.py | 6 +- 24 files changed, 4086 insertions(+), 4 deletions(-) create mode 100644 research/lab/alpha.py create mode 100644 research/lab/beat_qqq.py create mode 100644 research/lab/diversifiers.py create mode 100644 research/lab/earnings.py create mode 100644 research/lab/fetch_insider.py create mode 100644 research/lab/fetch_lazy_prices.py create mode 100644 research/lab/fetch_long_history.py create mode 100644 research/lab/fetch_us_mid.py create mode 100644 research/lab/industry.py create mode 100644 research/lab/insider.py create mode 100644 research/lab/lazy_prices.py create mode 100644 research/lab/letf_history.py create mode 100644 research/lab/levered_book.py create mode 100644 research/lab/panel_mid.py create mode 100644 research/lab/qqq_scan.py create mode 100644 research/lab/short_horizon.py create mode 100644 research/lab/us_target50.py create mode 100644 research/lab/vol_trend.py create mode 100644 strategies/qqq_trend.py diff --git a/CLAUDE.md b/CLAUDE.md index 310a91c..acb5ab4 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -52,6 +52,20 @@ PYTHONPATH=. uv run python -m research.lab.rotation --market us # walk- PYTHONPATH=. uv run python -m research.lab.persistence # does last year's best combo persist? PYTHONPATH=. uv run --with numba python -m research.lab.us_sweep --part overlays # construction / overlay sweeps PYTHONPATH=. uv run --with numba python -m research.lab.us_final # 3-sleeve US portfolio check +PYTHONPATH=. uv run --with numba python -m research.lab.short_horizon # short-horizon / residual sleeves, concentration +PYTHONPATH=. uv run python -m research.lab.levered_book # vol-managed leverage, beta-hedged alpha +PYTHONPATH=. uv run python -m research.lab.fetch_long_history [--no-proxy] # indices from 1927/1985, ETFs, BTC, FRED yields +PYTHONPATH=. uv run python -m research.lab.letf_history # synthetic 3x ETFs from 1990; TrendRider stress test +PYTHONPATH=. uv run python -m research.lab.vol_trend # vol-targeted leveraged trend, walk-forward +PYTHONPATH=. uv run python -m research.lab.diversifiers # multi-asset ETF trend, bitcoin trend +PYTHONPATH=. uv run python -m research.lab.us_target50 # the book 2006-2026 and the leverage 50%/yr needs +PYTHONPATH=. uv run --with numba python -m research.lab.qqq_scan # alpha vs QQQ of every sleeve / factor / mega-cap variant +PYTHONPATH=. uv run --with numba python -m research.lab.beat_qqq # sleeves picked on 2006-16, judged 2017-26; odds of beating QQQ +uv run python -m research.lab.fetch_insider && PYTHONPATH=. uv run python -m research.lab.insider # SEC Form 4 insider factors +PYTHONPATH=. uv run python -m research.lab.fetch_lazy_prices && PYTHONPATH=. uv run python -m research.lab.lazy_prices # 10-K/10-Q text similarity, 8-K earnings dates, SIC +PYTHONPATH=. uv run --with numba python -m research.lab.earnings # earnings-announcement return (needs sec_8k_earnings.csv) +PYTHONPATH=. uv run --with numba python -m research.lab.industry # industry momentum on SIC groups +PYTHONPATH=. uv run python -m research.lab.fetch_us_mid holdings|prices|status # PIT S&P 400/600 from iShares holdings (research/lab/panel_mid.py) # Setup uv sync # Install/sync dependencies @@ -76,7 +90,7 @@ uv run python -m unittest discover tests **Backtest engine** (`main.py`): Orchestrates data loading, strategy execution, and visualization. `backtest()` has two modes. Close mode (no `open_data`, used by `main.py` and old research scripts) is vectorized and trades at the close the signal saw. Live mode (`open_data` given, used by the web UI through `backtest_runner.py`) keeps the book like `trader.py`: same signal timing (`strategy_positions`), holdings that drift, the $50/0.1% order floor, `trader.RESIZE_BAND`, fees on the orders placed (`simulate_live_book`). -**Daily trader** (`trader.py`): Live/forward-testing system with persistent portfolio state in `data/trader_{market}_{strategy}.json`. The `auto` subcommand runs both signal generation and execution in a single invocation — designed for cron. The `simulate` subcommand replays a date range day-by-day with realistic portfolio tracking (fractional shares, cash, commissions). Orders: whole A-share lots are allocated nearest-lot within cash (`_allocate_lots`); a held name is resized only when ≥ `RESIZE_BAND` (10%) off target. Available strategies: `recovery_mom_top10`, `recovery_mom_top20`, `momentum`, `momentum_quality`, `dual_momentum`, `inverse_vol`, `trend_following`, `buy_and_hold`, 32 factor-combo strategies (`fc_{signal}_{freq}` — see `strategies/factor_combo.py`), the robust-momentum bar strategies `overnight_recovery_top10`, `lowrisk_momentum_top15`, `lowbeta_sharpe_mom_top20`, `lowrisk_top20`, and the US sleeve portfolio `us_three_sleeve`. +**Daily trader** (`trader.py`): Live/forward-testing system with persistent portfolio state in `data/trader_{market}_{strategy}.json`. The `auto` subcommand runs both signal generation and execution in a single invocation — designed for cron. The `simulate` subcommand replays a date range day-by-day with realistic portfolio tracking (fractional shares, cash, commissions). Orders: whole A-share lots are allocated nearest-lot within cash (`_allocate_lots`); a held name is resized only when ≥ `RESIZE_BAND` (10%) off target. Available strategies: `recovery_mom_top10`, `recovery_mom_top20`, `momentum`, `momentum_quality`, `dual_momentum`, `inverse_vol`, `trend_following`, `buy_and_hold`, 32 factor-combo strategies (`fc_{signal}_{freq}` — see `strategies/factor_combo.py`), the robust-momentum bar strategies `overnight_recovery_top10`, `lowrisk_momentum_top15`, `lowbeta_sharpe_mom_top20`, `lowrisk_top20`, and the US sleeve portfolios `us_three_sleeve` and `us_momentum_qqq_trend`. **Live books** (`trader.LIVE_BOOKS`): a book is a cohort of strategies, each with its own account, started on one date with one capital and a flat fee per order. The original book (every live strategy, `MONITOR_CAPITAL`, `MARKET_FEES`) keeps state in `data/trader_*.json`; a named book keeps it in `data/books//`. The monitor runs every book; `auto`/`status`/`log`/`compare`/`replay` and `bridge.py live` take `--book`, and `bridge.py books` lists them for the UI's book selector. `selected_2026_09` (29 strategies, US $10,000 / CN ¥100,000, $2 / ¥10 per order) starts on 2026-09-25 (US) and 2026-09-28 (CN). @@ -94,7 +108,8 @@ uv run python -m unittest discover tests - `adaptive_momentum.py` — Momentum weighted by inverse volatility - `recovery_momentum.py` — Recovery (price/63d low) + 12-1mo momentum composite. Best US performer. - `robust_momentum.py` — Bar strategies from the 2026-09 PIT research (report: https://claude.ai/artifact/Mo9nrzs5SWj5FfyshERyWS): US `overnight_recovery` (overnight-return momentum + distance from 126d low), CN `lowrisk_momentum` (amplitude-cut momentum + low beta + price-volume corr + momentum/vol) and close-only `lowbeta_sharpe_mom`. Top-N with rank hysteresis, rebalance every 10 days, weights drift between rebalances. `low_risk` (US `lowrisk_top20`: idiosyncratic vol + low beta, top 20, monthly) is the defensive sleeve of the 2026-09-29 US research (report: https://claude.ai/artifact/VGhZvKkS9UBYZkLeqGjcNC). -- `sleeve_portfolio.py` — Several strategies in one account with fixed capital shares. `us_three_sleeve` = 1/3 each of `overnight_recovery_top10`, `lowrisk_top20`, `trend_rider_v7_vt24` (a mixed stock+ETF strategy); sleeve shares drift with their returns and reset at a month start once any is >25% off. 2011-2026 research: Sharpe 1.21 vs 0.88-1.03 for the sleeves alone. Research finding behind it: choosing factors by their recent record (rotation, factor momentum, last year's winner) did not beat a fixed composite; construction knobs (N, weighting, dynamic N) barely matter. +- `sleeve_portfolio.py` — Several strategies in one account with fixed capital shares. `us_three_sleeve` = 1/3 each of `overnight_recovery_top10`, `lowrisk_top20`, `trend_rider_v7_vt24` (a mixed stock+ETF strategy); sleeve shares drift with their returns and reset at a month start once any is >25% off. 2011-2026 research: Sharpe 1.21 vs 0.88-1.03 for the sleeves alone. Research finding behind it: choosing factors by their recent record (rotation, factor momentum, last year's winner) did not beat a fixed composite; construction knobs (N, weighting, dynamic N) barely matter. On synthetic 3x ETFs from 1990 (`research/lab/letf_history.py`), TrendRider V7's 2016-2026 record does not hold: 1990-99 Sharpe 0.2-0.3, 2003-07 negative, 1990-2026 about 0.7 with MaxDD -58% (vt24) to -69% (vt36). +- `qqq_trend.py` — `QQQVolTrend`: QQQ exposure = 1{QQQ > 200-day average} × min(1.5, 30% / EWMA vol), held as TQQQ (E/3) + SHY. 1990-2026 on synthetic ETFs: 16.0% vs QQQ 14.9%, MaxDD -37% vs -83%; 2011-2026 on real ETFs it trails QQQ (17.5% vs 19.1%) with a smaller drawdown. `us_momentum_qqq_trend` (in `sleeve_portfolio.py`) holds it half and half with `overnight_recovery_top10`. Research behind it (`research/lab/qqq_scan.py`, `beat_qqq.py`): measured as alpha against QQQ with a 2006-16 design / 2017-26 holdout split, momentum top 10 is the only stock selection whose alpha held (+10%/yr, t 1.4); single factors, walk-forward ML, a mega-cap universe, SEC insider trades, 10-K/10-Q text similarity, earnings-announcement returns, industry momentum and a PIT S&P 400/600 universe (`panel_mid.py`) showed none. In the web engine (2016-09 to 2026-09, $2/order) `us_momentum_qqq_trend` earns 21.0% vs QQQ 20.9%, MaxDD -32%. - `factor_combo.py` — Configurable factor-combination strategies with daily/weekly/biweekly/monthly rebalancing. US champions: `rec_mfilt+deep_upvol`, `ma200+mom7m+rec126`, `rec_mfilt+ma200`, `mom7m+rec126`. CN champions: `up_cap+quality_mom`, `down_resil+qual_mom`, `rec63+mom_gap`, `up_cap+mom_gap` — these were mined on today's CSI 300 members; on PIT members they earn ~3-7% a year (2016-2026). All registered in trader.py as `fc_{signal}_{freq}` (e.g., `fc_rec_mfilt_deep_upvol_monthly`). 32 new strategies total. **Metrics** (`metrics.py`): Standalone functions for portfolio analytics (Sharpe, Sortino, Calmar, max drawdown, etc.). `summary()` prints a formatted report and returns a dict. diff --git a/bridge.py b/bridge.py index 1bbe93c..966f735 100644 --- a/bridge.py +++ b/bridge.py @@ -466,6 +466,12 @@ STRATEGY_META = { "params": {}, "markets": ["us"], }, + "us_momentum_qqq_trend": { + "label": "US Momentum + QQQ Trend (1/2 each)", + "category": "tactical_allocation", + "params": {}, + "markets": ["us"], + }, } diff --git a/research/lab/alpha.py b/research/lab/alpha.py new file mode 100644 index 0000000..08f7d0b --- /dev/null +++ b/research/lab/alpha.py @@ -0,0 +1,71 @@ +""" +Alpha against a benchmark (QQQ by default): the yardstick for "does this +strategy add anything a Nasdaq-100 index fund does not?" + + from research.lab.alpha import bench_returns, alpha_stats + b = bench_returns(r.index) # QQQ total return, daily + s = alpha_stats(r, b) # alpha, beta, t, IR, ... + +alpha/beta come from the daily regression of excess returns (over the +T-bill) on the benchmark's excess return; the alpha t-stat uses Newey-West +standard errors (10 lags). IR is the annualized mean over the standard +deviation of the active return r - b. QQQ comes from +data/research/long_history.pkl (research/lab/fetch_long_history.py). +""" + +from __future__ import annotations + +import os +import pickle + +import numpy as np +import pandas as pd + +from research.lab.levered_book import tbill + +LONG_HISTORY = os.path.join("data", "research", "long_history.pkl") + + +def bench_returns(index: pd.DatetimeIndex, ticker: str = "QQQ") -> pd.Series: + with open(LONG_HISTORY, "rb") as fh: + px = pickle.load(fh)["close"][ticker].dropna() + return px.reindex(index.union(px.index)).ffill().pct_change(fill_method=None).reindex(index) + + +def _nw_se(x: np.ndarray, e: np.ndarray, lags: int = 10) -> np.ndarray: + """Newey-West standard errors of OLS coefficients.""" + n = len(e) + xtx_inv = np.linalg.inv(x.T @ x) + u = x * e[:, None] + s = u.T @ u + for l in range(1, lags + 1): + w = 1 - l / (lags + 1) + g = u[l:].T @ u[:-l] + s += w * (g + g.T) + return np.sqrt(np.diag(xtx_inv @ s @ xtx_inv)) + + +def alpha_stats(r: pd.Series, b: pd.Series, rf: pd.Series | None = None) -> dict: + df = pd.DataFrame({"r": r, "b": b}).dropna() + if rf is None: + rf = tbill(df.index) + d = rf.reindex(df.index).fillna(0.0) / 252 + y, xb = (df["r"] - d).to_numpy(), (df["b"] - d).to_numpy() + X = np.column_stack([np.ones(len(y)), xb]) + coef, *_ = np.linalg.lstsq(X, y, rcond=None) + se = _nw_se(X, y - X @ coef) + active = df["r"] - df["b"] + yrs = len(df) / 252 + cagr = lambda x: (1 + x).prod() ** (1 / yrs) - 1 # noqa: E731 + yr_r = (1 + df["r"]).groupby(df.index.year).prod() + yr_b = (1 + df["b"]).groupby(df.index.year).prod() + return {"alpha": coef[0] * 252, "t": coef[0] / se[0], "beta": coef[1], + "ir": active.mean() / active.std() * np.sqrt(252), + "cagr": cagr(df["r"]), "bench_cagr": cagr(df["b"]), + "yrs_beat": float((yr_r > yr_b).mean())} + + +def fmt_alpha(name: str, s: dict) -> str: + return (f"{name:44s} alpha {s['alpha']*100:6.1f}%/yr t {s['t']:5.2f} beta {s['beta']:4.2f} " + f"IR {s['ir']:5.2f} cagr {s['cagr']*100:5.1f}% vs {s['bench_cagr']*100:5.1f}% " + f"yrs beat {s['yrs_beat']*100:3.0f}%") diff --git a/research/lab/beat_qqq.py b/research/lab/beat_qqq.py new file mode 100644 index 0000000..d7bdd5f --- /dev/null +++ b/research/lab/beat_qqq.py @@ -0,0 +1,139 @@ +""" +Can a US book beat QQQ with alpha that holds out of sample? Assembles the +sleeves that showed alpha against QQQ and checks them on 2017-2026 with +choices made on 2006-2016 only. + + PYTHONPATH=. uv run --with numba python -m research.lab.qqq_scan # factor alphas + PYTHONPATH=. uv run python -m research.lab.vol_trend # QQQ trend grid + PYTHONPATH=. uv run --with numba python -m research.lab.beat_qqq + +Sleeves (all long-only, 10 bp costs, next-open fills): + momentum overnight-return momentum + distance from 126d low, top 10, every + 10 days, buffer 2 (the live overnight_recovery_top10) + quality mean rank of the K fundamental factors with the highest design- + period (2006-2016) alpha t-stat against QQQ in qqq_scan_factors.csv, + top 20, monthly, buffer 2. K and the list are fixed before the + holdout is looked at. + qqq_trend QQQ exposure = 1{QQQ > 200-day average} x min(1.5, 30% / EWMA vol), + held as TQQQ + 2-year Treasuries (research/lab/vol_trend.py; the + 200-day rule and a vol target near QQQ's own are the textbook + choices, not fitted) +Books are equal capital, reset monthly. Everything is reported against QQQ: +alpha (Newey-West t), CAGR per half, MaxDD, and the calendar years won. +""" + +from __future__ import annotations + +import json +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab import fundamentals as fu +from research.lab.alpha import alpha_stats, bench_returns, fmt_alpha +from research.lab.factor_scan import RESULTS, load +from research.lab.rotation import factor_ranks +from research.lab.us_final import monthly_blend +from research.lab.us_sweep import _blend, run, sector_rank + +warnings.filterwarnings("ignore") +LONG = ("2006-09-27", "2026-09-24") +DESIGN_END = "2016-12-31" + + +def pick_fundamentals(k: int = 5) -> list[tuple[str, bool]]: + """(factor, sector-neutral) with the best design-period alpha t vs QQQ, + one variant per factor.""" + d = pd.read_csv(os.path.join(RESULTS, "qqq_scan_factors.csv")) + d = d[d["name"].str.startswith("f:")].sort_values("design_t", ascending=False) + out, seen = [], set() + for name in d["name"]: + base = name[2:].replace(" (sector)", "") + if base in seen: + continue + seen.add(base) + out.append((base, name.endswith("(sector)"))) + if len(out) == k: + break + return out + + +def sleeves(p) -> tuple[pd.DataFrame, list]: + pr = factor_ranks(p, "us") + F = fu.load(p) + picks = pick_fundamentals() + ranks = [sector_rank(F[f], p) if sec else F[f].where(p.member).rank(axis=1, pct=True) for f, sec in picks] + quality = _blend(*ranks) + mom = _blend(pr["overnight_mom"], pr["rec_126"]) + R = {"momentum": run(p, mom.where(p.member), [LONG], top_n=10, buffer=2.0, every=10)["_rets"], + "quality": run(p, quality.where(p.member), [LONG], top_n=20, buffer=2.0, every=21)["_rets"]} + with open(os.path.join(RESULTS, "vol_trend.pkl"), "rb") as fh: + R["qqq_trend"] = pickle.load(fh)["returns"]["QQQ ma200 vt30 cap1.5"].reindex(p.dates) + return pd.DataFrame(R).loc[LONG[0]:LONG[1]].fillna(0.0), picks + + +def show(name: str, r: pd.Series, b: pd.Series) -> dict: + a, d, h = alpha_stats(r, b), alpha_stats(r.loc[:DESIGN_END], b), alpha_stats(r.loc["2017":], b) + e = (1 + r).cumprod() + mdd = float((e / e.cummax() - 1).min()) + print(f"{fmt_alpha(name, a)} mdd {mdd*100:4.0f}% | 06-16 a {d['alpha']*100:5.1f}% t {d['t']:5.2f} " + f"cagr {d['cagr']*100:5.1f} vs {d['bench_cagr']*100:5.1f} | 17-26 a {h['alpha']*100:5.1f}% " + f"t {h['t']:5.2f} cagr {h['cagr']*100:5.1f} vs {h['bench_cagr']*100:5.1f}") + return {"full": a, "design": d, "holdout": h, "mdd": mdd} + + +def main() -> None: + p, _ = load("us") + R, picks = sleeves(p) + b = bench_returns(R.index) + print("fundamentals picked on 2006-2016 alpha t:", picks) + m = (1 + R).resample("ME").prod() - 1 + print("monthly correlation:\n", m.corr().round(2).to_string()) + out = {"picks": picks} + for c in R.columns: + out[c] = show(c, R[c], b) + books = {"momentum + qqq_trend": {"momentum": .5, "qqq_trend": .5}, + "momentum + quality": {"momentum": .5, "quality": .5}, + "momentum + quality + qqq_trend": {"momentum": 1 / 3, "quality": 1 / 3, "qqq_trend": 1 / 3}, + "momentum .5 + quality .25 + qqq_trend .25": {"momentum": .5, "quality": .25, "qqq_trend": .25}} + for name, w in books.items(): + r = monthly_blend(R, w) + out[name] = show(name, r, b) + r = monthly_blend(R, books["momentum + quality + qqq_trend"]) + yr = pd.DataFrame({"book": (1 + r).groupby(r.index.year).prod() - 1, + "QQQ": (1 + b).groupby(b.index.year).prod() - 1}) + print("\nyearly, momentum + quality + qqq_trend vs QQQ (%):") + print((yr * 100).round(1).assign(beat=yr["book"] > yr["QQQ"]).to_string()) + print("\nodds of beating QQQ over a horizon (historical rolling windows | stationary block bootstrap, 63-day blocks)") + cands = {"momentum": R["momentum"], "qqq_trend": R["qqq_trend"], + "momentum + qqq_trend": monthly_blend(R, books["momentum + qqq_trend"]), + "momentum + quality + qqq_trend": r} + rng = np.random.default_rng(0) + for name, x in cands.items(): + cells = [] + for years in (1, 3, 5): + h = 252 * years + lx, lb = np.log1p(x.to_numpy()), np.log1p(b.reindex(x.index).fillna(0.0).to_numpy()) + cx, cb = np.r_[0, np.cumsum(lx)], np.r_[0, np.cumsum(lb)] + hist = float(((cx[h:] - cx[:-h]) > (cb[h:] - cb[:-h])).mean()) + wins = 0 + for _ in range(2000): + idx = [] + while len(idx) < h: + s = rng.integers(0, len(lx)) + idx.extend(range(s, min(s + rng.geometric(1 / 63), len(lx)))) + idx = idx[:h] + wins += lx[idx].sum() > lb[idx].sum() + cells.append(f"{years}y {hist*100:3.0f}% | {wins / 20:3.0f}%") + print(f" {name:32s} " + " ".join(cells)) + + with open(os.path.join(RESULTS, "beat_qqq.json"), "w") as fh: + json.dump(out, fh, indent=1, default=float) + R.to_pickle(os.path.join(RESULTS, "beat_qqq_sleeves.pkl")) + + +if __name__ == "__main__": + main() diff --git a/research/lab/diversifiers.py b/research/lab/diversifiers.py new file mode 100644 index 0000000..b6f057d --- /dev/null +++ b/research/lab/diversifiers.py @@ -0,0 +1,140 @@ +""" +Return sources outside S&P 500 stocks for the US book: multi-asset trend on +ETFs, and trend on bitcoin. Does either add to the three current sleeves +(momentum, low risk, TrendRider)? + + PYTHONPATH=. uv run python -m research.lab.fetch_long_history # once (BTC) + PYTHONPATH=. uv run python -m research.lab.diversifiers + +Multi-asset trend (managed futures through ETFs): each asset held while above +its n-day average, sized to 10% vol (EWMA, 30-day half-life), gross capped at +1x, the rest in cash at the T-bill rate. ETFs from data/research/tr_long.pkl +(2008-) and data/us.csv (DBC). + +Bitcoin trend: long BTC or cash, decided on closes through t-1, traded at the +close of t, 10 bp per unit traded. Tradable in a US brokerage account through +spot bitcoin ETFs (IBIT, FBTC; 2024-), before that only GBTC (OTC, premium / +discount) or futures ETFs (BITO, 2021-, roll cost). The backtest uses BTC-USD +on US trading days (Yahoo's daily close is 00:00 UTC, ~4-5 hours after the +US close), so it is a proxy for the ETF. + sma{n} hold while BTC > its n-day average + tsmom{n} hold while the n-day return is positive + vt60 the same, sized to a 60% vol target (cap 1x) +The whole grid is printed; the median rule is the honest estimate. +""" + +from __future__ import annotations + +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factor_scan import RESULTS +from research.lab.levered_book import fmt, stats, tbill +from research.lab.us_final import monthly_blend + +warnings.filterwarnings("ignore") +SRC = os.path.join("data", "research", "long_history.pkl") +COST = 0.001 +SPAN = ("2011-01-03", "2026-09-24") +BTC_BLOCKS = [("2015-06-01", "2017-12-31"), ("2018-01-01", "2020-12-31"), ("2021-01-01", "2023-12-31"), + ("2024-01-01", "2026-09-25")] +THREE = {"momentum": 1 / 3, "lowrisk": 1 / 3, "trend": 1 / 3} + + +def multi_asset_trend(close: pd.DataFrame, rf: pd.Series, assets: list[str], n: int, + target: float = 0.10) -> tuple[pd.Series, pd.Series]: + r = close.pct_change(fill_method=None) + W = {} + for a in assets: + px = close[a] + on = (px > px.rolling(n).mean()).astype(float) + vol = r[a].ewm(halflife=30, min_periods=60).std() * np.sqrt(252) + W[a] = (on * (target / vol).clip(upper=3.0)).shift(1) + W = pd.DataFrame(W).fillna(0.0) + W = W.div(W.sum(axis=1).clip(lower=1.0), axis=0) + cash = 1 - W.sum(axis=1) + ret = (W * r[assets].fillna(0.0)).sum(axis=1) + cash * rf.reindex(W.index) / 252 \ + - W.diff().abs().sum(axis=1) * COST + return ret, W.sum(axis=1) + + +def btc_trend(px: pd.Series, rf: pd.Series, rule: str, n: int, + vol_target: float | None = None) -> tuple[pd.Series, pd.Series]: + r = px.pct_change(fill_method=None).fillna(0.0) + on = px > (px.rolling(n).mean() if rule == "sma" else px.shift(n)) + expo = on.astype(float) + if vol_target: + vol = r.ewm(halflife=20, min_periods=40).std() * np.sqrt(252) + expo = expo * (vol_target / vol).clip(upper=1.0).fillna(0.0) + expo = expo.shift(1).fillna(0.0) + out = expo * r + (1 - expo) * rf.reindex(r.index) / 252 - expo.diff().abs().fillna(0.0) * COST + return out, expo + + +def with_book(book: pd.DataFrame, name: str, r: pd.Series, rf: pd.Series, start: str) -> None: + span = book.loc[start:SPAN[1], ["momentum", "lowrisk", "trend"]].fillna(0.0) + span[name] = r.reindex(span.index).fillna(0.0) + m = (1 + span).resample("ME").prod() - 1 + print(f"monthly correlation of {name} with the sleeves:", + m.corr()[name].drop(name).round(2).to_dict()) + print(fmt(f"3 sleeves {start[:4]}-", stats(monthly_blend(span, THREE), rf))) + for w in (0.1, 0.25): + mix = {k: v * (1 - w) for k, v in THREE.items()} | {name: w} + print(fmt(f"3 sleeves + {int(w*100)}% {name}", stats(monthly_blend(span, mix), rf))) + + +def main() -> None: + book = pd.read_pickle(os.path.join(RESULTS, "us_final_sleeves.pkl")) + rf = tbill(book.index.union(pd.bdate_range("2014-01-01", "2026-09-30"))) + + print("=== multi-asset trend on ETFs, 2010-03-") + with open(os.path.join("data", "research", "tr_long.pkl"), "rb") as fh: + etf = pickle.load(fh)["close"].ffill() + etf = etf.join(pd.read_csv(os.path.join("data", "us.csv"), index_col=0, parse_dates=True)[["DBC"]], how="left") + best = None + for assets in (["SPY", "TLT", "GLD"], ["TQQQ", "TLT", "GLD"], ["SPY", "TLT", "GLD", "DBC"]): + for n in (100, 200): + r, g = multi_asset_trend(etf, rf, assets, n) + r = r.loc["2010-03-01":] + s = stats(r, rf) + print(fmt(f"{'+'.join(assets)} ma{n}", s, g.loc["2010-03-01":])) + if best is None or s["sharpe"] > best[1]: + best = (r, s["sharpe"], f"{'+'.join(assets)} ma{n}") + print(f"best of the grid: {best[2]} (optimistic: picked in sample)") + with_book(book, "xtrend", best[0], rf, SPAN[0]) + + if not os.path.exists(SRC): + print("\nno data/research/long_history.pkl: run research.lab.fetch_long_history for the BTC part") + return + print("\n=== bitcoin trend, 2015-06-") + with open(SRC, "rb") as fh: + src = pickle.load(fh) + days = src["close"]["SPY"].dropna().loc["2014-12-01":].index + btc = src["close"]["BTC-USD"].reindex(days, method="ffill").dropna() + hold = btc.pct_change().fillna(0.0).loc["2015-06-01":] + print(fmt("hold BTC", stats(hold, rf))) + rows = {} + for rule in ("sma", "tsmom"): + for n in (20, 50, 100, 150, 200): + for vt in (None, 0.6): + r, e = btc_trend(btc, rf, rule, n, vt) + r = r.loc["2015-06-01":] + key = f"{rule}{n}{' vt60' if vt else ''}" + rows[key] = r + bl = " ".join(f"{stats(r.loc[a:b], rf)['sharpe']:5.2f}" for a, b in BTC_BLOCKS) + print(fmt(key, stats(r, rf), e.loc["2015-06-01":]), f"| blocks {bl}") + R = pd.DataFrame(rows) + sh = R.apply(lambda x: stats(x, rf)["sharpe"]).sort_values() + med = sh.index[len(sh) // 2] + print(f"median-Sharpe rule: {med} ({sh[med]:.2f}); range {sh.min():.2f}-{sh.max():.2f}") + with_book(book, "btc", R[med], rf, "2015-06-01") + with open(os.path.join(RESULTS, "diversifiers.pkl"), "wb") as fh: + pickle.dump({"btc": R, "btc_median_rule": med, "xtrend": best[0]}, fh) + + +if __name__ == "__main__": + main() diff --git a/research/lab/earnings.py b/research/lab/earnings.py new file mode 100644 index 0000000..6bd7929 --- /dev/null +++ b/research/lab/earnings.py @@ -0,0 +1,124 @@ +""" +Earnings-announcement return (EAR) factor: the market's reaction to the +latest earnings release, which keeps drifting afterwards (Chan, Jegadeesh & +Lakonishok 1996; Brandt, Kishore, Santa-Clara & Venkatachalam 2008 find the +drift after the announcement *return* stronger and longer-lived than after +the earnings surprise). + + PYTHONPATH=. uv run --with numba python -m research.lab.earnings + +Announcement dates: 8-K filings with Item 2.02 (Results of Operations), +data/research/sec_8k_earnings.csv (written by research/lab/fetch_lazy_prices.py +from the SEC submissions JSON). A release after the close on day d is +reacted to on d+1, so the event window is the three sessions d-1..d+1: + + ear = sum of (stock return - equal-weight member return) over d-1, d, d+1 + +It is known at the close of d+1, the score is dated then (the harness trades +at the next open) and held until the next release, at most 100 trading days. +Also: ``ear_z`` = ear divided by the stock's 60-day residual vol (a surprise in +units of normal noise), and ``ear_count`` = the sum of the last four signs. +""" + +from __future__ import annotations + +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factors import _mkt, _ret + +warnings.filterwarnings("ignore") +SRC = os.path.join("data", "research", "sec_8k_earnings.csv") +CACHE = os.path.join("data", "research", "earnings_factors.pkl") +MAX_AGE = 100 + + +def announcement_days(p) -> pd.DataFrame: + """Bool panel: an Item 2.02 8-K was filed that day (mapped to the next + trading day when filed on a holiday or weekend).""" + ev = pd.read_csv(SRC, parse_dates=["filingDate"]) + ev = ev[ev["ticker"].isin(p.close.columns)] + pos = p.dates.searchsorted(ev["filingDate"].values) + ev = ev.assign(pos=pos)[pos < len(p.dates)] + A = np.zeros(p.close.shape, bool) + col = {c: i for i, c in enumerate(p.close.columns)} + A[ev["pos"].to_numpy(), ev["ticker"].map(col).to_numpy()] = True + return pd.DataFrame(A, index=p.dates, columns=p.close.columns) + + +def _hold(values: pd.DataFrame, known: pd.DataFrame, max_age: int) -> pd.DataFrame: + """Carry each event's value forward from the day it is known, for at most + ``max_age`` trading days.""" + idx = np.broadcast_to(np.arange(len(values))[:, None], values.shape) + last_pos = pd.DataFrame(np.where(known.values, idx, np.nan), index=values.index, + columns=values.columns).ffill() + return values.where(known).ffill().mask(idx - last_pos.values > max_age) + + +def build(p) -> dict[str, pd.DataFrame]: + A = announcement_days(p) + ab = _ret(p).sub(_mkt(p), axis=0) # abnormal daily return + car = ab.shift(1) + ab + ab.shift(-1) # d-1..d+1, indexed at d + known = A.shift(1, fill_value=False) # known at the close of d+1 + vol = ab.rolling(60, min_periods=40).std().shift(2) # pre-event noise + ear = _hold(car.shift(1), known, MAX_AGE) + ear_z = _hold((car / vol).shift(1), known, MAX_AGE) + sign = np.sign(car.shift(1)).where(known) + count = pd.DataFrame(np.nan, index=p.dates, columns=p.close.columns) + for c in count.columns: + x = sign[c].dropna() + count.loc[x.index, c] = x.rolling(4, min_periods=2).sum() + ear_count = _hold(count, known & count.notna(), MAX_AGE) + return {"ear": ear, "ear_z": ear_z, "ear_count": ear_count, "announce": A} + + +def load(p, rebuild: bool = False) -> dict[str, pd.DataFrame]: + if os.path.exists(CACHE) and not rebuild: + with open(CACHE, "rb") as fh: + return pickle.load(fh) + out = build(p) + with open(CACHE, "wb") as fh: + pickle.dump(out, fh, protocol=4) + return out + + +def main() -> None: + from research.lab.alpha import alpha_stats, bench_returns, fmt_alpha + from research.lab.factor_scan import load as load_panel + from research.lab.report import forward_returns + from research.lab.rotation import factor_ranks + from research.lab.us_sweep import _blend, run + + p, _ = load_panel("us") + E = load(p, rebuild=True) + A = E["announce"] & p.member + print("announcements among members per year:", A.sum(axis=1).groupby(p.dates.year).sum().to_dict()) + print("members with a live ear value (median per day by year):", + E["ear"].where(p.member).notna().sum(axis=1).groupby(p.dates.year).median().to_dict()) + fwd = forward_returns(p, 21).where(p.member) + bench = bench_returns(p.dates) + pr = factor_ranks(p, "us") + mom = _blend(pr["overnight_mom"], pr["rec_126"]) + for name in ("ear", "ear_z", "ear_count"): + rk = E[name].where(p.member).rank(axis=1, pct=True) + ic = rk.corrwith(fwd.rank(axis=1, pct=True), axis=1).iloc[::21].dropna() + for a, b in (("2006", "2016"), ("2017", "2026")): + x = ic.loc[a:b] + print(f"{name} IC21 {a}-{b}: {x.mean():+.4f} t {x.mean() / x.std() * np.sqrt(len(x)):+.2f}") + for label, sc, kw in ((f"{name} top20 e21", rk, dict(top_n=20, every=21)), + (f"{name} top20 e10", rk, dict(top_n=20, every=10)), + (f"momentum + {name} top10 e10", _blend(mom, rk), dict(top_n=10, every=10))): + st = run(p, sc.where(p.member), [("2006-09-27", "2026-09-24")], buffer=2.0, **kw) + r = st["_rets"] + print(fmt_alpha(label, alpha_stats(r, bench)), "| holdout", + fmt_alpha("", alpha_stats(r.loc["2017":], bench)).split("alpha", 1)[1][:40], flush=True) + st = run(p, mom.where(p.member), [("2006-09-27", "2026-09-24")], top_n=10, buffer=2.0, every=10) + print(fmt_alpha("momentum alone top10 e10", alpha_stats(st["_rets"], bench))) + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_insider.py b/research/lab/fetch_insider.py new file mode 100644 index 0000000..50016f1 --- /dev/null +++ b/research/lab/fetch_insider.py @@ -0,0 +1,372 @@ +""" +Fetch SEC "Insider Transactions Data Sets" (structured Form 3/4/5 data, +quarterly zips from 2006Q1) and keep the open-market trades of officers and +directors of S&P 500 ever-members in the research price panel. + + uv run python -m research.lab.fetch_insider + +Source: https://www.sec.gov/data-research/sec-markets-data/insider-transactions-data-sets +The zip links are read from that page (the 2026Q2 file moved to a new path). +Zips are cached in data/research/ext/insider/ and never downloaded twice. +Every request sends the SEC User-Agent, goes out directly (no proxy) and is +spaced >= 0.25 s from the previous one (SEC allows 10/s; we stay under 5/s). + +Kept rows (NONDERIV_TRANS joined to SUBMISSION and REPORTINGOWNER): + * DOCUMENT_TYPE "4". Amendments (4/A, ~3% of Form 4s) are dropped: they + mostly restate the original's trades, so keeping both double counts, and + the original is what the market saw on its filing date. + * TRANS_CODE P (open-market purchase, acquired) or S (open-market sale, + disposed) of a non-preferred equity security, no equity swap. + * The reporting owner is an officer and/or director. Pure 10% owners and + "other" are dropped. A joint filing (rare) is attributed to one officer/ + director owner (lowest CIK) so a trade is never counted twice. + * Exact duplicate lines within a filing (same date, code, shares, price, + shares owned after, ownership form) are dropped. + +Issuer -> panel ticker: the current SEC ticker map (company_tickers.json) by +issuer CIK; share classes of one CIK (GOOG/GOOGL, ...) all map to it. Issuers +missing from it (renamed, reorganized, taken private) fall back to the +filing's own ISSUERTRADINGSYMBOL, limited to the dates that symbol was used +and to dates before the ticker's current owner first filed, and kept only if +trade prices match the ticker's unadjusted closes (median error < 25%). + +Writes data/research/insider_trades.pkl: + {"trades": DataFrame, one row per kept transaction: ticker, cik, + insider_cik, role, title, ten_pct, filed, trans_date, code, + shares, price, value, price_ok, plan_10b51, plan_fn, ownership, + accession + "tickers": DataFrame[ticker, cik, source, first, last, n_filings, price_err] + (which issuer CIK covers each panel ticker, and when), + "unmapped": panel tickers with no issuer} + +Point in time: a trade is usable from the first trading day strictly after +``filed`` (FILING_DATE). ``trans_date`` is only used to classify traders +(calendar month and year), never to decide when a trade is known. +``value`` = shares x reported price; where the price is missing or more than +2x off the issuer's usual price / unadjusted-close ratio (``price_ok`` False: +mostly total proceeds typed as the price), shares x the unadjusted close. +``plan_10b51`` is the Rule 10b5-1 checkbox (Form 4 from April 2023, NaN +before); ``plan_fn`` is True when a footnote or remark of the filing mentions +Rule 10b5-1 (available in every year). +""" + +from __future__ import annotations + +import csv +import os +import pickle +import re +import time +import zipfile +from urllib.parse import urljoin + +import numpy as np +import pandas as pd +import requests + +from research.lab.fetch_us_fundamentals import UA, cik_map + +PAGE = "https://www.sec.gov/data-research/sec-markets-data/insider-transactions-data-sets" +ZIP_DIR = os.path.join("data", "research", "ext", "insider") +OUT = os.path.join("data", "research", "insider_trades.pkl") +FUND_RAW = os.path.join("data", "research", "us_fund_raw.pkl") + +SUB_COLS = ["ACCESSION_NUMBER", "FILING_DATE", "DOCUMENT_TYPE", "ISSUERCIK", + "ISSUERTRADINGSYMBOL", "REMARKS", "AFF10B5ONE"] +OWN_COLS = ["ACCESSION_NUMBER", "RPTOWNERCIK", "RPTOWNER_RELATIONSHIP", "RPTOWNER_TITLE"] +TRN_COLS = ["ACCESSION_NUMBER", "SECURITY_TITLE", "TRANS_DATE", "TRANS_CODE", "EQUITY_SWAP_INVOLVED", + "TRANS_SHARES", "TRANS_PRICEPERSHARE", "TRANS_ACQUIRED_DISP_CD", + "SHRS_OWND_FOLWNG_TRANS", "DIRECT_INDIRECT_OWNERSHIP"] +PLAN_RE = r"10b5\s*-\s*1|10b-5-1|10b5one" +NOT_COMMON = r"prefer|pfd|note|debenture|warrant|bond|option" + +_last = [0.0] + + +def _session() -> requests.Session: + s = requests.Session() + s.trust_env = False # direct connection: ignore proxy env vars + s.headers.update(UA) + return s + + +def _get(sess: requests.Session, url: str, **kw) -> requests.Response: + for attempt in range(5): + wait = 0.25 - (time.time() - _last[0]) + if wait > 0: + time.sleep(wait) + _last[0] = time.time() + try: + r = sess.get(url, timeout=300, **kw) + r.raise_for_status() + return r + except Exception as exc: # noqa: BLE001 — retry transient SEC errors + print(f" {url}: {exc}; retry {attempt + 1}", flush=True) + time.sleep(5 * (attempt + 1)) + raise RuntimeError(f"failed: {url}") + + +def zip_urls(sess: requests.Session) -> dict[str, str]: + html = _get(sess, PAGE).text + found = re.findall(r'href="([^"]*?(\d{4}q[1-4])_form345\.zip)"', html) + return {q: urljoin(PAGE, href) for href, q in sorted(found, key=lambda x: x[1])} + + +def download(sess: requests.Session) -> list[str]: + os.makedirs(ZIP_DIR, exist_ok=True) + urls = zip_urls(sess) + print(f"{len(urls)} quarterly zips listed ({min(urls)} .. {max(urls)})", flush=True) + paths = [] + for q, url in urls.items(): + path = os.path.join(ZIP_DIR, f"{q}_form345.zip") + if not (os.path.exists(path) and zipfile.is_zipfile(path)): + r = _get(sess, url) + with open(path + ".part", "wb") as fh: + fh.write(r.content) + os.replace(path + ".part", path) + print(f" {q}: {len(r.content) / 1e6:.1f} MB", flush=True) + paths.append(path) + return paths + + +# --------------------------------------------------------------------------- +# parse one quarter +# --------------------------------------------------------------------------- + +def _read(z: zipfile.ZipFile, name: str, cols: list[str]) -> pd.DataFrame: + with z.open(name) as fh: + return pd.read_csv(fh, sep="\t", dtype=str, quoting=csv.QUOTE_NONE, + usecols=lambda c: c in cols, encoding_errors="replace") + + +def _date(s: pd.Series) -> pd.Series: + return pd.to_datetime(s, format="%d-%b-%Y", errors="coerce") + + +def parse_quarter(path: str) -> tuple[pd.DataFrame, pd.DataFrame]: + """(kept trades, every Form 4 filing's issuer CIK / symbol / date).""" + z = zipfile.ZipFile(path) + sub = _read(z, "SUBMISSION.tsv", SUB_COLS) + sub = sub[sub["DOCUMENT_TYPE"].str.strip() == "4"].copy() + sub["filed"] = _date(sub["FILING_DATE"]) + sub["cik"] = pd.to_numeric(sub["ISSUERCIK"], errors="coerce") + sub = sub.dropna(subset=["filed", "cik"]) + filings = sub[["ACCESSION_NUMBER", "cik", "ISSUERTRADINGSYMBOL", "filed"]] + + own = _read(z, "REPORTINGOWNER.tsv", OWN_COLS) + own = own[own["ACCESSION_NUMBER"].isin(sub["ACCESSION_NUMBER"])] + rel = own["RPTOWNER_RELATIONSHIP"].fillna("").str.lower() + own = own.assign(officer=rel.str.contains("officer"), director=rel.str.contains("director"), + ten_pct=rel.str.contains("tenpercent"), + insider_cik=pd.to_numeric(own["RPTOWNERCIK"], errors="coerce")) + own = own[(own["officer"] | own["director"]) & own["insider_cik"].notna()] + own = own.sort_values(["ACCESSION_NUMBER", "insider_cik"]).drop_duplicates("ACCESSION_NUMBER") + own["role"] = np.where(own["officer"] & own["director"], "officer+director", + np.where(own["officer"], "officer", "director")) + + tr = _read(z, "NONDERIV_TRANS.tsv", TRN_COLS) + code = tr["TRANS_CODE"].str.strip() + ad = tr["TRANS_ACQUIRED_DISP_CD"].str.strip() + keep = (((code == "P") & (ad == "A")) | ((code == "S") & (ad == "D"))) + keep &= tr["ACCESSION_NUMBER"].isin(own["ACCESSION_NUMBER"]) + keep &= ~tr["SECURITY_TITLE"].fillna("").str.contains(NOT_COMMON, case=False, regex=True) + keep &= ~tr["EQUITY_SWAP_INVOLVED"].fillna("0").str.strip().isin(["1", "true"]) + tr = tr[keep].copy() + tr["code"] = tr["TRANS_CODE"].str.strip() + tr["trans_date"] = _date(tr["TRANS_DATE"]) + for src, dst in (("TRANS_SHARES", "shares"), ("TRANS_PRICEPERSHARE", "price"), + ("SHRS_OWND_FOLWNG_TRANS", "owned_after")): + tr[dst] = pd.to_numeric(tr[src], errors="coerce") + tr["ownership"] = tr["DIRECT_INDIRECT_OWNERSHIP"].str.strip() + tr = tr[tr["shares"] > 0] + tr = tr.drop_duplicates(["ACCESSION_NUMBER", "trans_date", "code", "shares", "price", + "owned_after", "ownership"]) + + # Rule 10b5-1: the 2023 checkbox, and any footnote / remark that cites the rule + fn = _read(z, "FOOTNOTES.tsv", ["ACCESSION_NUMBER", "FOOTNOTE_TXT"]) + fn = fn[fn["ACCESSION_NUMBER"].isin(tr["ACCESSION_NUMBER"])] + cited = set(fn.loc[fn["FOOTNOTE_TXT"].str.contains(PLAN_RE, case=False, na=False, regex=True), + "ACCESSION_NUMBER"]) + cited |= set(sub.loc[sub["REMARKS"].str.contains(PLAN_RE, case=False, na=False, regex=True), + "ACCESSION_NUMBER"]) + box = (sub["AFF10B5ONE"].str.strip().str.lower().map({"1": 1.0, "true": 1.0, "0": 0.0, "false": 0.0}) + if "AFF10B5ONE" in sub else pd.Series(np.nan, index=sub.index)) + sub = sub.assign(plan_10b51=box) + + out = (tr.merge(own[["ACCESSION_NUMBER", "insider_cik", "role", "RPTOWNER_TITLE", "ten_pct"]], + on="ACCESSION_NUMBER") + .merge(sub[["ACCESSION_NUMBER", "cik", "ISSUERTRADINGSYMBOL", "filed", "plan_10b51"]], + on="ACCESSION_NUMBER")) + out["plan_fn"] = out["ACCESSION_NUMBER"].isin(cited) + out = out.rename(columns={"ACCESSION_NUMBER": "accession", "RPTOWNER_TITLE": "title", + "ISSUERTRADINGSYMBOL": "symbol"}) + return (out[["accession", "cik", "symbol", "insider_cik", "role", "title", "ten_pct", "filed", + "trans_date", "code", "shares", "price", "plan_10b51", "plan_fn", "ownership"]], + filings) + + +# --------------------------------------------------------------------------- +# issuer CIK -> panel ticker +# --------------------------------------------------------------------------- + +def norm_symbols(sym) -> list[str]: + """'(WST)' -> [WST]; 'NYSE: KRC' -> [KRC]; 'CBS, CBS.A' -> [CBS, CBS-A].""" + if not isinstance(sym, str): + return [] + s = re.sub(r'["()\[\]]', "", sym.upper()) + s = re.sub(r"\b(NYSE|NASDAQ|NYSE AMERICAN|AMEX|NMS|OTC)\s*:", "", s) + parts = re.split(r"[,;/&]| AND ", s) + return [re.sub(r"\s+", "", x).replace(".", "-") for x in parts if x.strip()] + + +def raw_closes(tickers: list[str]) -> pd.DataFrame | None: + """Unadjusted closes (split-adjusted close x later split ratios), from the + fundamentals fetch; used to check that an issuer's trade prices belong to + the ticker, and to value trades whose reported price is off.""" + if not os.path.exists(FUND_RAW): + return None + with open(FUND_RAW, "rb") as fh: + raw = pickle.load(fh) + close_sa = raw["close_sa"].reindex(columns=tickers) + ratio = raw["splits"].reindex(index=close_sa.index, columns=tickers).fillna(0.0).replace(0.0, 1.0) + after = ratio.iloc[::-1].cumprod().iloc[::-1].shift(-1).fillna(1.0) + return close_sa * after + + +def _raw_at(raw: pd.DataFrame, tickers: pd.Series, dates: pd.Series) -> np.ndarray: + """raw close of each (ticker, date) pair as of that date (last close on or before).""" + out = np.full(len(tickers), np.nan) + ff = raw.ffill() + col = {t: i for i, t in enumerate(ff.columns)} + ci = tickers.map(col) + ri = ff.index.searchsorted(dates.values, side="right") - 1 + ok = ci.notna().values & (ri >= 0) & dates.notna().values + out[ok] = ff.values[ri[ok], ci[ok].astype(int)] + return out + + +def map_tickers(trades: pd.DataFrame, filings: pd.DataFrame, panel_tickers: list[str], + raw: pd.DataFrame | None, close: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + """Assign each trade to a panel ticker; returns (trades with ticker, ticker map).""" + cmap = cik_map() + pset = set(panel_tickers) + rev: dict[int, list[str]] = {} + for t, c in cmap.items(): + if t in pset: + rev.setdefault(c, []).append(t) + span = filings.groupby("cik")["filed"].agg(["min", "max", "count"]) + rows = [] + for c, ts in rev.items(): + if c in span.index: + for t in sorted(ts): + rows.append((t, c, "sec_map", span.at[c, "min"], span.at[c, "max"], span.at[c, "count"])) + # symbol fallback for issuers the current SEC map does not tie to a panel ticker + fb = filings[~filings["cik"].isin(list(rev))] + pairs = (fb.groupby(["cik", "ISSUERTRADINGSYMBOL"])["filed"] + .agg(first="min", last="max", n="count").reset_index()) + cand = [] + for r in pairs.itertuples(index=False): + for t in norm_symbols(r.ISSUERTRADINGSYMBOL): + if t in pset: + cand.append((t, int(r.cik), r.first, r.last, r.n)) + cand = (pd.DataFrame(cand, columns=["ticker", "cik", "first", "last", "n"]) + .groupby(["ticker", "cik"]).agg(first=("first", "min"), last=("last", "max"), n=("n", "sum")) + .reset_index()) + px_first = close.apply(lambda s: s.first_valid_index()) + px_last = close.apply(lambda s: s.last_valid_index()) + for r in cand.sort_values("n", ascending=False).itertuples(index=False): + first, last = r.first, r.last + owner = cmap.get(r.ticker) + if owner is not None and owner in span.index: # a predecessor only, never alongside the owner + last = min(last, span.at[owner, "min"] - pd.Timedelta(days=1)) + if r.ticker in px_first.index and pd.notna(px_first[r.ticker]): + first, last = max(first, px_first[r.ticker]), min(last, px_last[r.ticker]) + taken = [(a, b) for (t, _, _, a, b, _) in rows if t == r.ticker] + if last <= first or any(a <= last and first <= b for a, b in taken): + continue + rows.append((r.ticker, r.cik, "symbol", first, last, r.n)) + tm = pd.DataFrame(rows, columns=["ticker", "cik", "source", "first", "last", "n_filings"]) + + # assign trades: SEC-map issuers to the filing's own share class when it is + # one of theirs, fallback issuers only inside their window + prim = {c: ts[0] for c, ts in ((c, sorted(ts)) for c, ts in rev.items())} + t_sec = trades["cik"].map(prim) + sym0 = trades["symbol"].map(lambda s: (norm_symbols(s) or [None])[0]) + own_class = [s if (isinstance(s, str) and c in rev and s in rev[c]) else None + for s, c in zip(sym0.values, trades["cik"].values)] + trades = trades.assign(ticker=pd.Series(own_class, index=trades.index).fillna(t_sec)) + for r in tm[tm["source"] == "symbol"].itertuples(index=False): + m = (trades["cik"] == r.cik) & trades["filed"].between(r.first, r.last) & trades["ticker"].isna() + trades.loc[m, "ticker"] = r.ticker + trades = trades.dropna(subset=["ticker"]).copy() + + # price check: trade prices vs the ticker's unadjusted close on the trade date + trades["raw_close"] = _raw_at(raw, trades["ticker"], trades["trans_date"]) if raw is not None else np.nan + with np.errstate(divide="ignore", invalid="ignore"): + err = np.abs(np.log(trades["price"] / trades["raw_close"])) + trades["price_err"] = err.where(trades["price"] > 0) + chk = trades.groupby(["ticker", "cik"])["price_err"].median().rename("price_err") + tm = tm.merge(chk, left_on=["ticker", "cik"], right_index=True, how="left") + bad = tm[(tm["source"] == "symbol") & (tm["price_err"] > np.log(1.25))] + if len(bad): + print(" symbol matches rejected on price:", bad[["ticker", "cik", "price_err"]].to_dict("records")) + tm = tm.drop(bad.index) + keep = pd.MultiIndex.from_frame(trades[["ticker", "cik"]]).isin( + pd.MultiIndex.from_frame(bad[["ticker", "cik"]])) + trades = trades[~keep] + return trades, tm.reset_index(drop=True) + + +def main() -> None: + for k in ("http_proxy", "https_proxy", "HTTP_PROXY", "HTTPS_PROXY", "all_proxy", "ALL_PROXY"): + os.environ.pop(k, None) # direct connection for cik_map as well + sess = _session() + paths = download(sess) + parts, fil, t0 = [], [], time.time() + for path in paths: + tr, f = parse_quarter(path) + parts.append(tr) + fil.append(f) + print(f" {os.path.basename(path)}: {len(tr):,} P/S lines by officers/directors " + f"[{time.time() - t0:.0f}s]", flush=True) + trades = pd.concat(parts, ignore_index=True) + filings = pd.concat(fil, ignore_index=True).drop_duplicates("ACCESSION_NUMBER") + trades = trades.drop_duplicates(["accession", "trans_date", "code", "shares", "price", "ownership"]) + trades["cik"] = trades["cik"].astype("int64") + filings["cik"] = filings["cik"].astype("int64") + print(f"{len(trades):,} kept lines, {len(filings):,} Form 4 filings", flush=True) + + with open(os.path.join("data", "research", "us_ohlcv.pkl"), "rb") as fh: + close = pickle.load(fh)["close"] + tickers = list(close.columns) + raw = raw_closes(tickers) + trades, tm = map_tickers(trades, filings, tickers, raw, close) + + # value at the reported price unless it is > 2x off the issuer's usual + # price / close ratio (~1 normally, ~1500 for BRK-A trades mapped to BRK-B): + # such prices are unit errors (total proceeds typed as the price) and are + # replaced by the ticker's unadjusted close + with np.errstate(divide="ignore", invalid="ignore"): + lr = pd.Series(np.log(trades["price"] / trades["raw_close"]), index=trades.index) + lr = lr.where(trades["price"] > 0) + usual = lr.groupby([trades["ticker"], trades["cik"]]).transform("median") + trades["price_ok"] = (lr - usual).abs() < np.log(2) + px = np.where(trades["price_ok"], trades["price"], + np.where(trades["raw_close"] > 0, trades["raw_close"], trades["price"])) + trades["value"] = trades["shares"] * np.where(px > 0, px, np.nan) + trades["insider_cik"] = trades["insider_cik"].astype("int64") + cols = ["ticker", "cik", "insider_cik", "role", "title", "ten_pct", "filed", "trans_date", "code", + "shares", "price", "value", "price_ok", "plan_10b51", "plan_fn", "ownership", "accession"] + trades = trades[cols].sort_values(["filed", "ticker"]).reset_index(drop=True) + unmapped = sorted(set(tickers) - set(tm["ticker"])) + print(f"{len(trades):,} trades in {trades['ticker'].nunique()} panel tickers " + f"({(trades['code'] == 'P').sum():,} P / {(trades['code'] == 'S').sum():,} S); " + f"map: {(tm['source'] == 'sec_map').sum()} by CIK, {(tm['source'] == 'symbol').sum()} by symbol; " + f"unmapped {unmapped}", flush=True) + with open(OUT, "wb") as fh: + pickle.dump({"trades": trades, "tickers": tm, "unmapped": unmapped}, fh, protocol=4) + print(f"Saved {OUT}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_lazy_prices.py b/research/lab/fetch_lazy_prices.py new file mode 100644 index 0000000..a1d4150 --- /dev/null +++ b/research/lab/fetch_lazy_prices.py @@ -0,0 +1,527 @@ +""" +"Lazy Prices" text similarity (Cohen, Malloy & Nguyen 2020, JF) for every +S&P 500 ever-member in the US research panel, from SEC EDGAR. + + PYTHONPATH=. uv run python -m research.lab.fetch_lazy_prices # index, 10-K, 10-Q + PYTHONPATH=. uv run python -m research.lab.fetch_lazy_prices --forms 10-K # 10-Ks only + nohup env PYTHONPATH=. uv run python -m research.lab.fetch_lazy_prices \ + > logs/fetch_lazy_prices.log 2>&1 & + +Stages, each resumable (a rerun skips what is already on disk): + 1. index submissions JSON of each CIK (+ its older pages from 2004 on) → + data/research/lazy_prices_index.csv (original 10-K/10-K405/10-Q), + data/research/sec_sic.csv (ticker, cik, sic, sicDescription, name), + data/research/sec_8k_earnings.csv (8-K and 8-K/A with item 2.02). + 2. 10-K each 10-K vs the previous 10-K (period end 300-430 days earlier). + 3. 10-Q each 10-Q vs the 10-Q of the same fiscal quarter a year earlier + (period end 350-380 days earlier). + +Text: the filing's primary document only (no exhibits); HTML, hidden inline +XBRL, scripts and tables whose characters are >15% digits are dropped, then +lines that are >25% digits and "table of contents" headers; lowercase words +of 2+ letters minus English stop words. Similarities: cosine of term-count +vectors and Jaccard of word sets, for the whole document and (10-K only) for +Item 1A Risk Factors when a section of 200+ words can be located. + +Output: data/research/lazy_prices.csv, one row per (ticker, filing), appended +as filings finish: ticker, cik, form, filed, period, accession, ref_accession, +cos_sim, jac_sim, n_words, cos_rf, jac_rf, n_words_rf. No documents are kept; +data/research/lazy_prices_vec/ holds the term-count vectors of the filings a +later filing still needs as its reference (plus the latest 10-K / four 10-Qs). + +SEC fair access: one User-Agent, at most 5 requests/second across all worker +processes (a shared limiter), direct connection (proxy variables ignored). +Amendments (10-K/A, 10-Q/A) are left out as in the paper: they are often +partial documents and would read as "changers". +""" + +from __future__ import annotations + +import argparse +import csv +import fcntl +import html +import multiprocessing as mp +import os +import pickle +import re +import time +import traceback +import zlib +from collections import Counter +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +import pandas as pd +import requests + +from research.lab.fetch_us_fundamentals import cik_map +from research.lab.panel import US_ETFS + +DIR = os.path.join("data", "research") +INDEX = os.path.join(DIR, "lazy_prices_index.csv") +OUT = os.path.join(DIR, "lazy_prices.csv") +SIC = os.path.join(DIR, "sec_sic.csv") +EARN = os.path.join(DIR, "sec_8k_earnings.csv") +VEC = os.path.join(DIR, "lazy_prices_vec") +COLS = ["ticker", "cik", "form", "filed", "period", "accession", "ref_accession", + "cos_sim", "jac_sim", "n_words", "cos_rf", "jac_rf", "n_words_rf"] + +UA = {"User-Agent": "quant-research/0.1 gahow@example.com", "Accept-Encoding": "gzip, deflate"} +MIN_INTERVAL = 0.21 # seconds between request starts, all processes (≈4.8/s) +START = "2004-01-01" +GROUPS = {"10-K": ("10-K", "10-K405"), "10-Q": ("10-Q",)} +# ever-members missing from today's company_tickers.json; CIKs checked by name +EXTRA_CIKS = {"AVB": 915912, "EQR": 906107, "EA": 712515, "GENZ": 732485, "GR": 42542} + +STOP = frozenset(""" +a about above after again against all am an and any are as at be because been before being +below between both but by can cannot could did do does doing down during each few for from +further had has have having he her here hers herself him himself his how i if in into is it +its itself let me more most my myself no nor not of off on once only or other ought our ours +ourselves out over own same she should so some such than that the their theirs them themselves +then there these they this those through to too under until up very was we were what when where +which while who whom why with would you your yours yourself yourselves will shall +""".split()) + +# --------------------------------------------------------------------------- +# HTTP: a limiter shared by threads and worker processes +# --------------------------------------------------------------------------- + +_LOCK = None +_NEXT = None +_SESSION: requests.Session | None = None + + +def _init_worker(lock, nxt) -> None: + global _LOCK, _NEXT, _SESSION + _LOCK, _NEXT = lock, nxt + _SESSION = requests.Session() + _SESSION.trust_env = False # direct connection, never the shell's proxy + + +def _throttle() -> None: + with _LOCK: + t = max(time.time(), _NEXT.value) + _NEXT.value = t + MIN_INTERVAL + wait = t - time.time() + if wait > 0: + time.sleep(wait) + + +def _get(url: str) -> requests.Response | None: + """GET with the shared limiter; None for a missing document.""" + for attempt in range(7): + _throttle() + try: + r = _SESSION.get(url, headers=UA, timeout=120) + except requests.RequestException as exc: + print(f" ! {url}: {exc.__class__.__name__}; retry {attempt + 1}", flush=True) + time.sleep(5 * (attempt + 1)) + continue + if r.status_code == 200: + return r + if r.status_code == 404: + return None + rate = r.status_code == 429 or (r.status_code == 403 + and b"Request Rate Threshold" in r.content[:5000]) + if r.status_code == 403 and not rate: + return None + # rate threshold (SEC blocks for ~10 minutes) or server error: back off + wait = 60 * (attempt + 1) if rate else 10 * (attempt + 1) + print(f" ! {url}: HTTP {r.status_code}; sleeping {wait}s", flush=True) + time.sleep(wait) + return None + + +# --------------------------------------------------------------------------- +# stage 1: filing index, SIC codes, 8-K earnings dates +# --------------------------------------------------------------------------- + +def _columns(block: dict) -> pd.DataFrame: + keys = ["accessionNumber", "filingDate", "reportDate", "acceptanceDateTime", "form", + "items", "primaryDocument"] + return pd.DataFrame({k: block.get(k, [""] * len(block["accessionNumber"])) for k in keys}) + + +def _submissions(cik: int) -> tuple[dict, pd.DataFrame] | None: + r = _get(f"https://data.sec.gov/submissions/CIK{cik:010d}.json") + if r is None: + return None + js = r.json() + parts = [_columns(js["filings"]["recent"])] + for f in js["filings"].get("files", []): + if f.get("filingTo", "9999") < START: + continue + r2 = _get(f"https://data.sec.gov/submissions/{f['name']}") + if r2 is not None: + parts.append(_columns(r2.json())) + df = pd.concat(parts, ignore_index=True).drop_duplicates("accessionNumber") + meta = {"sic": js.get("sic", ""), "sicDescription": js.get("sicDescription", ""), + "name": js.get("name", "")} + return meta, df + + +def ticker_ciks(tickers: list[str]) -> dict[str, int]: + cm = cik_map() + out = {t: cm[t] for t in tickers if t in cm} + for t, c in EXTRA_CIKS.items(): + if t in tickers and t not in out: + out[t] = c + return out + + +def build_index(tickers: list[str]) -> None: + tc = ticker_ciks(tickers) + print(f"{len(tc)}/{len(tickers)} panel tickers map to a CIK " + f"({len(set(tc.values()))} CIKs; {len(tc) - sum(t in EXTRA_CIKS for t in tc)} " + f"from company_tickers.json)", flush=True) + by_cik: dict[int, list[str]] = {} + for t, c in tc.items(): + by_cik.setdefault(c, []).append(t) + res: dict[int, tuple] = {} + t0 = time.time() + with ThreadPoolExecutor(max_workers=5) as ex: + for k, (cik, out) in enumerate(zip(by_cik, ex.map(_submissions, by_cik)), 1): + if out is not None: + res[cik] = out + if k % 50 == 0: + print(f" submissions {k}/{len(by_cik)} [{time.time() - t0:.0f}s]", flush=True) + sic, idx, earn = [], [], [] + for cik, (meta, df) in res.items(): + for t in sorted(by_cik[cik]): + sic.append({"ticker": t, "cik": cik, **meta}) + e = df[df["form"].isin(["8-K", "8-K/A"]) + & df["items"].fillna("").str.split(",").apply(lambda xs: "2.02" in [x.strip() for x in xs])] + if len(e): + earn.append(e.assign(ticker=t, cik=cik)[ + ["ticker", "cik", "form", "filingDate", "acceptanceDateTime", "accessionNumber", "items"]]) + f = df[df["form"].isin(GROUPS["10-K"] + GROUPS["10-Q"]) & (df["filingDate"] >= START)] + idx.append(pd.DataFrame({"cik": cik, "tickers": "|".join(sorted(by_cik[cik])), "form": f["form"], + "filed": f["filingDate"], "period": f["reportDate"], + "accession": f["accessionNumber"], "doc": f["primaryDocument"]})) + pd.DataFrame(sic).sort_values("ticker").to_csv(SIC, index=False) + ear = pd.concat(earn, ignore_index=True).sort_values(["ticker", "filingDate"]) + ear.to_csv(EARN, index=False) + ix = pd.concat(idx, ignore_index=True).sort_values(["cik", "filed"]) + ix.to_csv(INDEX, index=False) + print(f"index: {len(ix)} filings ({(ix['form'] != '10-Q').sum()} 10-K) for {ix['cik'].nunique()} CIKs; " + f"{len(ear)} 8-K item 2.02 rows ({ear['filingDate'].min()} → {ear['filingDate'].max()}); " + f"sic for {len(sic)} tickers [{time.time() - t0:.0f}s]", flush=True) + + +def with_refs(ix: pd.DataFrame, group: str) -> pd.DataFrame: + """Filings of one form group with their comparable prior filing.""" + f = ix[ix["form"].isin(GROUPS[group])].copy() + f["filed_d"] = pd.to_datetime(f["filed"]) + f["period_d"] = pd.to_datetime(f["period"], errors="coerce") + # a period filed twice as an original (re-filing) keeps the first copy + f = f.sort_values("filed_d") + f = pd.concat([f[f["period_d"].isna()], + f[f["period_d"].notna()].drop_duplicates(["cik", "period_d"], keep="first")]) + lo, hi = (300, 430) if group == "10-K" else (350, 380) + refs = [] + for _, g in f.sort_values(["cik", "filed_d"]).groupby("cik", sort=False): + pe, fe, acc = g["period_d"].values, g["filed_d"].values, g["accession"].values + for i in range(len(g)): + best, best_gap = "", None + for j in range(i - 1, -1, -1): + if not pd.isna(pe[i]) and not pd.isna(pe[j]): + gap = (pe[i] - pe[j]).astype("timedelta64[D]").astype(int) + ok = lo <= gap <= hi + else: # no period on one side: filing dates, wider window + gap = (fe[i] - fe[j]).astype("timedelta64[D]").astype(int) + ok = lo - 20 <= gap <= hi + 20 + if ok and (best_gap is None or abs(gap - 365) < abs(best_gap - 365)): + best, best_gap = acc[j], gap + if gap > hi + 60: + break + refs.append((g.index[i], best)) + f.loc[[i for i, _ in refs], "ref"] = [r for _, r in refs] + return f.sort_values(["cik", "filed_d"]).drop(columns=["filed_d", "period_d"]) + + +# --------------------------------------------------------------------------- +# text → term counts +# --------------------------------------------------------------------------- + +_BLOCK = ("p", "div", "br", "tr", "li", "h1", "h2", "h3", "h4", "h5", "h6", "table", "pre", + "center", "ul", "ol", "dl", "dt", "dd", "blockquote", "hr") +_DROP_TAGS = {"script", "style", "head", "title", "ix:header", "ix:hidden", "xbrli:context"} +_DIG = re.compile(r"\d") +_ALPHA = re.compile(r"[A-Za-z]") +_WORD = re.compile(r"[a-z]{2,}") +_TOC = re.compile(r"^\s*(table\s+of\s+contents|index|page)\s*$", re.I) +_RF_START = re.compile(r"^[ \t]*item[\s.]*1a\b[\s.:\-–—]*risk\s+factors", re.I | re.M) +_RF_END = re.compile(r"^[ \t]*item[\s.]*(?:1b|2)\b", re.I | re.M) +_HTMLISH = re.compile(rb"<(html|body|div|p|font|br|table)[\s>/]", re.I) + + +def _digit_share(s: str) -> float: + d = len(_DIG.findall(s)) + a = len(_ALPHA.findall(s)) + return d / (a + d) if a + d else 0.0 + + +def _html_text(raw: bytes) -> str: + import lxml.html + parser = lxml.html.HTMLParser(huge_tree=True, remove_comments=True, remove_pis=True) + root = lxml.html.document_fromstring(raw, parser=parser) + drop = [] + for el in root.iter(): + if not isinstance(el.tag, str): + continue + if el.tag.lower() in _DROP_TAGS: + drop.append(el) + continue + st = el.get("style") + if st and "display:none" in st.replace(" ", "").lower(): + drop.append(el) + for el in drop: + el.drop_tree() + for tb in reversed(list(root.iter("table"))): # innermost first + if _digit_share(tb.text_content()) > 0.15: + tb.drop_tree() + for el in root.iter(*_BLOCK): + el.tail = "\n" + (el.tail or "") + for el in root.iter("td", "th"): + el.tail = " " + (el.tail or "") + return root.text_content() + + +def _plain_text(raw: bytes) -> str: + s = raw.decode("utf-8", errors="replace") + s = re.sub(r"<[^>\n]{1,40}>", " ", s) # EDGAR SGML tags (, , ) + return html.unescape(s) + + +def _primary_part(raw: bytes) -> bytes: + """The first of a full-submission text file (the form itself).""" + if b"" not in raw[:20000]: + return raw + m = re.search(rb"(.*?)(?:|\Z)", raw, re.S) + body = m.group(1) if m else raw + t = re.search(rb"(.*?)(?:|\Z)", body, re.S) + return t.group(1) if t else body + + +def clean_text(raw: bytes) -> str: + raw = _primary_part(raw) + try: + txt = _html_text(raw) if _HTMLISH.search(raw[:20000]) else _plain_text(raw) + except Exception: # noqa: BLE001 — malformed markup: fall back to tag stripping + txt = _plain_text(re.sub(rb"<[^>]*>", b" ", raw)) + txt = txt.replace("\xa0", " ").replace("\r", "\n") + keep = [] + for line in txt.split("\n"): + line = line.strip() + if not line or _TOC.match(line) or _digit_share(line) > 0.25: + continue + keep.append(line) + return "\n".join(keep) + + +def tokens(text: str) -> Counter: + return Counter(w for w in _WORD.findall(text.lower()) if w not in STOP and len(w) <= 30) + + +def risk_factors(text: str) -> str | None: + """Item 1A of a 10-K: the longest span from a line-leading "Item 1A. Risk + Factors" to the next line-leading Item 1B / Item 2 (the table of contents + gives short spans).""" + ends = [m.start() for m in _RF_END.finditer(text)] + best = None + for m in _RF_START.finditer(text): + e = next((x for x in ends if x > m.end()), None) + if e is None: + continue + if best is None or e - m.end() > best[1] - best[0]: + best = (m.end(), e) + if best is None: + return None + sec = text[best[0]:best[1]] + return sec if len(sec.split()) >= 200 else None + + +def cosine(a: Counter, b: Counter) -> float: + if not a or not b: + return np.nan + if len(a) > len(b): + a, b = b, a + dot = sum(v * b.get(k, 0) for k, v in a.items()) + na = np.sqrt(sum(v * v for v in a.values())) + nb = np.sqrt(sum(v * v for v in b.values())) + return float(dot / (na * nb)) + + +def jaccard(a: Counter, b: Counter) -> float: + if not a or not b: + return np.nan + sa, sb = set(a), set(b) + return len(sa & sb) / len(sa | sb) + + +def doc_vectors(cik: int, accession: str, doc: str, group: str) -> dict | None: + """Term counts of a filing's primary document (+ Item 1A for a 10-K). + None when the download failed (retry later); empty counts when the + document exists but holds no text.""" + acc = accession.replace("-", "") + name = doc if isinstance(doc, str) and doc.strip() else f"{accession}.txt" + r = _get(f"https://www.sec.gov/Archives/edgar/data/{cik}/{acc}/{name}") + if r is None: + return None + if name.lower().endswith(".pdf") or r.headers.get("content-type", "").startswith("application/pdf"): + return {"vec": Counter(), "rf": None} + text = clean_text(r.content) + rf = risk_factors(text) if group == "10-K" else None + return {"vec": tokens(text), "rf": tokens(rf) if rf else None} + + +# --------------------------------------------------------------------------- +# stages 2-3: per-CIK chains of filings +# --------------------------------------------------------------------------- + +def _cache_path(cik: int, group: str) -> str: + return os.path.join(VEC, f"{cik}_{group}.pkl.z") + + +def _load_cache(cik: int, group: str) -> dict: + path = _cache_path(cik, group) + if not os.path.exists(path): + return {} + try: + with open(path, "rb") as fh: + return pickle.loads(zlib.decompress(fh.read())) + except Exception: # noqa: BLE001 — a torn write: rebuild from downloads + return {} + + +def _save_cache(cik: int, group: str, cache: dict) -> None: + path = _cache_path(cik, group) + tmp = path + f".{os.getpid()}.tmp" + with open(tmp, "wb") as fh: + fh.write(zlib.compress(pickle.dumps(cache, protocol=5), 6)) + os.replace(tmp, path) + + +def _append_rows(rows: list[dict]) -> None: + with open(OUT, "a", newline="") as fh: + fcntl.flock(fh, fcntl.LOCK_EX) + w = csv.DictWriter(fh, fieldnames=COLS) + for r in rows: + w.writerow(r) + fh.flush() + fcntl.flock(fh, fcntl.LOCK_UN) + + +def work(job: tuple) -> tuple[int, int, int]: + """Process the pending filings of one CIK in filing order. Returns + (cik, filings written, failures).""" + cik, group, filings, done = job + keep_last = 1 if group == "10-K" else 4 + cache = _load_cache(cik, group) + accs = [f["accession"] for f in filings] + n_ok = n_fail = 0 + for i, f in enumerate(filings): + if f["accession"] in done: + continue + try: + cur = cache.get(f["accession"]) or doc_vectors(cik, f["accession"], f["doc"], group) + if cur is None: + n_fail += 1 + continue + ref = None + if f["ref"]: + ref = cache.get(f["ref"]) + if ref is None: + rf_row = next(g for g in filings if g["accession"] == f["ref"]) + ref = doc_vectors(cik, rf_row["accession"], rf_row["doc"], group) + if ref is None: + n_fail += 1 + continue + cache[f["ref"]] = ref + cache[f["accession"]] = cur + done.add(f["accession"]) + # keep what a pending filing still needs, plus the latest filings + needed = {g["ref"] for g in filings if g["accession"] not in done and g["ref"]} + needed |= set(accs[max(0, i + 1 - keep_last): i + 1]) + cache = {a: v for a, v in cache.items() if a in needed} + _save_cache(cik, group, cache) + v, rv = cur["vec"], (ref or {}).get("vec") + base = {"cik": cik, "form": f["form"], "filed": f["filed"], "period": f["period"], + "accession": f["accession"], "ref_accession": f["ref"] or "", + "n_words": sum(v.values()), + "cos_sim": cosine(v, rv) if rv else np.nan, + "jac_sim": jaccard(v, rv) if rv else np.nan, + "n_words_rf": sum(cur["rf"].values()) if cur.get("rf") else np.nan} + rr = (ref or {}).get("rf") + base["cos_rf"] = cosine(cur["rf"], rr) if cur.get("rf") and rr else np.nan + base["jac_rf"] = jaccard(cur["rf"], rr) if cur.get("rf") and rr else np.nan + for c in ("cos_sim", "jac_sim", "cos_rf", "jac_rf"): + if not np.isnan(base[c]): + base[c] = round(base[c], 6) + _append_rows([{"ticker": t, **base} for t in f["tickers"].split("|")]) + n_ok += 1 + except Exception: # noqa: BLE001 — log and leave it for a rerun + n_fail += 1 + print(f" ! {cik} {f['accession']}:\n{traceback.format_exc()}", flush=True) + return cik, n_ok, n_fail + + +def run_group(ix: pd.DataFrame, group: str, workers: int, lock, nxt) -> None: + f = with_refs(ix, group) + f["ref"] = f["ref"].fillna("") + f["doc"] = f["doc"].fillna("") + f["period"] = f["period"].fillna("") + done = set() + if os.path.exists(OUT): + done = set(pd.read_csv(OUT, usecols=["accession"])["accession"]) + jobs = [] + for cik, g in f.groupby("cik"): + recs = g[["form", "filed", "period", "accession", "ref", "doc", "tickers"]].to_dict("records") + d = {r["accession"] for r in recs if r["accession"] in done} + if len(d) < len(recs): + jobs.append((int(cik), group, recs, d)) + total = sum(len(j[2]) - len(j[3]) for j in jobs) + print(f"=== {group}: {len(f)} filings, {len(f) - total} done, {total} to fetch " + f"({len(jobs)} CIKs); {(f['ref'] != '').mean():.1%} have a comparable prior filing", + flush=True) + if not jobs: + return + jobs.sort(key=lambda j: -(len(j[2]) - len(j[3]))) # long chains first + t0, n, fails = time.time(), 0, 0 + with mp.get_context("fork").Pool(workers, initializer=_init_worker, initargs=(lock, nxt)) as pool: + for k, (cik, ok, bad) in enumerate(pool.imap_unordered(work, jobs), 1): + n += ok + fails += bad + el = time.time() - t0 + if k % 10 == 0 or k == len(jobs): + eta = (total - n) / max(n / el, 1e-9) / 3600 + print(f" {group} {k}/{len(jobs)} CIKs, {n}/{total} filings, {fails} failed, " + f"{n / el:.2f}/s, ETA {eta:.1f}h [{el / 60:.0f} min]", flush=True) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--forms", nargs="*", default=["10-K", "10-Q"], choices=list(GROUPS)) + ap.add_argument("--workers", type=int, default=8) + ap.add_argument("--refresh-index", action="store_true") + args = ap.parse_args() + os.makedirs(VEC, exist_ok=True) + lock, nxt = mp.Lock(), mp.Value("d", 0.0) + _init_worker(lock, nxt) + if args.refresh_index or not os.path.exists(INDEX): + with open(os.path.join(DIR, "us_ohlcv.pkl"), "rb") as fh: + tickers = [t for t in pickle.load(fh)["close"].columns if t not in US_ETFS] + build_index(tickers) + ix = pd.read_csv(INDEX, dtype={"period": str, "doc": str}) + if not os.path.exists(OUT): + with open(OUT, "w", newline="") as fh: + csv.DictWriter(fh, fieldnames=COLS).writeheader() + for group in args.forms: + run_group(ix, group, args.workers, lock, nxt) + print("done", flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_long_history.py b/research/lab/fetch_long_history.py new file mode 100644 index 0000000..f3d5773 --- /dev/null +++ b/research/lab/fetch_long_history.py @@ -0,0 +1,107 @@ +""" +Long daily history for the leveraged-ETF and cross-asset research +(research/lab/letf_history.py): stock indices from 1927/1985, the ETFs the +trend strategies trade, gold futures, bitcoin, and Treasury yields. + + PYTHONPATH=. uv run python -m research.lab.fetch_long_history [--no-proxy] + +Yahoo is fetched ticker by ticker (batched requests fail intermittently) +and each ticker is cached in data/research/ext/yahoo/ as it arrives, so a +re-run after rate limiting fetches only what is missing. The fetch stops +after three tickers fail in a row: Yahoo is rate limiting, and the live +monitor needs it too. FRED yields come from the public fredgraph CSV +endpoint. Writes data/research/long_history.pkl: {"close", "open"} +DataFrames (adjusted for splits and dividends) and "fred" (yields in +percent, and the Nasdaq-100 level). +""" + +from __future__ import annotations + +import argparse +import io +import os +import pickle +import time +import urllib.request +import warnings + +import pandas as pd +import yfinance as yf + +from research.lab.fetch_us import _proxy_env + +warnings.filterwarnings("ignore") +OUT = os.path.join("data", "research", "long_history.pkl") +YAHOO = [ + "^GSPC", "^NDX", "^VIX", # indices (price only) + "SPY", "QQQ", "IWM", "EFA", "EEM", "VNQ", "TLT", "IEF", "SHY", "GLD", "DBC", "GSG", "UUP", + "TQQQ", "UPRO", "QLD", "SSO", "SOXL", "TMF", "UGL", # leveraged ETFs (for validation) + "VFINX", "VUSTX", # mutual funds with dividends, 1980s- + "GC=F", "BTC-USD", "ETH-USD", "IBIT", +] +FRED = ["DTB3", "DGS2", "DGS10", "DGS30", "NASDAQ100"] + + +CACHE = os.path.join("data", "research", "ext", "yahoo") + + +def _cached(t: str) -> str: + return os.path.join(CACHE, t.replace("^", "_").replace("=", "_") + ".csv") + + +def fetch_yahoo(tickers: list[str], start: str = "1927-01-01") -> tuple[pd.DataFrame, pd.DataFrame]: + os.makedirs(CACHE, exist_ok=True) + fails = 0 + for t in tickers: + if os.path.exists(_cached(t)) or fails >= 3: + continue + for attempt in range(2): + try: + h = yf.Ticker(t).history(start=start, auto_adjust=True) + if len(h): + h.index = h.index.tz_localize(None).normalize() + h[["Open", "Close"]].to_csv(_cached(t)) + break + except Exception as exc: # noqa: BLE001 — retry transient Yahoo errors + print(f" {t}: {type(exc).__name__}", flush=True) + time.sleep(20) + fails = 0 if os.path.exists(_cached(t)) else fails + 1 + print(f"{t:8s} {'ok' if os.path.exists(_cached(t)) else 'MISSING'}", flush=True) + time.sleep(2) + close, opens = {}, {} + for t in tickers: + if os.path.exists(_cached(t)): + h = pd.read_csv(_cached(t), index_col=0, parse_dates=True) + close[t], opens[t] = h["Close"], h["Open"] + missing = [t for t in tickers if t not in close] + if missing: + print("missing (re-run later):", missing) + return pd.DataFrame(close).sort_index(), pd.DataFrame(opens).sort_index() + + +def fetch_fred(series: list[str]) -> pd.DataFrame: + out = {} + for s in series: + url = f"https://fred.stlouisfed.org/graph/fredgraph.csv?id={s}" + with urllib.request.urlopen(url, timeout=60) as fh: + df = pd.read_csv(io.BytesIO(fh.read()), index_col=0, parse_dates=True) + out[s] = pd.to_numeric(df.iloc[:, 0], errors="coerce") + return pd.DataFrame(out) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--no-proxy", action="store_true", + help="connect directly (when the proxy's IP is rate limited by Yahoo)") + if not ap.parse_args().no_proxy: + _proxy_env() + close, opens = fetch_yahoo(YAHOO) + fred = fetch_fred(FRED) + os.makedirs(os.path.dirname(OUT), exist_ok=True) + with open(OUT, "wb") as fh: + pickle.dump({"close": close, "open": opens, "fred": fred}, fh) + print(f"Saved {OUT}: {close.shape[1]} series, FRED {list(fred.columns)}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_us_mid.py b/research/lab/fetch_us_mid.py new file mode 100644 index 0000000..ac46a2d --- /dev/null +++ b/research/lab/fetch_us_mid.py @@ -0,0 +1,462 @@ +""" +Point-in-time S&P MidCap 400 / SmallCap 600 research universe: month-end +holdings of the iShares Core ETFs that track them (IJH, IJR), and Yahoo daily +bars for every ticker that was ever held. research/lab/panel_mid.py builds +the Panels. + + PYTHONPATH=. uv run python -m research.lab.fetch_us_mid holdings # ~45 min; then parses + PYTHONPATH=. uv run python -m research.lab.fetch_us_mid prices # ~3 h; resumable + PYTHONPATH=. uv run python -m research.lab.fetch_us_mid parse # re-parse cached CSVs + PYTHONPATH=. uv run python -m research.lab.fetch_us_mid status # counts, weight sums + +Holdings. BlackRock's product-page API serves past month-end snapshots +(``get-fund-document?...&asOfDate=YYYYMMDD&component=holdings``; only the +last trading day of a month works) from 2006-09-29 for both funds, except +2017-01..2017-06, which return nothing. Raw CSVs are cached in +data/research/ext/ishares//.csv; months that return nothing are +remembered in index.json. Requests go out directly (no proxy), one per ~1.2 s. + +Old snapshots list each security under its LAST ticker in BlackRock's +security master (Westwood One is "DIAL" in 2006), with internal codes +(``WM1`` is Williams-Sonoma, ``VRX.1``, ``PSYS.``), no weights (quantity and +price only), and a few rows without a ticker. ``_resolve`` maps codes and +missing tickers to the clean symbol the company name carries in other +snapshots, else strips the suffix; panel_mid's price check catches what is +still wrong. Rights, warrants and stub positions (under 2% of an equal +weight: shares left over after a deletion) are dropped and weights are +recomputed from quantity x price. Outputs: holdings.csv (one row per fund, +date, ticker) and snapshots.csv (count, equity value, shares outstanding). + +Prices. Yahoo daily bars from 2005 (auto_adjust=False so the as-traded close +can be checked against the holdings file's price), one ticker per request, +cached in data/research/ext/yahoo_mid/.csv; tickers Yahoo does not +serve (mostly delisted names) are listed in missing.json. The connection is +the shell's default one (Yahoo blocks this host's own IP); the proxy in +QUANT_MARKET_PROXY is removed from the environment, because the live monitor +fetches through it. Pauses 1.5 s between tickers, 60 s after an error, and +10 minutes after five errors in a row; a run of 25 misses is re-checked +against SPY in case Yahoo is throttling. +""" + +from __future__ import annotations + +import argparse +import csv +import io +import json +import os +import re +import time +import urllib.request +from datetime import date, datetime, timedelta + +import numpy as np +import pandas as pd + +ROOT = os.path.join("data", "research", "ext") +HOLD_DIR = os.path.join(ROOT, "ishares") +PX_DIR = os.path.join(ROOT, "yahoo_mid") +HOLDINGS = os.path.join(HOLD_DIR, "holdings.csv") +FUNDS = {"IJH": 239763, "IJR": 239774} # S&P MidCap 400, S&P SmallCap 600 +ETFS = ["IJH", "IJR", "MDY", "SLY", "IWM", "SPY", "QQQ", "BIL", "SHY"] +API = ("https://www.ishares.com/varnish-api/blk-one01-product-data/product-data/api/" + "v1/get-fund-document?appType=PRODUCT_PAGE&appSubType=ISHARES&targetSite=us-ishares" + "&locale=en_US&portfolioId={pid}&userType=individual&asOfDate={d}&component=holdings") +DATES_API = ("https://www.ishares.com/varnish-api/blk-one01-product-data/product-data/api/" + "v2/get-product-data?appSubType=ISHARES&appType=PRODUCT_PAGE&component=holdings" + "&locale=en_US&portfolioId={pid}&targetSite=us-ishares&userType=individual" + "&excludeContent=true&includeConfig=true") +UA = "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/124 Safari/537.36" +FIRST_MONTH = "2006-01" +NOT_STOCK = r"\bRIGHTS?\b|\bRTS\b|WARRANT|\bWTS?\b|\bCVR\b|ESCROW|CONTRA\b|CASH ACCRUAL|VESTING" +STUB_REL = 0.02 # positions below 2% of an equal weight are leftovers +CASH_SECTORS = {"S-T Securities", "Cash and/or Derivatives", "Cash", "Money Market"} + + +# --------------------------------------------------------------------------- +# holdings +# --------------------------------------------------------------------------- + +_direct = urllib.request.build_opener(urllib.request.ProxyHandler({})) + + +def _get(url: str, headers: dict | None = None) -> bytes: + req = urllib.request.Request(url, headers={"User-Agent": UA, **(headers or {})}) + for attempt in range(4): + try: + with _direct.open(req, timeout=60) as fh: + return fh.read() + except Exception as exc: # noqa: BLE001 — retry transient network errors + print(f" {type(exc).__name__}: {exc}; retry {attempt + 1}", flush=True) + time.sleep(10 * (attempt + 1)) + raise RuntimeError(f"giving up on {url}") + + +def _month_end_candidates(month: pd.Period) -> list[date]: + """The last four weekdays of a month, latest first (the last trading day + is one of them).""" + d = month.end_time.date() + out = [] + while len(out) < 4: + if d.weekday() < 5: + out.append(d) + d -= timedelta(days=1) + return out + + +def _as_of(text: str) -> str | None: + m = re.search(r'Fund Holdings as of,"([^"]+)"', text) + if not m or m.group(1).strip() == "-": + return None + return datetime.strptime(m.group(1), "%b %d, %Y").strftime("%Y%m%d") + + +def fetch_holdings() -> None: + os.makedirs(HOLD_DIR, exist_ok=True) + idx_path = os.path.join(HOLD_DIR, "index.json") + index = json.load(open(idx_path)) if os.path.exists(idx_path) else {} + today = date.today() + for fund, pid in FUNDS.items(): + fdir = os.path.join(HOLD_DIR, fund) + os.makedirs(fdir, exist_ok=True) + seen = index.setdefault(fund, {}) + months = pd.period_range(FIRST_MONTH, pd.Period(today, "M") - 1, freq="M") + # recent dates the product page lists (the latest may be mid-month) + meta = json.loads(_get(DATES_API.format(pid=pid), {"Accept": "application/json", + "x-application-id": "pp-ui-csr"})) + listed = meta["componentsByNameMap"]["holdings"]["containersByNameMap"]["all"][ + "dataPointsByNameMap"]["dateList"]["value"] + time.sleep(1.2) + todo = [(str(m), _month_end_candidates(m)) for m in months] + todo += [(f"listed-{d}", [datetime.strptime(str(d), "%Y%m%d").date()]) for d in listed] + for key, cands in todo: + if key in seen: + continue + got = None + for d in cands: + if d >= today: + continue + ds = d.strftime("%Y%m%d") + path = os.path.join(fdir, f"{ds}.csv") + if os.path.exists(path): + got = ds + break + text = _get(API.format(pid=pid, d=ds)).decode("utf-8-sig", "replace") + time.sleep(1.2) + asof = _as_of(text) + if asof is None: + continue + if asof != ds: + print(f" {fund} {ds}: file says {asof}", flush=True) + with open(os.path.join(fdir, f"{asof}.csv"), "w") as fh: + fh.write(text) + got = asof + break + seen[key] = got + print(f"{fund} {key}: {got or 'none'}", flush=True) + json.dump(index, open(idx_path, "w"), indent=1, sort_keys=True) + json.dump(index, open(idx_path, "w"), indent=1, sort_keys=True) + parse_all() + + +def _num(x: str) -> float: + x = (x or "").replace(",", "").strip() + try: + return float(x) + except ValueError: + return np.nan + + +def parse_snapshot(path: str) -> pd.DataFrame: + """All rows of one holdings CSV, numbers parsed.""" + rows = list(csv.reader(io.open(path, encoding="utf-8-sig"))) + h = next(i for i, r in enumerate(rows) if r and r[0] == "Ticker") + shares_out = next((_num(r[1]) for r in rows[:h] if len(r) > 1 and r[0] == "Shares Outstanding"), + np.nan) + head = rows[h] + body = [r for r in rows[h + 1:] if len(r) == len(head)] + df = pd.DataFrame(body, columns=head) + wcol = next((c for c in ("Weight (%)", "Market Weight") if c in df.columns), None) + return pd.DataFrame({ + "shares_out": shares_out, + "ticker_raw": df["Ticker"].str.strip(), + "name": df["Name"].str.strip().str.upper(), + "sector": df["Sector"].str.strip(), + "asset_class": df["Asset Class"].str.strip(), + "type": df["Type"].str.strip() if "Type" in df else "-", + "quantity": df["Quantity"].map(_num) if "Quantity" in df else np.nan, + "price": df["Price"].map(_num) if "Price" in df else np.nan, + "market_value": df["Market Value"].map(_num) if "Market Value" in df else np.nan, + "weight_reported": df[wcol].map(_num) if wcol else np.nan, + }) + + +def _is_equity(df: pd.DataFrame) -> pd.Series: + """Stocks only: drops cash, money-market funds, futures, swaps and + collateral. A few old rows have no asset class (and no ticker) but a GICS + sector and a share price; they are stocks.""" + typ_ok = df["type"].isin(["EQUITY", "-", ""]) + listed = df["asset_class"] == "Equity" + unclassed = (df["asset_class"].isin(["-", ""]) & ~df["sector"].isin(CASH_SECTORS | {"-", ""}) + & (df["price"] > 1.0) & (df["price"] != 100.0)) + not_stock = df["name"].str.contains(NOT_STOCK, regex=True, na=False) + return typ_ok & (listed | unclassed) & ~not_stock + + +_CLEAN = re.compile(r"[A-Z]{1,5}([ .][A-Z])?") + + +def _yahoo_symbol(raw: str) -> str | None: + """BlackRock ticker → Yahoo symbol. Class shares use a dash (BRK.B → + BRK-B, "MOG A" → MOG-A); BlackRock's disambiguation suffixes are dropped + (VRX.1 → VRX, PSYS. → PSYS, LKQ1 → LKQ, "CC WI" → CC). For suffixed + codes this is the last resort: ``_resolve`` first looks the company name + up (WM1 is Williams-Sonoma, WSM).""" + t = re.sub(r"\s+", " ", raw.strip().upper()) + if t in ("", "-") or re.search(r"[ ./](WS|RT|RTS)$", t): # rights and warrants + return None + t = re.sub(r"[ .](WI|US)$", "", t) # when-issued lines, "EHTH US" + m = re.fullmatch(r"([A-Z]+)[ .]([A-Z])", t) + if m: + return f"{m.group(1)}-{m.group(2)}" + t = re.sub(r"\.\d*$", "", t) # VRX.1, PSYS. + m = re.fullmatch(r"([A-Z]+)\d+([A-Z]?)", t) # LKQ1, RS6, FR1A + if m: + t = m.group(1) + (f"-{m.group(2)}" if m.group(2) else "") + t = t.replace(".", "-").replace(" ", "-").replace("/", "-") + return t or None + + +def _norm_name(n: str) -> str: + n = str(n).upper().replace("&", " AND ") + n = re.sub(r"\bCL(ASS)?\s+([A-Z])\b", r"CL\2", n) + n = re.sub(r"[^A-Z0-9 ]", " ", n) + n = re.sub(r"\b(INC|INCORPORATED|CORP|CORPORATION|CO|COMPANY|LTD|LIMITED|PLC|THE|NEW|" + r"HOLDINGS?|HLDGS?|GROUP|GRP|INTL|INTERNATIONAL)\b", " ", n) + return " ".join(n.split()) + + +def _resolve(rows: pd.DataFrame) -> pd.Series: + """Yahoo symbol per row. Rows whose ticker is a clean symbol keep it; + suffixed codes and missing tickers take the clean symbol the same + company name carries in any snapshot (exact normalized name, else a + unique prefix match of at least 10 characters), else the stripped code.""" + raw = rows["ticker_raw"].fillna("").str.strip().str.upper() + clean = raw.str.fullmatch(_CLEAN) + nm = rows["name"].map(_norm_name) + sym = raw.map(_yahoo_symbol) + by_name = (pd.DataFrame({"nm": nm[clean], "sym": sym[clean]}) + .groupby("nm")["sym"].agg(lambda s: s.mode().iloc[0])) + names = sorted(by_name.index) + out = sym.where(clean) + for n in nm[~clean].unique(): + hit = by_name.get(n) + if hit is None and len(n) >= 10: + cand = {by_name[k] for k in names if k.startswith(n) or (len(k) >= 10 and n.startswith(k))} + hit = cand.pop() if len(cand) == 1 else None + if hit is not None: + out[(~clean) & (nm == n)] = hit + left = out.isna() + out[left] = sym[left] + return out + + +def parse_all() -> pd.DataFrame: + parts = [] + for fund in FUNDS: + fdir = os.path.join(HOLD_DIR, fund) + for f in sorted(os.listdir(fdir)) if os.path.isdir(fdir) else []: + if f.endswith(".csv"): + df = parse_snapshot(os.path.join(fdir, f)) + df.insert(0, "date", pd.Timestamp(f[:-4])) + df.insert(0, "fund", fund) + parts.append(df) + allrows = pd.concat(parts, ignore_index=True) + eq = allrows[_is_equity(allrows)].copy() + eq["mv"] = eq["quantity"] * eq["price"] + eq["mv"] = eq["market_value"].where(eq["market_value"] > 0, eq["mv"]) + eq["weight"] = eq["mv"] / eq.groupby(["fund", "date"])["mv"].transform("sum") + eq["ticker"] = _resolve(eq) + unresolved = eq[eq["ticker"].isna()] + if len(unresolved): + print(f"{len(unresolved)} equity rows without a ticker (dropped):", + sorted(unresolved["name"].unique())[:20]) + eq = eq.dropna(subset=["ticker"]).sort_values("mv", ascending=False) + # one row per ticker and snapshot (a security held in two lots is summed; + # the price is the largest lot's) + eq = (eq.groupby(["fund", "date", "ticker"], as_index=False) + .agg(name=("name", "first"), ticker_raw=("ticker_raw", "first"), + sector=("sector", "first"), quantity=("quantity", "sum"), price=("price", "first"), + mv=("mv", "sum"), weight=("weight", "sum"), + weight_reported=("weight_reported", "sum"), shares_out=("shares_out", "first"))) + # stub positions: a few hundred shares left after a deletion or a + # corporate action, not index members (crashed members stay above 2%) + rel = eq["weight"] * eq.groupby(["fund", "date"])["ticker"].transform("size") + stub = rel < STUB_REL + print(f"dropped {stub.sum()} stub positions (< {STUB_REL:.0%} of an equal weight)") + eq = eq[~stub].copy() + eq["weight"] = eq["weight"] / eq.groupby(["fund", "date"])["weight"].transform("sum") + eq.drop(columns="shares_out").to_csv(HOLDINGS, index=False) + summ = eq.groupby(["fund", "date"]).agg(n=("ticker", "size"), equity_mv=("mv", "sum"), + weight_reported=("weight_reported", "sum"), + shares_out=("shares_out", "first")) + summ.to_csv(os.path.join(HOLD_DIR, "snapshots.csv")) + print(f"Saved {HOLDINGS}: {len(eq)} rows, {eq['ticker'].nunique()} tickers") + return eq + + +def load_holdings() -> pd.DataFrame: + return pd.read_csv(HOLDINGS, parse_dates=["date"]) + + +# --------------------------------------------------------------------------- +# prices +# --------------------------------------------------------------------------- + +def _px_path(t: str) -> str: + return os.path.join(PX_DIR, t.replace("^", "_").replace("=", "_").replace("/", "_") + ".csv") + + +def _no_market_proxy() -> None: + """Never fetch Yahoo through QUANT_MARKET_PROXY (the live monitor's IP).""" + mp = os.environ.get("QUANT_MARKET_PROXY", "") + if not mp and os.path.exists(".env"): + for line in open(".env"): + if line.startswith("QUANT_MARKET_PROXY="): + mp = line.split("=", 1)[1].strip().strip("'\"") + for k in ("http_proxy", "https_proxy", "HTTP_PROXY", "HTTPS_PROXY", "all_proxy", "ALL_PROXY"): + if mp and os.environ.get(k, "").rstrip("/") == mp.rstrip("/"): + del os.environ[k] + + +def fetch_prices(start: str = "2005-01-01", retry_missing: bool = False) -> None: + import yfinance as yf + from yfinance.exceptions import YFPricesMissingError, YFTzMissingError + yf.config.debug.hide_exceptions = False + _no_market_proxy() + os.makedirs(PX_DIR, exist_ok=True) + miss_path = os.path.join(PX_DIR, "missing.json") + missing = json.load(open(miss_path)) if os.path.exists(miss_path) else {} + tickers = ETFS + sorted(set(load_holdings()["ticker"]) - set(ETFS)) + todo = [t for t in tickers if not os.path.exists(_px_path(t)) + and (retry_missing or t not in missing)] + print(f"{len(tickers)} tickers, {len(todo)} to fetch", flush=True) + errors = 0 + miss_streak: list[str] = [] + i = 0 + while i < len(todo): + t = todo[i] + try: + h = yf.Ticker(t).history(start=start, auto_adjust=False, actions=True) + if len(h) == 0: + raise YFPricesMissingError(t, "empty") + h.index = pd.DatetimeIndex(h.index.tz_localize(None).normalize(), name="Date") + h.to_csv(_px_path(t)) + missing.pop(t, None) + print(f"[{i + 1}/{len(todo)}] {t:8s} ok {len(h)} bars {h.index[0].date()}", flush=True) + errors, miss_streak = 0, [] + except Exception as exc: # noqa: BLE001 — sorted below + msg = f"{type(exc).__name__}: {exc}" + gone = (isinstance(exc, (YFTzMissingError, YFPricesMissingError)) + or "404" in msg or "Not Found" in msg or "delisted" in msg) + if not gone: # rate limit or network error: wait, retry the same ticker + errors += 1 + print(f"[{i + 1}/{len(todo)}] {t:8s} error {msg[:140]}", flush=True) + json.dump(missing, open(miss_path, "w"), indent=0, sort_keys=True) + time.sleep(600 if errors >= 5 else 60) + if errors >= 5: + errors = 0 + continue + missing[t] = datetime.now().isoformat(timespec="seconds") + miss_streak.append(t) + print(f"[{i + 1}/{len(todo)}] {t:8s} missing", flush=True) + if len(miss_streak) >= 25: # a long run of misses: is Yahoo throttling? + try: + ok = len(yf.Ticker("SPY").history(period="5d")) > 0 + except Exception: # noqa: BLE001 + ok = False + if not ok: + print(" canary SPY failed: throttled; redoing the streak after 10 min", flush=True) + for m in miss_streak: + missing.pop(m, None) + i -= len(miss_streak) - 1 + miss_streak = [] + time.sleep(600) + continue + miss_streak = [] + errors = 0 + if (i + 1) % 25 == 0: + json.dump(missing, open(miss_path, "w"), indent=0, sort_keys=True) + i += 1 + time.sleep(1.5) + json.dump(missing, open(miss_path, "w"), indent=0, sort_keys=True) + print("done", flush=True) + + +def load_bars(t: str) -> pd.DataFrame | None: + """Adjusted OHLCV (as yfinance's auto_adjust=True) plus the unadjusted close.""" + path = _px_path(t) + if not os.path.exists(path): + return None + h = pd.read_csv(path, index_col=0, parse_dates=True) + h = h[~h.index.duplicated(keep="last")].sort_index() + ratio = (h["Adj Close"] / h["Close"]).where(h["Close"] > 0) + out = pd.DataFrame({k.lower(): h[k] * ratio for k in ("Open", "High", "Low", "Close")}) + out["volume"] = h["Volume"] + # Yahoo's Close is split-adjusted: undo later splits for the as-traded price + splits = h["Stock Splits"].replace(0, 1.0).fillna(1.0) + later = splits[::-1].cumprod()[::-1].shift(-1).fillna(1.0) + out["raw_close"] = h["Close"] * later + return out + + +def status() -> None: + eq = load_holdings() + tick = set(eq["ticker"]) + have = {t for t in tick if os.path.exists(_px_path(t))} + miss_path = os.path.join(PX_DIR, "missing.json") + missing = set(json.load(open(miss_path))) if os.path.exists(miss_path) else set() + print(f"tickers {len(tick)}: prices {len(have)}, missing {len(missing & tick)}, " + f"todo {len(tick - have - missing)}") + print(snapshot_summary().groupby("fund") + .agg(first=("date", "min"), last=("date", "max"), snaps=("date", "size"), + n_min=("n", "min"), n_max=("n", "max"), + wsum_min=("weight_sum", "min"), wsum_med=("weight_sum", "median"), + wsum_max=("weight_sum", "max"))) + + +def snapshot_summary() -> pd.DataFrame: + """Per snapshot: equity holdings count and their weight sum, the equity + market value over the fund's assets (the files' own cash rows are + unusable: cash is sometimes priced at 100). Assets = shares outstanding + x the ETF's close, both on today's split basis (the header's share count + is split-adjusted, like Yahoo's unadjusted-for-dividends Close).""" + s = pd.read_csv(os.path.join(HOLD_DIR, "snapshots.csv"), parse_dates=["date"]) + s["weight_sum"] = np.nan + for fund in FUNDS: + if not os.path.exists(_px_path(fund)): + continue + close = pd.read_csv(_px_path(fund), index_col=0, parse_dates=True)["Close"] + m = s["fund"] == fund + px = close.reindex(close.index.union(s.loc[m, "date"])).ffill().reindex(s.loc[m, "date"]) + s.loc[m, "weight_sum"] = s.loc[m, "equity_mv"].values / (s.loc[m, "shares_out"].values * px.values) + return s + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("step", choices=["holdings", "parse", "prices", "status"]) + ap.add_argument("--retry-missing", action="store_true") + args = ap.parse_args() + if args.step == "holdings": + fetch_holdings() + elif args.step == "parse": + parse_all() + elif args.step == "prices": + fetch_prices(retry_missing=args.retry_missing) + else: + status() + + +if __name__ == "__main__": + main() diff --git a/research/lab/industry.py b/research/lab/industry.py new file mode 100644 index 0000000..6ea277b --- /dev/null +++ b/research/lab/industry.py @@ -0,0 +1,83 @@ +""" +Industry momentum (Moskowitz & Grinblatt 1999) on S&P 500 members, with +industries from SEC SIC codes (data/research/sec_sic.csv, written by +research/lab/fetch_lazy_prices.py) grouped to 2-digit SIC. + + PYTHONPATH=. uv run --with numba python -m research.lab.industry + +Industry return = equal-weight mean of its PIT members' daily returns (groups +with fewer than 3 members that day are left out). Scores, known at each close: + ind_mom_6_1 industry return over months 6..1 (skip the last month) + ind_mom_12_1 same over 12..1 +Uses: (a) a sleeve holding the top-20 stocks by their industry's momentum, +(b) the momentum sleeve restricted to stocks in the top half of industries. +Alpha against QQQ, 2006-2016 vs 2017-2026. +""" + +from __future__ import annotations + +import dataclasses +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factors import _ret + +warnings.filterwarnings("ignore") +SIC = os.path.join("data", "research", "sec_sic.csv") + + +def industries(p) -> pd.Series: + s = pd.read_csv(SIC).drop_duplicates("ticker").set_index("ticker")["sic"] + s = pd.to_numeric(s, errors="coerce").dropna().astype(int) // 100 + return s.reindex(p.close.columns) + + +def industry_momentum(p, lookback: int, skip: int = 21) -> pd.DataFrame: + """Each stock gets its industry's cumulative return over [t-lookback, t-skip].""" + ind = industries(p) + r = _ret(p).where(p.member) + out = pd.DataFrame(np.nan, index=p.dates, columns=p.close.columns) + for g in ind.dropna().unique(): + cols = list(ind.index[ind == g]) + n = r[cols].notna().sum(axis=1) + gr = r[cols].mean(axis=1).where(n >= 3) + cum = np.log1p(gr.fillna(0.0)).rolling(lookback - skip, min_periods=(lookback - skip) // 2).sum().shift(skip) + out[cols] = np.repeat(cum.where(n >= 3).to_numpy()[:, None], len(cols), axis=1) + return out.where(p.member) + + +def main() -> None: + from research.lab.alpha import alpha_stats, bench_returns, fmt_alpha + from research.lab.factor_scan import load + from research.lab.rotation import factor_ranks + from research.lab.us_sweep import _blend, run + + p, _ = load("us") + ind = industries(p) + print(f"{ind.notna().sum()}/{len(ind)} tickers with SIC; {ind.nunique()} 2-digit groups") + b = bench_returns(p.dates) + pr = factor_ranks(p, "us") + mom = _blend(pr["overnight_mom"], pr["rec_126"]) + span = [("2006-09-27", "2026-09-24")] + + def show(name, r): + a, d, h = alpha_stats(r, b), alpha_stats(r.loc[:"2016"], b), alpha_stats(r.loc["2017":], b) + print(f"{fmt_alpha(name, a)} | 06-16 a {d['alpha']*100:5.1f}% t {d['t']:5.2f} | 17-26 a {h['alpha']*100:5.1f}% t {h['t']:5.2f}", flush=True) + + for lb in (126, 252): + im = industry_momentum(p, lb) + rk = im.rank(axis=1, pct=True) + show(f"industry mom {lb}d top20 e21", run(p, rk, span, top_n=20, buffer=2.0, every=21)["_rets"]) + tie = rk + 0.001 * pr["mom_12_1"].fillna(0) # stocks of the best industries, by own momentum + show(f"industry mom {lb}d (+stock tiebreak) top20", run(p, tie.where(p.member), span, top_n=20, buffer=2.0, every=21)["_rets"]) + keep = rk >= 0.5 + pm = dataclasses.replace(p, member=p.member & keep.fillna(False)) + show(f"momentum top10, top-half industries ({lb}d)", run(pm, mom.where(pm.member), span, top_n=10, buffer=2.0, every=10)["_rets"]) + show("momentum top10 (reference)", run(p, mom.where(p.member), span, top_n=10, buffer=2.0, every=10)["_rets"]) + + +if __name__ == "__main__": + main() diff --git a/research/lab/insider.py b/research/lab/insider.py new file mode 100644 index 0000000..6e8a424 --- /dev/null +++ b/research/lab/insider.py @@ -0,0 +1,432 @@ +""" +Point-in-time insider-trading factors for the US panel, from officers' and +directors' open-market Form 4 purchases (P) and sales (S) +(research/lab/fetch_insider.py, SEC Insider Transactions Data Sets, 2006 →). + +A trade is usable from the first trading day strictly after its filing date, +so a factor on date t only counts filings dated before t (the engine then +trades at the open of t+1). A window of W calendar days counts a filing on t +when t - W <= filed < t. + +Factors (dates x panel tickers). NaN where no issuer CIK covers the ticker on +that date (``covered``), 0 where it is covered but nothing was filed: + insider_buyers_90d distinct insiders with a purchase filed in the last + 90 days minus distinct insiders with a sale + insider_net_value_180d (purchase $ - sale $) filed in the last 180 days / + market cap; NaN where fundamentals' mcap is missing + (all names before 2009, ~10% after) + insider_opportunistic opportunistic buyers minus opportunistic sellers, + last 90 days (Cohen, Malloy & Pomorski 2012); NaN + before 2009-04-01 (needs three prior years of trades) + insider_cluster_buy 1 if >= 2 distinct officers/directors bought in the + last 30 days, else 0 + +Routine vs opportunistic (CMP 2012, "Decoding Inside Information", p. 11-12): +an insider who traded in each of the three preceding calendar years is +"routine" for year Y if some calendar month m had a trade in each of Y-1, Y-2 +and Y-3, else "opportunistic"; all her trades in Y get that label. Insiders +without a trade in each of the three years are unclassified. The paper does +not condition on direction; ``same_direction=True`` is the variant that +counts only same-direction trades for the routine test. Classification is per +(insider, issuer) and, stricter than the paper, only uses history filed +before the classified trade's filing date. + + from research.lab import insider + I = insider.load(p) # cached + PYTHONPATH=. uv run python -m research.lab.insider # evaluation +""" + +from __future__ import annotations + +import json +import os +import pickle +import time +import warnings + +import numpy as np +import pandas as pd + +from research.lab import fundamentals as fu +from research.lab.factor_scan import RESULTS +from research.lab.factor_scan import load as load_panel +from research.lab.report import stats + +warnings.filterwarnings("ignore") + +TRADES = os.path.join("data", "research", "insider_trades.pkl") +CACHE = os.path.join("data", "research", "insider_factors.pkl") +OPP_START = "2009-04-01" # first 90-day window made only of classifiable (2009+) trades +COVER_SLACK = 120 # days after an issuer's last Form 4 that still count as covered +NAMES = ["insider_buyers_90d", "insider_net_value_180d", "insider_opportunistic", "insider_cluster_buy"] + + +def _raw() -> tuple[pd.DataFrame, pd.DataFrame]: + with open(TRADES, "rb") as fh: + raw = pickle.load(fh) + return raw["trades"], raw["tickers"] + + +# --------------------------------------------------------------------------- +# trades → events +# --------------------------------------------------------------------------- + +def classify_cmp(tr: pd.DataFrame, same_direction: bool = False) -> pd.Series: + """'routine' / 'opportunistic' / None for every trade (index of ``tr``).""" + d = tr[["insider_cik", "cik", "code", "trans_date", "filed"]].copy() + d["y"] = d["trans_date"].dt.year + d["m"] = d["trans_date"].dt.month + ok = d["y"].notna() + d = d[ok].astype({"y": int, "m": int}) + d["fd"] = (d["filed"] - pd.Timestamp("1970-01-01")).dt.days.astype(float) # NaN compares False + who = ["insider_cik", "cik"] + mkey = who + (["code"] if same_direction else []) + # first filing day of any trade per year, and of a trade per (year, month) + any_y = d.groupby(who + ["y"])["fd"].min() + mon = d.groupby(mkey + ["y", "m"])["fd"].min().unstack("m").reindex(columns=range(1, 13)) + f = d["fd"].values + hist = np.ones(len(d), bool) + routine = np.ones((len(d), 12), bool) + for k in (1, 2, 3): + a = any_y.reindex(pd.MultiIndex.from_arrays([d[c] for c in who] + [d["y"] - k])).values + hist &= a < f + M = mon.reindex(pd.MultiIndex.from_arrays([d[c] for c in mkey] + [d["y"] - k])).values + routine &= M < f[:, None] + lab = np.where(hist & routine.any(axis=1), "routine", np.where(hist, "opportunistic", None)) + return pd.Series(lab, index=d.index).reindex(tr.index) + + +def _explode(tr: pd.DataFrame, tm: pd.DataFrame) -> pd.DataFrame: + """One row per (trade, panel ticker): every share class of an issuer CIK + (GOOG/GOOGL, FOX/FOXA, ...) gets that issuer's trades.""" + alias = tm[tm["source"] == "sec_map"].groupby("cik")["ticker"].apply(sorted) + multi = alias[alias.str.len() > 1] + parts = [tr[~tr["cik"].isin(multi.index)]] + for c, ts in multi.items(): + sub = tr[tr["cik"] == c] + parts += [sub.assign(ticker=t) for t in ts] + return pd.concat(parts, ignore_index=True) + + +def share_class_group(tm: pd.DataFrame) -> dict[str, str]: + alias = tm[tm["source"] == "sec_map"].groupby("cik")["ticker"].apply(sorted) + return {t: ts[0] for ts in alias for t in ts} + + +def covered(tm: pd.DataFrame, dates: pd.DatetimeIndex, tickers: list[str]) -> pd.DataFrame: + """True where an issuer CIK mapped to the ticker was filing Form 4s: from + the session after its first filing in the data to COVER_SLACK days after + its last one.""" + col = {t: i for i, t in enumerate(tickers)} + out = np.zeros((len(dates), len(tickers)), bool) + for r in tm.itertuples(index=False): + if r.ticker in col: + a = dates.searchsorted(r.first, side="right") + b = dates.searchsorted(r.last + pd.Timedelta(days=COVER_SLACK), side="right") + out[a:b, col[r.ticker]] = True + return pd.DataFrame(out, index=dates, columns=tickers) + + +def window(ev: pd.DataFrame, dates: pd.DatetimeIndex, n_cols: int, days: int, + weights: np.ndarray | None = None) -> np.ndarray: + """(dates x columns) count of distinct ``insider_cik`` with a filing in the + last ``days`` calendar days before each date; with ``weights``, the sum of + the weights of those filings instead.""" + filed = ev["filed"].values + s = dates.searchsorted(filed, side="right") # first session after filing + e = dates.searchsorted(filed + np.timedelta64(days, "D"), side="right") + D = np.zeros((len(dates) + 1, n_cols)) + c = ev["col"].values + if weights is not None: + np.add.at(D, (s, c), weights) + np.add.at(D, (e, c), -weights) + out = np.cumsum(D, axis=0)[:-1] + out[np.abs(out) < 1.0] = 0.0 # float residue of the running sum + return out + # merge each insider's overlapping [s, e) spans so she counts once + df = pd.DataFrame({"c": c, "k": ev["insider_cik"].values, "f": filed, "s": s, "e": e}) + df = df.sort_values(["c", "k", "f"]) + new = ((df["c"] != df["c"].shift()) | (df["k"] != df["k"].shift()) + | (df["s"] > df["e"].shift())).values + run = np.cumsum(new) + g = df.groupby(run).agg(c=("c", "first"), s=("s", "min"), e=("e", "max")) + np.add.at(D, (g["s"].values, g["c"].values), 1.0) + np.add.at(D, (g["e"].values, g["c"].values), -1.0) + return np.cumsum(D, axis=0)[:-1] + + +# --------------------------------------------------------------------------- +# build +# --------------------------------------------------------------------------- + +def build(p) -> dict[str, pd.DataFrame]: + t0 = time.time() + tr, tm = _raw() + dates, tickers = p.dates, list(p.close.columns) + tr = tr.assign(cmp=classify_cmp(tr)) + ex = _explode(tr, tm) + col = pd.Series(np.arange(len(tickers)), index=tickers) + ex = ex[ex["ticker"].isin(col.index)].copy() + ex["col"] = col.reindex(ex["ticker"]).values + cov = covered(tm, dates, tickers) + n = len(tickers) + buy, sell = ex[ex["code"] == "P"], ex[ex["code"] == "S"] + print(f" {len(ex):,} trades on {ex['ticker'].nunique()} tickers [{time.time() - t0:.0f}s]", flush=True) + + def frame(x): + return pd.DataFrame(x, index=dates, columns=tickers).where(cov).astype("float32") + + F = {} + F["insider_buyers_90d"] = frame(window(buy, dates, n, 90) - window(sell, dates, n, 90)) + signed = np.where(ex["code"] == "P", 1.0, -1.0) * ex["value"].fillna(0.0).values + net = window(ex, dates, n, 180, weights=signed) + mcap = fu.load(p)["mcap"].reindex(index=dates, columns=tickers) + F["insider_net_value_180d"] = frame(net / mcap.values) + ob = buy[buy["cmp"] == "opportunistic"] + osl = sell[sell["cmp"] == "opportunistic"] + opp = window(ob, dates, n, 90) - window(osl, dates, n, 90) + opp[dates < pd.Timestamp(OPP_START)] = np.nan + F["insider_opportunistic"] = frame(opp) + F["insider_cluster_buy"] = frame((window(buy, dates, n, 30) >= 2).astype(float)) + print(f" factors built [{time.time() - t0:.0f}s]", flush=True) + return F + + +def load(p, rebuild: bool = False) -> dict[str, pd.DataFrame]: + if os.path.exists(CACHE) and not rebuild: + with open(CACHE, "rb") as fh: + return pickle.load(fh) + F = build(p) + with open(CACHE, "wb") as fh: + pickle.dump(F, fh, protocol=4) + return F + + +# --------------------------------------------------------------------------- +# evaluation +# --------------------------------------------------------------------------- + +END = "2026-09-24" +# 2010-2016 = the first period without the 2008-09 crash and rebound +PERIODS = {"2006-2016": ("2006-01-01", "2016-12-31"), "2010-2016": ("2010-01-01", "2016-12-31"), + "2017-2026": ("2017-01-01", END), "all": ("2006-01-01", END)} +HORIZONS = (5, 21, 63) +# first date each portfolio is evaluated: full 180-day window / classifiable +# history / market caps for >= ~80% of members +START = {"insider_buyers_90d": "2006-07-03", "insider_cluster_buy": "2006-07-03", + "insider_opportunistic": OPP_START, "insider_net_value_180d": "2010-07-01"} +LONG_START = "2006-07-03" + + +def _avail_rows(dates: pd.DatetimeIndex, filed: pd.Series) -> np.ndarray: + return dates.searchsorted(filed.values, side="right") + + +def trade_events(p, ex: pd.DataFrame, group: dict[str, str]) -> pd.DataFrame: + """Ticker-days on which a filing became usable, for PIT members only; + share classes of one issuer count once.""" + i = _avail_rows(p.dates, ex["filed"]) + e = ex.assign(row=i)[i < len(p.dates)] + e = e[p.member.values[e["row"].values, e["col"].values]] + e = e.assign(grp=e["ticker"].map(lambda t: group.get(t, t))) + return e.drop_duplicates(["grp", "row"])[["row", "col"]] + + +def onset_events(p, flag: pd.DataFrame, group: dict[str, str]) -> pd.DataFrame: + """Days a 0/1 panel switches from 0 to 1, for PIT members.""" + v = flag.values + prev = np.vstack([np.full((1, v.shape[1]), np.nan), v[:-1]]) + on = (v == 1) & (prev == 0) & p.member.values + r, c = np.nonzero(on) + e = pd.DataFrame({"row": r, "col": c}) + e["grp"] = [group.get(t, t) for t in flag.columns[c]] + return e.drop_duplicates(["grp", "row"])[["row", "col"]] + + +def excess_fwd(p, h: int) -> np.ndarray: + """Open(t) → open(t+h) return minus the equal-weight PIT-member average.""" + o = p.open.where(p.open > 0).ffill() + f = o.shift(-h) / o - 1 + ew = f.where(p.member).mean(axis=1) + return f.sub(ew, axis=0).values + + +def clustered_t(x: np.ndarray, month: np.ndarray) -> tuple[float, float, int]: + ok = ~np.isnan(x) + x, month = x[ok], month[ok] + n = len(x) + if n < 10: + return float("nan"), float("nan"), n + m = x.mean() + g = pd.Series(x - m).groupby(month).sum() + G = len(g) + se = np.sqrt((g ** 2).sum() * G / max(G - 1, 1)) / n + return float(m), float(m / se), n + + +def event_study(p, events: dict[str, pd.DataFrame]) -> dict: + X = {h: excess_fwd(p, h) for h in HORIZONS} + ym = np.asarray(p.dates.year * 100 + p.dates.month) + out = {} + for kind, e in events.items(): + out[kind] = {} + for per, (a, b) in PERIODS.items(): + ra, rb = p.dates.searchsorted(pd.Timestamp(a)), p.dates.searchsorted(pd.Timestamp(b), side="right") + ee = e[(e["row"] >= ra) & (e["row"] < rb)] + out[kind][per] = {} + for h in HORIZONS: + x = X[h][ee["row"].values, ee["col"].values] + m, t, n = clustered_t(x, ym[ee["row"].values]) + out[kind][per][f"+{h}"] = {"mean": m, "t": t, "n": n} + return out + + +def coverage(p, I: dict, ex: pd.DataFrame) -> pd.DataFrame: + yr = p.dates.year + mem = p.member + cov = I["insider_buyers_90d"].notna() & mem + i = _avail_rows(p.dates, ex["filed"]) + e = ex.assign(row=i)[i < len(p.dates)] + e = e[mem.values[e["row"].values, e["col"].values]] + e = e.assign(year=p.dates[e["row"].values].year).drop_duplicates(["accession", "code"]) + return pd.DataFrame({ + "members": mem.sum(axis=1).groupby(yr).mean().round(0), + "mapped": cov.sum(axis=1).groupby(yr).mean().round(0), + "buy_filings": e[e["code"] == "P"].groupby("year").size(), + "sell_filings": e[e["code"] == "S"].groupby("year").size(), + "opp_buy_filings": e[(e["code"] == "P") & (e["cmp"] == "opportunistic")].groupby("year").size(), + "names_bought": e[e["code"] == "P"].groupby("year")["ticker"].nunique(), + }).loc[2006:].fillna(0).astype(int) + + +def _active(r: pd.Series, bench: pd.Series) -> tuple[float, float]: + """Annualized mean daily return over ``bench`` and its t-stat.""" + act = (r - bench.reindex(r.index)).dropna() + ir = act.mean() / act.std() * np.sqrt(252) + return act.mean() * 252, ir * np.sqrt(len(act) / 252) + + +def _period_stats(r: pd.Series, ew: pd.Series, tie: pd.Series, a: str, b: str) -> dict: + r = r.loc[a:b].dropna() + s = stats((1 + r).cumprod()) + act, t = _active(r, ew) + act_tie, t_tie = _active(r, tie) + return {"sharpe": s["sharpe"], "cagr": s["cagr"], "mdd": s["mdd"], + "active_ann": act, "active_t": t, "vs_tiebreak_ann": act_tie, "vs_tiebreak_t": t_tie} + + +def main() -> None: + from research.lab.rotation import factor_ranks + from research.lab.us_sweep import _blend, run + + p, _ = load_panel("us") + I = load(p) + tr, tm = _raw() + tr = tr.assign(cmp=classify_cmp(tr), cmp_dir=classify_cmp(tr, same_direction=True)) + ex = _explode(tr, tm) + col = pd.Series(np.arange(len(p.close.columns)), index=p.close.columns) + ex = ex[ex["ticker"].isin(col.index)].copy() + ex["col"] = col.reindex(ex["ticker"]).values + group = share_class_group(tm) + out: dict = {} + + # ---- mapping and coverage + src = tm.groupby("source")["ticker"].nunique().to_dict() + unmapped = sorted(set(p.close.columns) - set(tm["ticker"])) + unm_mem = [t for t in unmapped if p.member[t].any()] + print(f"ticker map: {src} ; unmapped ever-members: {unm_mem}") + print("by symbol:", sorted(tm.loc[tm["source"] == "symbol", "ticker"])) + cv = coverage(p, I, ex) + print("\n=== coverage among PIT S&P 500 members (per year; filings = Form 4s with a P or S line)") + print(cv.to_string()) + out["mapping"] = {"by_source": src, "unmapped_ever_members": unm_mem, + "by_symbol": sorted(tm.loc[tm["source"] == "symbol", "ticker"])} + out["coverage"] = {int(k): {c: int(v) for c, v in r.items()} for k, r in cv.iterrows()} + + # ---- event study + B, S = ex["code"] == "P", ex["code"] == "S" + ev = { + "all purchases": trade_events(p, ex[B], group), + "opportunistic purchases": trade_events(p, ex[B & (ex["cmp"] == "opportunistic")], group), + "cluster buys (onset)": onset_events(p, I["insider_cluster_buy"], group), + "opportunistic sales": trade_events(p, ex[S & (ex["cmp"] == "opportunistic")], group), + # references + "routine purchases": trade_events(p, ex[B & (ex["cmp"] == "routine")], group), + "all sales": trade_events(p, ex[S], group), + "opp. purchases, same-direction rule": trade_events(p, ex[B & (ex["cmp_dir"] == "opportunistic")], group), + "opp. sales, same-direction rule": trade_events(p, ex[S & (ex["cmp_dir"] == "opportunistic")], group), + } + es = event_study(p, ev) + out["event_study"] = es + print("\n=== event study: excess return vs EW members from the next open after filing, % " + "(t clustered by month) [n]") + print(f"{'':38s}" + "".join(f"{per:>27s}" for per in PERIODS)) + for kind, d in es.items(): + for h in HORIZONS: + cells = "".join(f" {d[per][f'+{h}']['mean']*100:+6.2f} ({d[per][f'+{h}']['t']:+5.1f}) " + f"[{d[per][f'+{h}']['n']:5d}]" for per in PERIODS) + print(f"{(kind if h == 5 else ''):34s} +{h:<3d}{cells}") + + # ---- portfolios: top 20 among members, factor first, 12-1 momentum breaks ties + pr = factor_ranks(p, "us", ["mom_12_1", "overnight_mom", "rec_126"]) + tie = pr["mom_12_1"].fillna(0.0) + + def lexi(f): + r = f.where(p.member).rank(axis=1, method="dense") + return (r + 0.5 * tie).where(r.notna()) + + kw = dict(top_n=20, buffer=2.0, every=21) + spans = lambda a: [(a, "2016-12-30"), ("2017-01-03", END)] # noqa: E731 + R, runs = {}, {} + ew_sc = p.close.notna().astype(float).where(p.member) + runs["ew_universe"] = run(p, ew_sc, spans(LONG_START), top_n=10000, every=21) + runs["mom_12_1 alone (tie-breaker)"] = run(p, pr["mom_12_1"].where(p.member), spans(LONG_START), **kw) + mom = _blend(pr["overnight_mom"], pr["rec_126"]) + runs["momentum sleeve"] = run(p, mom.where(p.member), spans(LONG_START), top_n=10, buffer=2.0, every=10) + for name in NAMES: + runs[name] = run(p, lexi(I[name]), spans(START[name]), **kw) + print(f" ran {name}", flush=True) + for k, v in runs.items(): + R[k] = v["_rets"] + ew, ms, tb = R["ew_universe"], R["momentum sleeve"], R["mom_12_1 alone (tie-breaker)"] + + print("\n=== long-only top 20 (every 21 days, buffer 2, next-open fills, 10 bp); " + "periods: start-2016 | 2017-2026 | whole") + print(" sharpe / cagr / mdd per period; active = vs EW universe, vs-tie = vs the 12-1 momentum " + "tie-breaker alone (ann. %, t); phase-avg block Sharpe; corr with momentum sleeve (raw, active vs EW)") + port = {} + for k, v in runs.items(): + start = START.get(k, LONG_START) + pers = {"p1": (start, "2016-12-30"), "p2": ("2017-01-03", END), "all": (start, END)} + st = {pk: _period_stats(R[k], ew, tb, a, b) for pk, (a, b) in pers.items()} + r = R[k].loc[start:END] + c_raw = r.corr(ms.loc[start:END]) + c_act = (r - ew.loc[start:END]).corr(ms.loc[start:END] - ew.loc[start:END]) + n_sig = None + if k in I: + f = I[k].where(p.member).loc[start:END] + n_sig = float((f > 0).sum(axis=1).mean()) + port[k] = {"start": start, **st, "blocks": v["blocks"], "phase_sharpe": v["sharpe"], + "phase_cagr": v["cagr"], "phase_mdd": v["mdd"], "turn": v["turn"], + "corr_mom": c_raw, "corr_mom_active": c_act, "avg_signal_names": n_sig} + cells = " | ".join(f"{st[pk]['sharpe']:.2f} {st[pk]['cagr']*100:5.1f}% {st[pk]['mdd']*100:5.1f}%" + for pk in ("p1", "p2", "all")) + act = " | ".join(f"{st[pk]['active_ann']*100:+5.1f}% ({st[pk]['active_t']:+4.1f})" for pk in ("p1", "p2", "all")) + vst = " | ".join(f"{st[pk]['vs_tiebreak_ann']*100:+5.1f}% ({st[pk]['vs_tiebreak_t']:+4.1f})" + for pk in ("p1", "p2", "all")) + sig = f" sig-names {n_sig:5.1f}" if n_sig is not None else "" + print(f"{k:30s} from {start} {cells}\n{'':30s} active {act} blocks " + f"[{' '.join(f'{b:.2f}' for b in v['blocks'])}] corr {c_raw:.2f}/{c_act:+.2f} " + f"turn {v['turn']:.1f}x{sig}") + if k in I: + print(f"{'':30s} vs-tie {vst}") + out["portfolio"] = port + os.makedirs(RESULTS, exist_ok=True) + path = os.path.join(RESULTS, "insider_eval.json") + with open(path, "w") as fh: + json.dump(out, fh, indent=1, default=float) + print(f"\nSaved {path}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/lazy_prices.py b/research/lab/lazy_prices.py new file mode 100644 index 0000000..f794239 --- /dev/null +++ b/research/lab/lazy_prices.py @@ -0,0 +1,336 @@ +""" +"Lazy Prices" factors for the US panel (Cohen, Malloy & Nguyen 2020, JF): +firms whose 10-K / 10-Q text barely changes from the comparable prior filing +("non-changers") outperformed "changers" in 1995-2014. + +Source: data/research/lazy_prices.csv from research/lab/fetch_lazy_prices.py +(SEC EDGAR primary documents; each 10-K vs the previous 10-K, each 10-Q vs the +10-Q of the same fiscal quarter a year earlier). + +Point-in-time rule: a filing's similarity is usable from the first trading +day strictly after its EDGAR filing date, and the most recent filing's value +is held until the next filing (a newer filing without a comparable prior one +resets it to NaN). Values older than 400 days are dropped. Filings whose own +or reference text has fewer than 500 words (exhibit-only shells, cover pages +of documents incorporated by reference) get no value. + + from research.lab import lazy_prices as lp + F = lp.load(panel) # {name: dates x tickers}, cached + # lazy_cos cosine similarity of term counts, latest 10-K or 10-Q (higher = fewer changes) + # lazy_jac Jaccard similarity of word sets, latest 10-K or 10-Q + # lazy_cos_10k cosine similarity, latest 10-K only + # lazy_rf cosine similarity of Item 1A Risk Factors, latest 10-K only + # lazy_cos_pct lazy_cos as a percentile among same-form filings of the past + # 365 days (10-Qs are less similar than 10-Ks: median cosine + # 0.97 vs 0.99, so the raw mix partly sorts on form type) + # lazy_jac_pct the same for lazy_jac (median 0.78 for 10-Qs vs 0.88) + +Evaluation (coverage, monthly quintile sorts and rank IC, the same sorts on +within-sector ranks, long-only top 20, correlation with the momentum sleeve) +for 2006-2014, 2015-2026 and the whole span; writes data/research/results/lazy_prices_eval.json: + + PYTHONPATH=. uv run python -m research.lab.lazy_prices +""" + +from __future__ import annotations + +import json +import os +import pickle +import warnings + +import numpy as np +import pandas as pd +from scipy import stats as sps + +from research.lab.fetch_lazy_prices import GROUPS, OUT as RAW +from research.lab.report import stats + +warnings.filterwarnings("ignore") + +CACHE = os.path.join("data", "research", "lazy_prices_factors.pkl") +MAX_AGE = 400 +MIN_WORDS = 500 +FACTORS = { # name: (similarity column, forms) + "lazy_cos": ("cos_sim", GROUPS["10-K"] + GROUPS["10-Q"]), + "lazy_jac": ("jac_sim", GROUPS["10-K"] + GROUPS["10-Q"]), + "lazy_cos_10k": ("cos_sim", GROUPS["10-K"]), + "lazy_rf": ("cos_rf", GROUPS["10-K"]), + "lazy_cos_pct": ("cos_pct", GROUPS["10-K"] + GROUPS["10-Q"]), + "lazy_jac_pct": ("jac_pct", GROUPS["10-K"] + GROUPS["10-Q"]), +} +PERIODS = {"2006-2014": ("2006-01-01", "2014-12-31"), "2015-2026": ("2015-01-01", "2026-09-24"), + "2006-2026": ("2006-01-01", "2026-09-24")} +BLOCKS = [("2006-01-03", "2010-12-31"), ("2011-01-01", "2014-12-31"), ("2015-01-01", "2018-12-31"), + ("2019-01-01", "2022-12-31"), ("2023-01-01", "2026-09-24")] + + +# --------------------------------------------------------------------------- +# factor panels +# --------------------------------------------------------------------------- + +def filings() -> pd.DataFrame: + """One row per (ticker, filing) with the word-count screen applied.""" + df = pd.read_csv(RAW, dtype={"period": str, "ref_accession": str}) + df["filed"] = pd.to_datetime(df["filed"]) + words = df.drop_duplicates("accession").set_index("accession")["n_words"] + ref_words = df["ref_accession"].map(words) + short = (df["n_words"] < MIN_WORDS) | (ref_words < MIN_WORDS) + df.loc[short, ["cos_sim", "jac_sim", "cos_rf", "jac_rf"]] = np.nan + df["cos_pct"] = _trailing_pct(df, "cos_sim") + df["jac_pct"] = _trailing_pct(df, "jac_sim") + return df + + +def _trailing_pct(df: pd.DataFrame, col: str, days: int = 365) -> pd.Series: + """A similarity as the share of same-form filings filed in the past + ``days`` days (through the same filing date) with a value at or below it.""" + out = pd.Series(np.nan, index=df.index) + u = df.drop_duplicates("accession").dropna(subset=[col]) + pct = {} + for forms in GROUPS.values(): + g = u[u["form"].isin(forms)].sort_values("filed") + f = g["filed"].values + v = g[col].values + lo = np.searchsorted(f, f - np.timedelta64(days, "D"), side="left") + hi = np.searchsorted(f, f, side="right") + for acc, a, b, x in zip(g["accession"].values, lo, hi, v): + pct[acc] = (v[a:b] <= x).mean() + ok = df[col].notna() + out[ok] = df.loc[ok, "accession"].map(pct) + return out + + +def _asof(ev: pd.DataFrame, col: str, dates: pd.DatetimeIndex, tickers: list[str]) -> pd.DataFrame: + """Latest filing's value from the session after its filing date; NaN + values reset, and values older than MAX_AGE days drop out.""" + out = np.full((len(dates), len(tickers)), np.nan, dtype="float32") + cols = {t: i for i, t in enumerate(tickers)} + day = dates.values.astype("datetime64[D]").astype("int64") + for t, g in ev.groupby("ticker"): + if t not in cols: + continue + g = g.sort_values("filed") + pos = dates.searchsorted(g["filed"].values, side="right") + ok = pos < len(dates) + pos, vals = pos[ok], g[col].values[ok] + fday = g["filed"].values[ok].astype("datetime64[D]").astype("int64") + if not len(pos): + continue + # several filings before one session: the last one wins + last = np.r_[pos[1:] != pos[:-1], True] + pos, vals, fday = pos[last], vals[last], fday[last] + k = np.full(len(dates), -1) + k[pos] = np.arange(len(pos)) + k = np.maximum.accumulate(k) + have = k >= 0 + v = np.full(len(dates), np.nan) + v[have] = vals[k[have]] + v[have & (day - np.where(have, fday[np.maximum(k, 0)], 0) > MAX_AGE)] = np.nan + out[:, cols[t]] = v + return pd.DataFrame(out, index=dates, columns=tickers) + + +def build(p) -> dict[str, pd.DataFrame]: + df = filings() + dates, tickers = p.dates, list(p.close.columns) + return {name: _asof(df[df["form"].isin(forms)], col, dates, tickers) + for name, (col, forms) in FACTORS.items()} + + +def load(p, rebuild: bool = False) -> dict[str, pd.DataFrame]: + if os.path.exists(CACHE) and not rebuild: + with open(CACHE, "rb") as fh: + return pickle.load(fh) + F = build(p) + with open(CACHE, "wb") as fh: + pickle.dump(F, fh, protocol=4) + return F + + +# --------------------------------------------------------------------------- +# evaluation +# --------------------------------------------------------------------------- + +def nw_t(x, lags: int | None = None) -> float: + """Newey-West t-statistic of the mean (Bartlett kernel).""" + x = np.asarray(pd.Series(x).dropna(), dtype=float) + T = len(x) + if T < 12: + return np.nan + L = lags if lags is not None else int(4 * (T / 100) ** (2 / 9)) + e = x - x.mean() + s = e @ e / T + for lag in range(1, L + 1): + s += 2 * (1 - lag / (L + 1)) * (e[lag:] @ e[:-lag]) / T + return float(x.mean() / np.sqrt(s / T)) + + +def monthly_returns(p) -> tuple[pd.DatetimeIndex, pd.DataFrame]: + """Signal on each month's last session, filled at the next open, held to + the open after the next month's last session.""" + o = p.open.where(p.open > 0).ffill() + me = pd.DatetimeIndex(p.dates.to_series().groupby(p.dates.to_period("M")).max().values) + pos = p.dates.searchsorted(me, side="right") + ok = pos < len(p.dates) + me, pos = me[ok], pos[ok] + fill = o.iloc[pos].values + R = pd.DataFrame(fill[1:] / fill[:-1] - 1, index=me[:-1], columns=p.close.columns) + return me[:-1], R + + +def quintiles(score: pd.DataFrame, R: pd.DataFrame, member: pd.DataFrame, + min_names: int = 50) -> pd.DataFrame: + """Per month: EW return of each quintile minus the EW member average, + Q5-Q1, and the rank IC (Q1 = lowest similarity = most changes).""" + rows = {} + for d in R.index: + m = member.loc[d] + r = R.loc[d].where(m) + univ = r.mean() + s = score.loc[d].where(m & r.notna()).dropna() + if len(s) < min_names: + continue + q = np.floor((s.rank(method="first") - 1) / len(s) * 5).astype(int) + row = {f"Q{k + 1}": r[q.index[q == k]].mean() - univ for k in range(5)} + row["Q5-Q1"] = row["Q5"] - row["Q1"] + row["ic"] = sps.spearmanr(s.values, r[s.index].values).statistic + row["n"] = len(s) + rows[d] = row + cols = ["Q1", "Q2", "Q3", "Q4", "Q5", "Q5-Q1", "ic", "n"] + return pd.DataFrame(rows, index=cols).T.astype(float) + + +def _qtable(Q: pd.DataFrame) -> dict: + out = {} + for c in ["Q1", "Q2", "Q3", "Q4", "Q5", "Q5-Q1"]: + out[c] = {"mean_bp": float(Q[c].mean() * 1e4), "t_nw": nw_t(Q[c]), + "t": float(Q[c].mean() / Q[c].std() * np.sqrt(Q[c].count()))} + out["ic"] = {"mean": float(Q["ic"].mean()), "t_nw": nw_t(Q["ic"]), + "t": float(Q["ic"].mean() / Q["ic"].std() * np.sqrt(Q["ic"].count())), + "hit": float((Q["ic"] > 0).mean())} + out["months"] = int(len(Q)) + out["avg_names"] = float(Q["n"].mean()) + return out + + +def coverage(F: pd.DataFrame, member: pd.DataFrame) -> dict: + have = (F.notna() & member).sum(axis=1) + mem = member.sum(axis=1) + y = pd.DataFrame({"have": have, "mem": mem}).loc["2005":].groupby(lambda d: d.year).mean() + return {int(k): {"members": round(float(r["mem"]), 1), "with_value": round(float(r["have"]), 1), + "share": round(float(r["have"] / r["mem"]), 3)} for k, r in y.iterrows()} + + +def _monthly(r: pd.Series) -> pd.Series: + return (1 + r).resample("ME").prod() - 1 + + +def main() -> None: + from research.lab.factor_scan import RESULTS, load as load_panel + from research.lab.rotation import factor_ranks + from research.lab.us_final import describe + from research.lab.us_sweep import _blend, run, sector_rank + + p, _ = load_panel("us") + F = load(p, rebuild=True) + df = filings() + out: dict = {"n_filings": {}, "similarity_by_form": {}} + for g, forms in GROUPS.items(): + sub = df[df["form"].isin(forms)].drop_duplicates("accession") + out["n_filings"][g] = {"total": int(len(sub)), "with_cos": int(sub["cos_sim"].notna().sum()), + "with_rf": int(sub["cos_rf"].notna().sum())} + by_y = sub.groupby(sub["filed"].dt.year)[["cos_sim", "jac_sim", "cos_rf"]].median() + out["similarity_by_form"][g] = {int(y): {k: round(float(v), 4) for k, v in r.items()} + for y, r in by_y.iterrows()} + print("filings:", out["n_filings"]) + + # coverage + out["coverage"] = {n: coverage(f, p.member) for n, f in F.items()} + print("\n=== coverage: PIT members with a value (yearly mean of daily counts)") + print("year members " + " ".join(f"{n:>13s}" for n in F)) + for y in out["coverage"]["lazy_cos"]: + c = out["coverage"] + print(f"{y} {c['lazy_cos'][y]['members']:7.0f} " + + " ".join(f"{c[n][y]['with_value']:6.0f} ({c[n][y]['share']:4.0%})" for n in F)) + + # quintile sorts and IC + _, R = monthly_returns(p) + out["quintiles"] = {} + print("\n=== monthly quintiles, EW among PIT members: mean excess over member EW, bp/month " + "[NW t]; Q1 = most changes, Q5 = fewest") + for n, f in F.items(): + Q = quintiles(f, R, p.member) + out["quintiles"][n] = {} + for per, (a, b) in PERIODS.items(): + tb = _qtable(Q.loc[a:b]) + out["quintiles"][n][per] = tb + print(f"{n:13s} {per}: " + " ".join( + f"{c} {tb[c]['mean_bp']:+5.1f} [{tb[c]['t_nw']:+.1f}]" for c in ["Q1", "Q2", "Q3", "Q4", "Q5"]) + + f" | Q5-Q1 {tb['Q5-Q1']['mean_bp']:+5.1f} [{tb['Q5-Q1']['t_nw']:+.2f}]" + + f" IC {tb['ic']['mean']:+.4f} [{tb['ic']['t_nw']:+.2f}] n={tb['avg_names']:.0f}") + yr = Q.groupby(Q.index.year)[["Q5-Q1", "ic"]].mean() + out["quintiles"][n]["by_year"] = {int(y): {"q5_q1_bp": float(r["Q5-Q1"] * 1e4), "ic": float(r["ic"])} + for y, r in yr.iterrows()} + + # robustness: the same sorts on within-sector ranks (GICS, PIT members) + out["quintiles_sector_neutral"] = {} + print("\n=== sector-neutral ranks: Q5-Q1 bp/month [NW t], IC [NW t]") + for n, f in F.items(): + Q = quintiles(sector_rank(f, p), R, p.member) + out["quintiles_sector_neutral"][n] = {per: _qtable(Q.loc[a:b]) for per, (a, b) in PERIODS.items()} + print(f"{n:13s} " + " ".join( + f"{per}: {tb['Q5-Q1']['mean_bp']:+5.1f} [{tb['Q5-Q1']['t_nw']:+.2f}] IC {tb['ic']['mean']:+.4f} " + f"[{tb['ic']['t_nw']:+.2f}]" for per, tb in out["quintiles_sector_neutral"][n].items())) + + # long-only top 20 vs the equal-weight universe + span = [(BLOCKS[0][0], BLOCKS[-1][1])] + ew = p.close.notna().astype(float).where(p.member) + runs = {n: run(p, f.where(p.member), BLOCKS, top_n=20, buffer=2.0, every=21) for n, f in F.items()} + runs["ew_members"] = run(p, ew, BLOCKS, top_n=10000, every=21) + pr = factor_ranks(p, "us", ["overnight_mom", "rec_126"]) + mom = _blend(pr["overnight_mom"], pr["rec_126"]) + runs["momentum_sleeve"] = run(p, mom.where(p.member), span, top_n=10, buffer=2.0, every=10) + ew_r = runs["ew_members"]["_rets"] + out["top20"] = {} + print("\n=== long-only top 20 (every 21, buffer 2, next-open fills, 10 bp), EW members for reference") + print(f"{'':16s} " + " ".join(f"{per:^34s}" for per in PERIODS) + " blocks (sharpe)") + for n, st in runs.items(): + r = st["_rets"] + row = {"turn": float(st["turn"]), + "block_sharpe": {f"{a[:4]}-{b[:4]}": float(v) for (a, b), v in zip(BLOCKS, st["blocks"])} + if n != "momentum_sleeve" else {}} + cells = [] + for per, (a, b) in PERIODS.items(): + d = describe(r.loc[a:b], blocks=[]) + xs = (r - ew_r).loc[a:b] + ir = float(xs.mean() / xs.std() * np.sqrt(252)) if xs.std() > 0 else np.nan + row[per] = {"cagr": d["cagr"], "sharpe": d["sharpe"], "mdd": d["mdd"], "vol": d["vol"], + "excess_ann": float(xs.mean() * 252), "ir_vs_ew": ir} + cells.append(f"cagr {d['cagr']*100:5.1f}% sh {d['sharpe']:.2f} IR {ir:+.2f} ") + out["top20"][n] = row + bl = " ".join(f"{v:.2f}" for v in row["block_sharpe"].values()) + print(f"{n:16s} " + " ".join(cells) + f" [{bl}] turn {st['turn']:.1f}x") + + # correlation with the momentum sleeve (monthly) + M = pd.DataFrame({n: _monthly(st["_rets"]) for n, st in runs.items()}) + X = M.sub(M["ew_members"], axis=0) + out["corr_with_momentum"] = {} + print("\n=== monthly correlation with the momentum sleeve (raw / excess over EW members)") + for n in F: + row = {} + for per, (a, b) in PERIODS.items(): + row[per] = {"raw": float(M.loc[a:b, n].corr(M.loc[a:b, "momentum_sleeve"])), + "excess": float(X.loc[a:b, n].corr(X.loc[a:b, "momentum_sleeve"]))} + out["corr_with_momentum"][n] = row + print(f"{n:13s} " + " ".join(f"{per}: {v['raw']:+.2f} / {v['excess']:+.2f}" for per, v in row.items())) + + os.makedirs(RESULTS, exist_ok=True) + path = os.path.join(RESULTS, "lazy_prices_eval.json") + with open(path, "w") as fh: + json.dump(out, fh, indent=1, default=float) + print(f"\nsaved {path}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/letf_history.py b/research/lab/letf_history.py new file mode 100644 index 0000000..8be875e --- /dev/null +++ b/research/lab/letf_history.py @@ -0,0 +1,185 @@ +""" +Synthetic leveraged ETFs back to 1986, and the trend strategies run on them. + +The live ETF panel starts in 2016 and the TrendRider research in 2010, so +neither saw the 2000-2002 Nasdaq crash or 2008 with 3x ETFs. This rebuilds +the TrendRider basket from its underlying indices and replays it. + + PYTHONPATH=. uv run python -m research.lab.fetch_long_history # once + PYTHONPATH=. uv run python -m research.lab.letf_history + +Synthetic series (daily): + SPY S&P 500 total return: SPY (1993-), VFINX before that + QQQ Nasdaq-100 total return: QQQ (1999-03-), ^NDX (or FRED NASDAQ100) + + 0.5%/yr dividends before + UPRO / TQQQ 3 x underlying - 2 x (T-bill + 0.4%) - 0.95%/yr expense, daily reset + SHY 2-year Treasury from DGS2 (duration 1.9); TLT from DGS30 (duration 17; + DGS10 changes fill the 2002-2006 gap in DGS30) + GLD gold futures GC=F from 2000-08 (- 0.4%/yr); cash (the SHY series) before + DBC DBC from 2006-02; cash before (TrendRider's risk-off pick then has no + commodity leg) + +The synthetic 3x series are checked against the real TQQQ/UPRO. Strategies +are run in two fill modes that bracket the live book: + close weights from closes through t-1 traded at the close of t-1 (optimistic) + lag the same weights traded one session later (pessimistic; live trades + the same session on the open's information) +Book: $10k, 5 bp, $2 per order, 10% resize band (as trend_rider_long.py). +""" + +from __future__ import annotations + +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from main import simulate_live_book, strategy_positions +from strategies.base import Strategy +from trader import RESIZE_BAND, STRATEGY_REGISTRY + +warnings.filterwarnings("ignore") +SRC = os.path.join("data", "research", "long_history.pkl") +OUT = os.path.join("data", "research", "letf_history.pkl") +FIN_SPREAD, EXPENSE = 0.004, 0.0095 +BLOCKS = [("1990-01-01", "1999-12-31"), ("2000-01-01", "2002-12-31"), ("2003-01-01", "2007-12-31"), + ("2008-01-01", "2009-12-31"), ("2010-01-01", "2015-12-31"), ("2016-01-01", "2026-09-25"), + ("1990-01-01", "2026-09-25")] + + +def _splice(primary: pd.Series, backfill_ret: pd.Series, index: pd.DatetimeIndex) -> pd.Series: + """Daily returns: the primary series where it exists, backfill before it.""" + rp = primary.reindex(index).pct_change(fill_method=None) + rb = backfill_ret.reindex(index) + first = primary.first_valid_index() + return rb.where(index <= first, rp).fillna(0.0) if first is not None else rb.fillna(0.0) + + +def build(src: dict) -> tuple[pd.DataFrame, pd.DataFrame]: + """(synthetic closes, real closes for validation) on the S&P 500 calendar from 1986.""" + c, fred = src["close"], src["fred"].ffill() + idx = c["^GSPC"].dropna().loc["1986-01-01":].index + rf = fred["DTB3"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100 + y2 = fred["DGS2"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100 + y30 = fred["DGS30"].reindex(idx.union(fred.index)).ffill(limit=5).reindex(idx) / 100 + y10 = fred["DGS10"].reindex(idx.union(fred.index)).ffill().reindex(idx) / 100 + dy30 = y30.diff().where(y30.notna() & y30.shift(1).notna(), y10.diff()).fillna(0.0) + y30f = y30.fillna(y10) + + spx_bf = c["VFINX"].reindex(idx).pct_change(fill_method=None) + spx_bf = spx_bf.fillna(c["^GSPC"].reindex(idx).pct_change(fill_method=None) + 0.03 / 252) + spx = _splice(c["SPY"], spx_bf, idx) + ndx_px = c["^NDX"] if "^NDX" in c else fred["NASDAQ100"] + ndx_bf = ndx_px.reindex(idx).pct_change(fill_method=None) + 0.005 / 252 + ndx = _splice(c["QQQ"], ndx_bf, idx) + + def letf(r: pd.Series, lev: float) -> pd.Series: + return lev * r - (lev - 1) * (rf + FIN_SPREAD) / 252 - EXPENSE / 252 + + shy = y2.shift(1) / 252 - 1.9 * y2.diff().fillna(0.0) + tlt = y30f.shift(1) / 252 - 17.0 * dy30 + gold = c["GC=F"].reindex(idx).ffill().pct_change(fill_method=None) - 0.004 / 252 + gld = gold.where(idx > c["GC=F"].first_valid_index(), shy) + dbc = c["DBC"].reindex(idx).ffill().pct_change(fill_method=None) + dbc = dbc.where(idx > c["DBC"].first_valid_index(), shy) + rets = pd.DataFrame({"SPY": spx, "QQQ": ndx, "UPRO": letf(spx, 3), "TQQQ": letf(ndx, 3), + "SSO": letf(spx, 2), "QLD": letf(ndx, 2), "SHY": shy, "TLT": tlt, + "GLD": gld, "DBC": dbc}).fillna(0.0) + syn = 100 * (1 + rets).cumprod() + real = c.reindex(idx)[["SPY", "QQQ", "UPRO", "TQQQ", "SSO", "QLD", "SHY", "TLT", "GLD", "DBC"]] + return syn, real + + +def validate(syn: pd.DataFrame, real: pd.DataFrame) -> None: + print("synthetic vs real (overlap): daily corr, annual return real / synthetic") + for t in ["UPRO", "TQQQ", "SSO", "QLD", "SHY", "TLT", "GLD"]: + r = real[t].pct_change(fill_method=None) + s = syn[t].pct_change(fill_method=None) + ok = r.notna() & s.notna() + r, s = r[ok], s[ok] + yrs = len(r) / 252 + print(f" {t:5s} {r.index[0].date()}- corr {r.corr(s):.4f} real {((1 + r).prod() ** (1 / yrs) - 1) * 100:6.2f}%" + f" syn {((1 + s).prod() ** (1 / yrs) - 1) * 100:6.2f}%") + + +class HoldRotation(Strategy): + """Gayed & Bilello (2016) style: hold ``on`` while the signal is above its + ``ma``-day average, else ``off``.""" + + def __init__(self, on: str, off: str = "SHY", signal: str = "SPY", ma: int = 200): + self.on, self.off, self.signal, self.ma = on, off, signal, ma + + def generate_signals(self, data: pd.DataFrame) -> pd.DataFrame: + s = data[self.signal] + up = s > s.rolling(self.ma).mean() + w = pd.DataFrame(0.0, index=data.index, columns=data.columns) + w[self.on] = up.astype(float) + w[self.off] = (~up).astype(float) + w.iloc[:self.ma] = 0.0 + return w.shift(1).fillna(0.0) + + +class Hold(Strategy): + def __init__(self, asset: str): + self.asset = asset + + def generate_signals(self, data: pd.DataFrame) -> pd.DataFrame: + w = pd.DataFrame(0.0, index=data.index, columns=data.columns) + w[self.asset] = 1.0 + return w + + +def strategies() -> dict[str, Strategy]: + out = {k: STRATEGY_REGISTRY[k]() for k in ("trend_rider_v7_vt24", "trend_rider_v7_vt36")} + out.update({"ma200 TQQQ/SHY": HoldRotation("TQQQ"), "ma200 UPRO/SHY": HoldRotation("UPRO"), + "ma200 SSO/SHY": HoldRotation("SSO"), "hold SPY": Hold("SPY"), "hold QQQ": Hold("QQQ"), + "hold TQQQ": Hold("TQQQ"), "hold UPRO": Hold("UPRO")}) + return out + + +def run(strategy: Strategy, close: pd.DataFrame, mode: str) -> pd.Series: + held = strategy_positions(strategy, close) + if mode == "lag": + held = held.shift(1).fillna(0.0) + eq, _ = simulate_live_book(held, close, 10_000, 0.0005, 2.0, 0, 0, 50.0, RESIZE_BAND) + return eq + + +def block_stats(eq: pd.Series) -> list[str]: + r = eq.pct_change().fillna(0.0) + out = [] + for a, b in BLOCKS: + x = r.loc[a:b] + if len(x) < 100: + out.append(" n/a") + continue + e = (1 + x).cumprod() + cagr = e.iloc[-1] ** (252 / len(x)) - 1 + sh = x.mean() / x.std() * np.sqrt(252) if x.std() > 0 else 0.0 + out.append(f"{cagr*100:6.1f}%/{sh:4.2f}/{((e / e.cummax()) - 1).min()*100:4.0f}%") + return out + + +def main() -> None: + with open(SRC, "rb") as fh: + src = pickle.load(fh) + syn, real = build(src) + validate(syn, real) + close = syn.loc["1989-01-01":] + head = " ".join(f"{a[:4]}-{b[2:4]:>2s}".rjust(16) for a, b in BLOCKS) + print(f"\nCAGR/Sharpe/MaxDD by block (synthetic ETFs; GLD = cash before 2000-09, DBC = cash before 2006-02)") + print(f"{'':30s}{head}") + eqs = {} + for name, strat in strategies().items(): + for mode in ("close", "lag"): + eq = run(strat, close, mode) + eqs[f"{name} [{mode}]"] = eq + print(f"{name + ' [' + mode + ']':30s}" + " ".join(s.rjust(16) for s in block_stats(eq)), flush=True) + with open(OUT, "wb") as fh: + pickle.dump({"close": syn, "equity": eqs}, fh) + + +if __name__ == "__main__": + main() diff --git a/research/lab/levered_book.py b/research/lab/levered_book.py new file mode 100644 index 0000000..fe56af8 --- /dev/null +++ b/research/lab/levered_book.py @@ -0,0 +1,172 @@ +""" +How much return can risk-managed leverage add to the US sleeves, and at what +drawdown? A small account can borrow (Reg-T margin, at most 2x) where a fund +with a mandate cannot; the question is whether scaling exposure by recent +volatility lets that leverage be taken without the tail that static leverage +brings. + + PYTHONPATH=. uv run python -m research.lab.levered_book + +Input: daily sleeve returns from research/lab/us_final.py (and, when present, +short_horizon.py). Leverage is applied to the sleeve's return stream: + + exposure_t = clip(target / vol_hat_{t-1}, 0, cap) (vol_hat from returns through t-1) + r_t = exposure_t * sleeve_t - max(exposure_t - 1, 0) * (T-bill + spread) + - |exposure_t - exposure_{t-1}| * cost + +Exposure only moves when it is more than ``band`` away from the current one +(fewer orders). Stock sleeves may go to 2x (Reg-T); the trend sleeve already +holds 3x ETFs, which brokers margin at 75-100%, so it is capped at 1x. + +Vol estimators: realized over the last 126 days (Barroso & Santa-Clara 2015) +and EWMA with a 20-day half-life. Targets are round numbers, not fitted. +""" + +from __future__ import annotations + +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factor_scan import RESULTS +from research.lab.us_final import monthly_blend + +warnings.filterwarnings("ignore") + +SPAN = ("2011-01-03", "2026-09-24") +LONG = ("2006-09-27", "2026-09-24") +SPREAD = 0.015 # IBKR benchmark + 1.5% under $100k, +1% to $1M +COST = 0.001 # one-way, per unit of exposure changed +RF_PATH = os.path.join("data", "research", "ext", "fred_dtb3.csv") + + +def tbill(index: pd.DatetimeIndex) -> pd.Series: + """Daily 3-month T-bill rate (annual, decimal) from FRED DTB3.""" + rf = pd.read_csv(RF_PATH, index_col=0, parse_dates=True).iloc[:, 0] + rf = pd.to_numeric(rf, errors="coerce").ffill() / 100 + return rf.reindex(index.union(rf.index)).ffill().reindex(index) + + +def vol_hat(r: pd.Series, kind: str) -> pd.Series: + if kind == "rv126": + v = np.sqrt((r ** 2).rolling(126, min_periods=63).mean() * 252) + else: + v = np.sqrt((r ** 2).ewm(halflife=20, min_periods=63).mean() * 252) + return v.shift(1) + + +def banded(target: pd.Series, band: float) -> pd.Series: + out = np.empty(len(target)) + cur = np.nan + for i, t in enumerate(target.to_numpy()): + if np.isnan(t): + out[i] = cur if not np.isnan(cur) else 1.0 + continue + if np.isnan(cur) or abs(t - cur) > band * max(cur, 0.25): + cur = t + out[i] = cur + return pd.Series(out, index=target.index) + + +def lever(r: pd.Series, rf: pd.Series, exposure: pd.Series) -> pd.Series: + borrow = (exposure - 1).clip(lower=0) * (rf.reindex(r.index) + SPREAD) / 252 + trade = exposure.diff().abs().fillna(0.0) * COST + return exposure * r - borrow - trade + + +def managed(r: pd.Series, rf: pd.Series, target: float, cap: float, kind: str = "rv126", + band: float = 0.1) -> tuple[pd.Series, pd.Series]: + exp = banded((target / vol_hat(r, kind)).clip(upper=cap), band) + return lever(r, rf, exp), exp + + +def stats(r: pd.Series, rf: pd.Series) -> dict: + r = r.dropna() + e = (1 + r).cumprod() + yrs = len(r) / 252 + yr = (1 + r).groupby(r.index.year).prod() - 1 + roll3 = e.pct_change(756).dropna() + return {"cagr": e.iloc[-1] ** (1 / yrs) - 1, + "sharpe": (r - rf.reindex(r.index) / 252).mean() / r.std() * np.sqrt(252), + "vol": r.std() * np.sqrt(252), "mdd": (e / e.cummax() - 1).min(), + "worst_yr": yr.min(), "yrs_below_20": int((yr.iloc[1:-1] < 0.20).sum()), + "n_yrs": len(yr) - 2, "worst_3y": (1 + roll3.min()) ** (1 / 3) - 1} + + +def fmt(name: str, s: dict, exp: pd.Series | None = None) -> str: + x = f" avg exp {exp.mean():.2f}" if exp is not None else "" + return (f"{name:46s} cagr {s['cagr']*100:5.1f}% xsh {s['sharpe']:.2f} vol {s['vol']*100:4.1f}% " + f"mdd {s['mdd']*100:6.1f}% worst yr {s['worst_yr']*100:6.1f}% worst 3y/yr {s['worst_3y']*100:5.1f}% " + f"yrs<20% {s['yrs_below_20']}/{s['n_yrs']}{x}") + + +def main() -> None: + R = pd.read_pickle(os.path.join(RESULTS, "us_final_sleeves.pkl")) + extra = os.path.join(RESULTS, "short_horizon_sleeves.pkl") + if os.path.exists(extra): + R = R.join(pd.read_pickle(extra), how="left") + rf = tbill(R.index) + + print("=== single sleeves, static vs vol-managed (2006- for stock sleeves, 2010- trend)") + for c in ["momentum", "lowrisk", "trend", "spy"]: + r = R[c].dropna() + cap = 1.0 if c == "trend" else 2.0 + print(fmt(f"{c} as is", stats(r, rf))) + for target in (0.15, 0.25): + for kind in ("rv126", "ewma20"): + m, e = managed(r, rf, target, cap, kind) + print(fmt(f" {c} vt{int(target*100)} {kind} cap {cap:.0f}x", stats(m, rf), e)) + + span = R.loc[SPAN[0]:SPAN[1]].fillna(0.0) + three = {"momentum": 1 / 3, "lowrisk": 1 / 3, "trend": 1 / 3} + book = monthly_blend(span, three) + print("\n=== 3-sleeve book 2011-2026: static leverage vs managed") + for L in (1.0, 1.5, 2.0, 2.5): + e = pd.Series(L, index=book.index) + print(fmt(f"static {L:.1f}x", stats(lever(book, rf, e), rf), e)) + for target in (0.20, 0.25, 0.30, 0.35, 0.40): + for kind in ("rv126", "ewma20"): + m, e = managed(book, rf, target, 2.5, kind) + print(fmt(f"book vt{int(target*100)} {kind} cap 2.5x", stats(m, rf), e)) + + print("\n=== sleeve-level management, then equal capital (stock sleeves cap 2x, trend 1x)") + for target in (0.15, 0.20, 0.25, 0.30): + parts = {} + for c in ("momentum", "lowrisk", "trend"): + cap = 1.0 if c == "trend" else 2.0 + parts[c], _ = managed(R[c].dropna(), rf, target, cap) + P = pd.DataFrame(parts).loc[SPAN[0]:SPAN[1]].fillna(0.0) + print(fmt(f"sleeves vt{int(target*100)} rv126, 1/3 each", stats(monthly_blend(P, three), rf))) + + print("\n=== the same leverage through 2008 (momentum + lowrisk 50/50, 2006-; no trend sleeve before 2010)") + two = monthly_blend(R.loc[LONG[0]:LONG[1]].fillna(0.0), {"momentum": .5, "lowrisk": .5}) + for L in (1.0, 2.0): + e = pd.Series(L, index=two.index) + print(fmt(f"static {L:.1f}x", stats(lever(two, rf, e), rf), e)) + for target in (0.25, 0.35): + m, e = managed(two, rf, target, 2.5, "ewma20") + print(fmt(f"vt{int(target*100)} ewma20 cap 2.5x", stats(m, rf), e)) + + print("\n=== alpha left after hedging market beta (rolling 126d beta, futures hedge), 2011-") + spy = R["spy"] + hedged = {} + for c in ["momentum", "lowrisk", "trend"] + [c for c in R.columns if c.endswith("e21")]: + r = R[c] + cov = (r * spy).rolling(126, min_periods=63).mean() - r.rolling(126, min_periods=63).mean() * spy.rolling(126, min_periods=63).mean() + beta = (cov / spy.rolling(126, min_periods=63).var()).shift(1).clip(0, 2) + hedged[c] = r - beta * (spy - rf / 252) + print(fmt(f"hedged {c} (avg beta {beta.loc[SPAN[0]:].mean():.2f})", stats(hedged[c].loc[SPAN[0]:].dropna(), rf))) + + print("\n=== yearly, book vt30 ewma20 cap 2.5x vs 3-sleeve 1x vs SPY") + m, e = managed(book, rf, 0.30, 2.5, "ewma20") + yr = pd.DataFrame({"vt30": (1 + m).groupby(m.index.year).prod() - 1, + "1x": (1 + book).groupby(book.index.year).prod() - 1, + "spy": (1 + span["spy"]).groupby(span.index.year).prod() - 1, + "avg_exp": e.groupby(e.index.year).mean()}) + print((yr * [100, 100, 100, 1]).round(1).to_string()) + + +if __name__ == "__main__": + main() diff --git a/research/lab/panel_mid.py b/research/lab/panel_mid.py new file mode 100644 index 0000000..ceb1de7 --- /dev/null +++ b/research/lab/panel_mid.py @@ -0,0 +1,304 @@ +""" +Point-in-time US mid- and small-cap research Panels (see engine.Panel), built +from iShares IJH / IJR month-end holdings and Yahoo bars fetched by +research/lab/fetch_us_mid.py. + + from research.lab.panel_mid import load_us_mid + p = load_us_mid("mid") # S&P MidCap 400 members, 2006-09-29 → + p = load_us_mid("small") # S&P SmallCap 600 + p = load_us_mid("midsmall") # the union (S&P 400 + 600) + +``member`` is True when the ticker is in the fund's latest month-end snapshot +on or before that date AND Yahoo has a close that day, so index changes show +up at the next month-end (up to a month late; no snapshots exist for +2017-01..2017-06, so the Dec-2016 list stands until Jul-2017). + +Same-security check. The holdings files carry each stock's as-traded price; +Yahoo's close with later splits undone reproduces it to ~1e-8. A snapshot +row counts only when the two agree (``price_check``), which rejects symbols +Yahoo now serves for another company (Barnes Group's "B" is Barrick, +CoreSite's "COR" is Cencora). Most delisted names are missing on Yahoo: +``survivorship_coverage`` measures how much of each snapshot has prices, and +``member_vs_etf`` compares the covered members with the ETF. + +Cleaning: Yahoo's stray ticks (``bad_prints``) are removed, and a series that +turns into another security mid-membership (``_splices``) is cut there. +can_buy = can_sell = has a price. Costs are 20 bp a side (mid/small spreads +are wider than the S&P 500's). Benchmark: IJH (mid), IJR (small), their +daily equal blend (midsmall); extra["etf"] holds IJH, IJR, MDY, IWM, SPY, +QQQ, BIL, SHY closes, extra["index_weight"] each name's weight in its fund +at the latest snapshot. + +Caches: data/research/us_mid_bars.pkl (all bars) and +data/research/us_{mid,small,midsmall}.panel.pkl; delete them after refetching. +""" + +from __future__ import annotations + +import os +import pickle + +import numpy as np +import pandas as pd + +from research.lab.engine import Costs, Panel +from research.lab.fetch_us_mid import ETFS, load_bars, load_holdings + +RAW_DIR = os.path.join("data", "research") +FUNDS = {"mid": ("IJH",), "small": ("IJR",), "midsmall": ("IJH", "IJR")} +START = "2005-01-01" +PRICE_TOL = 0.02 # |log(iShares price / Yahoo as-traded close)|; true matches are ~1e-8 +JUMP = np.log(3) # a one-day move this large is checked (bad print, splice) +STABLE_TOL = 0.005 # ... or an offset (merger ratio, spin-off) this steady over 3 snapshots +MAX_OFFSET = 0.5 # largest steady offset accepted +SECTOR_ALIASES = {"Communication": "Communication Services", + "Telecommunications": "Communication Services", + "Telecommunication Services": "Communication Services"} + + +def _cached(name: str, build): + path = os.path.join(RAW_DIR, f"{name}.panel.pkl") + if os.path.exists(path): + with open(path, "rb") as fh: + return pickle.load(fh) + p = build() + with open(path, "wb") as fh: + pickle.dump(p, fh, protocol=pickle.HIGHEST_PROTOCOL) + return p + + +def bad_prints(close: np.ndarray, volume: np.ndarray, window: int = 5) -> np.ndarray: + """Bars more than JUMP (3x) away from the last good close that either did + not trade (zero volume, or under 1% of the trailing 60-day median) or + come back within 1.5x of it inside ``window`` days: Yahoo's stray ticks.""" + n = len(close) + bad = np.zeros(n, bool) + medv = pd.Series(volume).rolling(60, min_periods=10).median().shift(1).to_numpy() + fin = np.isfinite(close) & (close > 0) + last = np.nan + for i in range(n): + if not fin[i]: + continue + if last == last and abs(np.log(close[i] / last)) > JUMP: + thin = volume[i] == 0 or (medv[i] == medv[i] and volume[i] < 0.01 * medv[i]) + nxt = close[i + 1:i + 1 + window] + nxt = nxt[np.isfinite(nxt) & (nxt > 0)] + if thin or (len(nxt) and np.any(np.abs(np.log(nxt / last)) < np.log(1.5))): + bad[i] = True + continue + last = close[i] + return bad + + +def _bars_all() -> dict: + """Wide adjusted OHLCV + as-traded close for every ever-member and ETF, + on SPY's trading calendar from 2005 (cached: rebuilding reads ~3k CSVs). + ``bad_prints`` are removed from OHLCV (``bad`` marks them); ``raw_close`` + is left as Yahoo serves it for the price check.""" + path = os.path.join(RAW_DIR, "us_mid_bars.pkl") + if os.path.exists(path): + with open(path, "rb") as fh: + return pickle.load(fh) + hold = load_holdings() + tickers = sorted(set(hold["ticker"]) | set(ETFS)) + spy = load_bars("SPY") + dates = spy.index[spy.index >= START] + fields = ("open", "high", "low", "close", "volume", "raw_close") + arr = {k: np.full((len(dates), len(tickers)), np.nan) for k in fields} + have = [] + for j, t in enumerate(tickers): + b = load_bars(t) + if b is None: + continue + b = b.reindex(dates) + for k in fields: + arr[k][:, j] = b[k].to_numpy(float) + have.append(t) + out = {k: pd.DataFrame(v, index=dates, columns=tickers)[have] for k, v in arr.items()} + bad = np.column_stack([bad_prints(out["close"][t].to_numpy(), out["volume"][t].to_numpy()) + for t in have]) + out["bad"] = pd.DataFrame(bad, index=dates, columns=have) + for k in ("open", "high", "low", "close", "volume"): + out[k] = out[k].mask(out["bad"]) + with open(path, "wb") as fh: + pickle.dump(out, fh, protocol=pickle.HIGHEST_PROTOCOL) + return out + + +def price_check(hold: pd.DataFrame, raw_close: pd.DataFrame) -> pd.DataFrame: + """Per holding row: log of the iShares share price over Yahoo's as-traded + close on the snapshot date (last close within 5 days). A row matches when + they agree within PRICE_TOL, or when the ratio holds within STABLE_TOL + over three consecutive snapshots of the ticker: Yahoo scales some + histories by a constant (a spin-off: VSH before 2010 is off by 10%; a + merger's exchange ratio: DCOM before 2021 is old Dime x 0.648). The + result is then the majority of each snapshot and its two neighbours, + which drops a lone coincidental match of another company and keeps a + lone event-month glitch. ``has_px``: Yahoo has a close.""" + rc = raw_close.ffill(limit=5) + h = hold[["fund", "date", "ticker", "price", "weight"]].sort_values(["fund", "ticker", "date"]) + rows = rc.index.get_indexer(h["date"], method="ffill") + cols = rc.columns.get_indexer(h["ticker"]) + ok = (rows >= 0) & (cols >= 0) + yc = np.full(len(h), np.nan) + yc[ok] = rc.to_numpy()[rows[ok], cols[ok]] + h["yahoo_close"] = yc + h["has_px"] = np.isfinite(yc) & (yc > 0) + with np.errstate(divide="ignore", invalid="ignore"): + lr = np.log(h["price"] / h["yahoo_close"]) + h["log_ratio"] = lr + key = [h["fund"], h["ticker"]] + g = lr.groupby(key) + near_prev = (lr - g.shift(1)).abs() <= STABLE_TOL + near_next = (lr - g.shift(-1)).abs() <= STABLE_TOL + # an offset held over three snapshots in a row (two unrelated prices rarely + # do that), and the run's first and last snapshot + core = near_prev & near_next & (lr.abs() <= MAX_OFFSET) + cg = core.groupby(key) + steady = (core | (near_next & cg.shift(-1, fill_value=False).astype(bool)) + | (near_prev & cg.shift(1, fill_value=False).astype(bool))) + raw = h["has_px"] & ((lr.abs() <= PRICE_TOL) | steady | h["price"].isna()) + hp = h["has_px"].groupby(key) + has_prev = hp.shift(1, fill_value=False).astype(bool) + has_next = hp.shift(-1, fill_value=False).astype(bool) + rg = raw.groupby(key) + votes = (raw.astype(int) + (rg.shift(1, fill_value=False).astype(bool) & has_prev).astype(int) + + (rg.shift(-1, fill_value=False).astype(bool) & has_next).astype(int)) + both = has_prev & has_next + h["match_raw"] = raw + h["match"] = h["has_px"] & np.where(both, votes >= 2, raw) + return h.sort_index() + + +def _splices(chk: pd.DataFrame, bars: dict) -> pd.DataFrame: + """Dates whose Yahoo prices belong to another security spliced onto the + series (Intermec's "IN" turns into a penny stock in Nov 2012): the + ticker matches at one snapshot, has a price that does not match at the + fund's next snapshot, and the series jumps (or prints badly) in between. + Prices are dropped from that day until the ticker matches again.""" + close = bars["close"] + jump = (np.log(close.ffill()).diff().abs() > JUMP) | bars["bad"] + drop = pd.DataFrame(False, index=close.index, columns=close.columns) + matched = chk[chk["match"]].groupby("ticker")["date"].apply(lambda d: np.sort(d.unique())) + for fund, gf in chk.groupby("fund"): + snaps = np.sort(gf["date"].unique()) + nxt = dict(zip(snaps[:-1], snaps[1:])) + for t, g in gf.sort_values("date").groupby("ticker"): + if t not in close.columns: + continue + d, m, px = g["date"].to_numpy(), g["match"].to_numpy(), g["has_px"].to_numpy() + for i in range(len(g) - 1): + if not (m[i] and px[i + 1] and not m[i + 1] and nxt.get(d[i]) == d[i + 1]): + continue + win = jump.loc[(jump.index > d[i]) & (jump.index <= d[i + 1]), t] + if not win.any(): + continue + start = win.idxmax() + later = matched.get(t, np.array([], dtype="datetime64[ns]")) + later = later[later > d[i + 1]] + end = later[0] if len(later) else close.index[-1] + pd.Timedelta(days=1) + drop.loc[(drop.index >= start) & (drop.index < end), t] = True + return drop + + +def _membership(chk: pd.DataFrame, dates: pd.DatetimeIndex, tickers: list[str]) -> pd.DataFrame: + """In the fund's latest snapshot on or before each date, with a matching price.""" + col = {t: i for i, t in enumerate(tickers)} + arr = np.zeros((len(dates), len(tickers)), bool) + for fund, g in chk.groupby("fund"): + snaps = sorted(g["date"].unique()) + for i, s in enumerate(snaps): + e = snaps[i + 1] if i + 1 < len(snaps) else dates[-1] + pd.Timedelta(days=1) + rows = (dates >= s) & (dates < e) + names = g.loc[(g["date"] == s) & g["match"], "ticker"] + cols = [col[t] for t in names if t in col] + arr[np.ix_(rows, cols)] = True + return pd.DataFrame(arr, index=dates, columns=tickers) + + +def _index_weight(chk: pd.DataFrame, dates, tickers) -> pd.DataFrame: + """Weight in its fund at the latest snapshot (held constant until the next).""" + w = chk.pivot_table(index="date", columns="ticker", values="weight", aggfunc="sum") + w = w.reindex(columns=tickers) + return w.reindex(dates.union(w.index)).ffill().reindex(dates).fillna(0.0).astype("float32") + + +def _build(which: str) -> Panel: + funds = FUNDS[which] + bars = _bars_all() + hold = load_holdings() + hold = hold[hold["fund"].isin(funds)] + chk = price_check(hold, bars["raw_close"]) + close_all = bars["close"] + dates = close_all.index + tickers = sorted(set(chk.loc[chk["match"], "ticker"]) & set(close_all.columns) - set(ETFS)) + spliced = _splices(chk, bars)[tickers] + f = {k: bars[k][tickers].mask(spliced) for k in ("open", "high", "low", "close", "volume")} + for k in ("open", "high", "low"): + f[k] = f[k].where(f[k] > 0) + f["close"] = f["close"].where(f["close"] > 0) + has_px = f["close"].notna() + member = _membership(chk, dates, tickers) & has_px + etf = close_all[[e for e in ETFS if e in close_all.columns]] + bil = etf["BIL"].pct_change(fill_method=None) + shy = etf["SHY"].pct_change(fill_method=None) + cash_ret = bil.where(bil.notna(), shy).fillna(0.0) + if which == "midsmall": # S&P 400 + 600 → equal blend of the two funds + r = etf[["IJH", "IJR"]].pct_change(fill_method=None).mean(axis=1).fillna(0.0) + bench = (1 + r).cumprod() + else: + bench = etf[funds[0]] + # the ticker's latest GICS sector (as panel.py, today's sector for all dates) + sector = (hold[~hold["sector"].isin(["-", "", "Other"])].sort_values("date") + .groupby("ticker")["sector"].last().replace(SECTOR_ALIASES).reindex(tickers)) + return Panel( + market=f"us_{which}", open=f["open"], high=f["high"], low=f["low"], close=f["close"], + volume=f["volume"], member=member, can_buy=has_px, can_sell=has_px, + cash_ret=cash_ret, benchmark=bench, + costs=Costs(buy=0.002, sell=0.002, fixed_fee=0.0, capital=10_000.0), + extra={"etf": etf, "sector": sector, + "index_weight": _index_weight(chk, dates, tickers)}, + ) + + +def load_us_mid(which: str = "mid") -> Panel: + if which not in FUNDS: + raise ValueError(f"which must be one of {list(FUNDS)}") + return _cached(f"us_{which}", lambda: _build(which)) + + +# --------------------------------------------------------------------------- +# diagnostics +# --------------------------------------------------------------------------- + +def survivorship_coverage(which: str = "midsmall") -> pd.DataFrame: + """Per fund and snapshot: members, and the share (by count and by index + weight) with a Yahoo price for the same security on that date.""" + bars = _bars_all() + hold = load_holdings() + chk = price_check(hold[hold["fund"].isin(FUNDS[which])], bars["raw_close"]) + chk["w_ok"] = chk["weight"] * chk["match"] + chk["w_px"] = chk["weight"] * chk["has_px"] + g = chk.groupby(["fund", "date"]) + return pd.DataFrame({ + "n": g.size(), "n_ok": g["match"].sum(), "count_cov": g["match"].mean(), + "weight_cov": g["w_ok"].sum() / g["weight"].sum(), + "count_any_px": g["has_px"].mean(), + "weight_any_px": g["w_px"].sum() / g["weight"].sum(), + }) + + +def member_vs_etf(p: Panel) -> pd.DataFrame: + """Calendar-year return of the equal-weight member portfolio (daily + rebalanced), the index-weighted covered members, and the ETF.""" + r = p.close.pct_change(fill_method=None) + prev = p.member.shift(1, fill_value=False) + ew = r.where(prev).mean(axis=1) + w = p.extra["index_weight"].shift(1).where(prev).astype(float) + cw = (r.where(prev) * w).sum(axis=1) / w.where(r.where(prev).notna()).sum(axis=1) + bench = p.benchmark.pct_change(fill_method=None) + df = pd.DataFrame({"ew_members": ew, "cw_members": cw, "etf": bench}).loc[prev.any(axis=1)] + yr = (1 + df.fillna(0.0)).groupby(df.index.year).prod() - 1 + yr["ew_minus_etf"] = yr["ew_members"] - yr["etf"] + yr["cw_minus_etf"] = yr["cw_members"] - yr["etf"] + return yr diff --git a/research/lab/qqq_scan.py b/research/lab/qqq_scan.py new file mode 100644 index 0000000..62c9290 --- /dev/null +++ b/research/lab/qqq_scan.py @@ -0,0 +1,143 @@ +""" +Which strategies and factors carry alpha against QQQ, before and after 2017? + + PYTHONPATH=. uv run python -m research.lab.fetch_long_history # once (QQQ) + PYTHONPATH=. uv run --with numba python -m research.lab.qqq_scan [--part sleeves|factors|mega] + +sleeves the current sleeves (us_final, short_horizon) and TrendRider +factors every price factor (ranks_us) and every EDGAR fundamental, long-only + top 20 among PIT S&P 500 members, every 21 days, buffer 2 +mega the same idea inside a QQQ-like universe: the 100 members with the + highest 126-day dollar volume (market caps are missing for many + delisted names, dollar volume is not), top 10 / 20 + +Alpha is the intercept of daily excess returns on QQQ's excess return, with +a Newey-West t-stat (research/lab/alpha.py). Periods: design 2006-09..2016, +holdout 2017..2026-09. A factor whose alpha is positive in both, with t > 2 +over the whole span, is a candidate; everything else is noise or beta. +Writes data/research/results/qqq_scan_{part}.pkl (daily returns) and .csv. +""" + +from __future__ import annotations + +import argparse +import dataclasses +import os +import pickle +import warnings + +import pandas as pd + +from research.lab import fundamentals as fu +from research.lab.alpha import alpha_stats, bench_returns, fmt_alpha +from research.lab.factor_scan import RESULTS, load +from research.lab.factors import _dollar_volume +from research.lab.rotation import factor_ranks +from research.lab.us_sweep import _blend, run, sector_rank + +warnings.filterwarnings("ignore") +LONG = ("2006-09-27", "2026-09-24") +DESIGN = ("2006-09-27", "2016-12-31") +HOLDOUT = ("2017-01-01", "2026-09-24") + + +def report(name: str, r: pd.Series, rows: list, bench: pd.Series, spy: pd.Series) -> None: + r = r.dropna() + full = alpha_stats(r, bench.reindex(r.index)) + d = alpha_stats(r.loc[DESIGN[0]:DESIGN[1]], bench) + h = alpha_stats(r.loc[HOLDOUT[0]:HOLDOUT[1]], bench) + s = alpha_stats(r, spy.reindex(r.index)) + rows.append({"name": name, **{f"full_{k}": v for k, v in full.items()}, + **{f"design_{k}": v for k, v in d.items()}, **{f"holdout_{k}": v for k, v in h.items()}, + "spy_alpha": s["alpha"], "spy_t": s["t"]}) + print(f"{fmt_alpha(name, full)} | design a {d['alpha']*100:5.1f}% t {d['t']:5.2f} | " + f"holdout a {h['alpha']*100:5.1f}% t {h['t']:5.2f} cagr {h['cagr']*100:5.1f}% vs {h['bench_cagr']*100:5.1f}% " + f"| vs SPY a {s['alpha']*100:5.1f}% t {s['t']:5.2f}", flush=True) + + +def part_sleeves(bench, spy, rows) -> dict: + R = pd.read_pickle(os.path.join(RESULTS, "us_final_sleeves.pkl")) + extra = os.path.join(RESULTS, "short_horizon_sleeves.pkl") + if os.path.exists(extra): + R = R.join(pd.read_pickle(extra)) + with open(os.path.join("data", "research", "letf_history.pkl"), "rb") as fh: + tr = pickle.load(fh)["equity"] + R["trend_rider_v7_vt24 (synthetic)"] = tr["trend_rider_v7_vt24 [lag]"].pct_change().reindex(R.index) + for c in R.columns: + report(c, R[c], rows, bench, spy) + return {} + + +def part_factors(p, bench, spy, rows) -> dict: + pr = factor_ranks(p, "us") + F = fu.load(p) + scores = {k: v for k, v in pr.items()} + for k, v in F.items(): + if k == "mcap": + continue + scores[f"f:{k}"] = v.where(p.member).rank(axis=1, pct=True) + scores[f"f:{k} (sector)"] = sector_rank(v, p) + out = {} + for name, sc in scores.items(): + st = run(p, sc.where(p.member), [LONG], top_n=20, buffer=2.0, every=21) + out[name] = st["_rets"] + report(name, st["_rets"], rows, bench, spy) + return out + + +def mega_panel(p, n: int = 100): + dv = _dollar_volume(p).rolling(126, min_periods=63).mean().where(p.member) + top = dv.rank(axis=1, ascending=False) <= n + return dataclasses.replace(p, member=p.member & top) + + +def part_mega(p, bench, spy, rows) -> dict: + pm = mega_panel(p) + pr = factor_ranks(p, "us") + F = fu.load(p) + ew = pm.close.notna().astype(float).where(pm.member) + scores = { + "mega EW (100)": (ew, 10000), + "mega overnight+rec": (_blend(pr["overnight_mom"], pr["rec_126"]), None), + "mega mom_12_1": (pr["mom_12_1"], None), "mega resid_mom": (pr["resid_mom"], None), + "mega high_52w": (pr["high_52w"], None), "mega sharpe_mom": (pr["sharpe_mom"], None), + "mega overnight_mom": (pr["overnight_mom"], None), "mega rec_126": (pr["rec_126"], None), + "mega mom_6_1": (pr["mom_6_1"], None), "mega trend_hzz": (pr["trend_hzz"], None), + "mega low_beta": (pr["low_beta"], None), "mega ivol": (pr["ivol"], None), + "mega gp_a": (F["gp_a"].rank(axis=1, pct=True), None), "mega roe": (F["roe"].rank(axis=1, pct=True), None), + "mega mom_12_1+gp_a": (_blend(pr["mom_12_1"], F["gp_a"].rank(axis=1, pct=True)), None), + } + out = {} + for name, (sc, n_all) in scores.items(): + for n in ([n_all] if n_all else (10, 20)): + st = run(pm, sc.where(pm.member), [LONG], top_n=n, buffer=2.0, every=21) + key = name if n_all else f"{name} top{n}" + out[key] = st["_rets"] + report(key, st["_rets"], rows, bench, spy) + return out + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--part", nargs="*", default=["sleeves", "mega", "factors"]) + args = ap.parse_args() + p, _ = load("us") + bench = bench_returns(p.dates) + spy = bench_returns(p.dates, "SPY") + report("QQQ itself", bench.loc[LONG[0]:], [], bench, spy) + for part in args.part: + print(f"\n=== {part}: alpha vs QQQ (full 2006-2026) | design 2006-16 | holdout 2017-26 | vs SPY", flush=True) + rows = [] + if part == "sleeves": + out = part_sleeves(bench, spy, rows) + elif part == "factors": + out = part_factors(p, bench, spy, rows) + else: + out = part_mega(p, bench, spy, rows) + pd.DataFrame(rows).to_csv(os.path.join(RESULTS, f"qqq_scan_{part}.csv"), index=False) + if out: + pd.DataFrame(out).to_pickle(os.path.join(RESULTS, f"qqq_scan_{part}.pkl")) + + +if __name__ == "__main__": + main() diff --git a/research/lab/short_horizon.py b/research/lab/short_horizon.py new file mode 100644 index 0000000..d1715cb --- /dev/null +++ b/research/lab/short_horizon.py @@ -0,0 +1,106 @@ +""" +Short-horizon and residual stock sleeves for the US book: do they add a +return source the three current sleeves (momentum, low risk, trend) lack? + + PYTHONPATH=. uv run --with numba python -m research.lab.short_horizon + +Candidates, price/volume only (S&P 500 PIT, next-open fills, 10 bp): + resid_mom residual 12-1 momentum / residual vol (Blitz, Huij & Martens 2011; + Blitz, Hanauer & Vidojevic 2020) + stmom_hiturn 1-month momentum among high-turnover names (Medhat & Schmeling 2022) + rev_sector_21 industry-adjusted 1-month reversal (Da, Liu & Schaumburg 2014) + rev_resid_21 residual 1-month reversal + seasonality same-calendar-month return, years 1-10 (Heston & Sadka 2008) + short_combo price-only subset of the Robeco short-horizon composite (Blitz et + al. 2023): short-term momentum for high-turnover names, residual + reversal for the rest, seasonality, low 1-month idiosyncratic vol + +Each is run as top 20 with a rank buffer (hold while rank <= 40), at two +rebalance speeds. Reports stand-alone stats, monthly correlation with the +current sleeves (data/research/results/us_final_sleeves.pkl, from us_final), +and the 2011-2026 book with the candidate as a fourth equal sleeve. + +Last, concentration: the momentum sleeve (overnight + rec_126) with 3 or 5 +names instead of 10, the other thing a small account can do that a fund +cannot. +""" + +from __future__ import annotations + +import os +import warnings + +import pandas as pd + +from research.lab.factor_scan import RESULTS, load +from research.lab.rotation import factor_ranks +from research.lab.us_final import describe, line, monthly_blend +from research.lab.us_sweep import _blend, fmt, log_trials, run + +warnings.filterwarnings("ignore") + +SPAN = ("2011-01-03", "2026-09-24") +LONG = ("2006-09-27", "2026-09-24") + + +def candidates(p) -> dict[str, pd.DataFrame]: + pr = factor_ranks(p, "us") + hi = pr["stmom_hiturn"].notna() + directional = pr["stmom_hiturn"].where(hi, pr["rev_resid_21"]) + return { + "resid_mom": pr["resid_mom"], + "stmom_hiturn": pr["stmom_hiturn"], + "rev_sector_21": pr["rev_sector_21"], + "rev_resid_21": pr["rev_resid_21"], + "seasonality": pr["seasonality"], + "short_combo": _blend(directional, pr["seasonality"], pr["ivol_20"]), + "resid_mom+short_combo": _blend(pr["resid_mom"], directional, pr["seasonality"], pr["ivol_20"]), + } + + +def main() -> None: + p, _ = load("us") + base = pd.read_pickle(os.path.join(RESULTS, "us_final_sleeves.pkl")) + cur = base[["momentum", "lowrisk", "trend"]] + rows, R = [], {} + for name, score in candidates(p).items(): + for every in (5, 21): + st = run(p, score.where(p.member), [LONG], top_n=20, buffer=2.0, every=every) + key = f"{name} e{every}" + R[key] = st["_rets"] + d = describe(st["_rets"].loc[SPAN[0]:]) + dl = describe(st["_rets"]) + rows.append({"score": key, **{k: v for k, v in st.items() if not k.startswith("_")}}) + print(f"{line(key + ' 2011-', d)} turn {st['turn']:4.1f}x | 2006- sh {dl['sharpe']:.2f} " + f"mdd {dl['mdd']*100:5.1f}%", flush=True) + log_trials(rows, "short_horizon") + R = pd.DataFrame(R) + R.to_pickle(os.path.join(RESULTS, "short_horizon_sleeves.pkl")) + + span = pd.concat([cur, R], axis=1).loc[SPAN[0]:SPAN[1]].fillna(0.0) + monthly = (1 + span).resample("ME").prod() - 1 + print("\nmonthly correlation with current sleeves, 2011-2026:") + print(monthly.corr().loc[R.columns, ["momentum", "lowrisk", "trend"]].round(2).to_string()) + + three = {"momentum": 1 / 3, "lowrisk": 1 / 3, "trend": 1 / 3} + print("\n2011-2026 book, monthly rebalanced between sleeves:") + print(line("current 3 sleeves", describe(monthly_blend(span, three)))) + for c in R.columns: + four = {"momentum": .25, "lowrisk": .25, "trend": .25, c: .25} + swap = {"momentum": 1 / 3, c: 1 / 3, "trend": 1 / 3} + print(line(f"+ {c} (4 x 1/4)", describe(monthly_blend(span, four)))) + print(line(f" {c} replaces lowrisk", describe(monthly_blend(span, swap)))) + + print("\nconcentration: overnight + rec_126, every 10 days, 2006-2026") + pr = factor_ranks(p, "us") + mom = _blend(pr["overnight_mom"], pr["rec_126"]).where(p.member) + blocks = [(LONG[0], "2010-12-31"), ("2011-01-01", "2015-12-31"), ("2016-01-01", "2020-12-31"), + ("2021-01-01", LONG[1])] + for n in (3, 5, 10): + for buf in (1.0, 2.0): + st = run(p, mom, blocks, top_n=n, buffer=buf, every=10) + print(fmt(f"top {n} buffer {buf}", st), flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/us_target50.py b/research/lab/us_target50.py new file mode 100644 index 0000000..cf86e6f --- /dev/null +++ b/research/lab/us_target50.py @@ -0,0 +1,99 @@ +""" +What would it take to earn 50% a year, and at what drawdown? The US book +over 2006-2026, including 2008, with the trend sleeve rebuilt from synthetic +3x ETFs (research/lab/letf_history.py), then levered. + + PYTHONPATH=. uv run python -m research.lab.us_final # sleeves (once) + PYTHONPATH=. uv run python -m research.lab.letf_history # TrendRider on synthetic ETFs + PYTHONPATH=. uv run python -m research.lab.vol_trend # simple vol-targeted trend + PYTHONPATH=. uv run python -m research.lab.us_target50 + +Books (sleeves at equal capital, reset monthly): + current momentum + lowrisk + TrendRider V7 vt24 (the live us_three_sleeve) + simple momentum + lowrisk + QQQ trend (ma200, vol target 20%, cap 3x; a + mid-grid configuration, not the best one) +Then the leverage each book needs for 50% (static, and vol-managed as in +levered_book.py) and what that leverage did in 2008 and in the worst years. +""" + +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 +from research.lab.levered_book import fmt, lever, managed, stats, tbill +from research.lab.us_final import monthly_blend + +warnings.filterwarnings("ignore") +LONG = ("2006-09-27", "2026-09-24") +BLOCKS = [("2006-09-27", "2010-12-31"), ("2011-01-01", "2015-12-31"), ("2016-01-01", "2020-12-31"), + ("2021-01-01", "2026-09-24")] + + +def sleeves() -> pd.DataFrame: + R = pd.read_pickle(os.path.join(RESULTS, "us_final_sleeves.pkl"))[["momentum", "lowrisk", "spy"]] + with open(os.path.join("data", "research", "letf_history.pkl"), "rb") as fh: + tr = pickle.load(fh)["equity"]["trend_rider_v7_vt24 [lag]"] + with open(os.path.join(RESULTS, "vol_trend.pkl"), "rb") as fh: + vt = pickle.load(fh)["returns"]["QQQ ma200 vt20 cap3"] + R["trend_rider"] = tr.pct_change().reindex(R.index) + R["qqq_trend"] = vt.reindex(R.index) + return R.loc[LONG[0]:LONG[1]] + + +def blocks(r: pd.Series, rf: pd.Series) -> str: + return " ".join(f"{stats(r.loc[a:b], rf)['cagr']*100:5.1f}%/{stats(r.loc[a:b], rf)['sharpe']:4.2f}" + for a, b in BLOCKS) + + +def main() -> None: + R = sleeves() + rf = tbill(R.index) + out = {} + print("=== sleeves 2006-2026 (trend sleeves on synthetic 3x ETFs, lag fills)") + print(f"{'':46s}{'':104s} blocks 06-10 / 11-15 / 16-20 / 21-26") + for c in R.columns: + s = stats(R[c].dropna(), rf) + print(fmt(c, s), "|", blocks(R[c].fillna(0.0), rf)) + books = {"current": {"momentum": 1 / 3, "lowrisk": 1 / 3, "trend_rider": 1 / 3}, + "simple": {"momentum": 1 / 3, "lowrisk": 1 / 3, "qqq_trend": 1 / 3}} + m = (1 + R.fillna(0.0)).resample("ME").prod() - 1 + print("\nmonthly correlation:\n", m.corr().round(2).to_string()) + print("\n=== books 2006-2026 and the leverage for 50%") + for name, w in books.items(): + b = monthly_blend(R.fillna(0.0), w) + s = stats(b, rf) + out[name] = {"1x": s} + print(fmt(f"{name} 1x", s), "|", blocks(b, rf)) + for L in (1.5, 2.0, 3.0): + e = pd.Series(L, index=b.index) + x = lever(b, rf, e) + out[name][f"{L}x"] = stats(x, rf) + print(fmt(f" {name} static {L}x", out[name][f"{L}x"], e), "|", blocks(x, rf)) + for target in (0.25, 0.35, 0.45): + x, e = managed(b, rf, target, 3.0, "ewma20") + out[name][f"vt{int(target*100)}"] = stats(x, rf) + print(fmt(f" {name} vt{int(target*100)} ewma20 cap 3x", out[name][f"vt{int(target*100)}"], e), + "|", blocks(x, rf)) + print("\n=== yearly returns: current book 1x / 2x static / vt35, SPY") + b = monthly_blend(R.fillna(0.0), books["current"]) + x2 = lever(b, rf, pd.Series(2.0, index=b.index)) + x35, _ = managed(b, rf, 0.35, 3.0, "ewma20") + yr = pd.DataFrame({k: (1 + v).groupby(v.index.year).prod() - 1 + for k, v in {"1x": b, "2x": x2, "vt35": x35, "spy": R["spy"].fillna(0.0)}.items()}) + print((yr * 100).round(1).to_string()) + for k in out: + for kk in out[k]: + out[k][kk] = {a: float(v) for a, v in out[k][kk].items()} + with open(os.path.join(RESULTS, "us_target50.json"), "w") as fh: + json.dump(out, fh, indent=1) + + +if __name__ == "__main__": + main() diff --git a/research/lab/vol_trend.py b/research/lab/vol_trend.py new file mode 100644 index 0000000..d606bb3 --- /dev/null +++ b/research/lab/vol_trend.py @@ -0,0 +1,125 @@ +""" +Leveraged trend without TrendRider's machinery: a trend filter times a +volatility target, implemented with 3x ETFs, tested on synthetic ETFs from +1990 and with parameters picked walk-forward. + + PYTHONPATH=. uv run python -m research.lab.fetch_long_history # once + PYTHONPATH=. uv run python -m research.lab.vol_trend + +Exposure to the index (SPY or QQQ, total return) decided on closes through +t-1: + E_t = 1{index > its n-day average} x min(cap, target / sigma_t) +sigma_t is the index's EWMA vol (20-day half-life). E/3 is held in the 3x +ETF (UPRO / TQQQ, synthetic: research/lab/letf_history.py), the rest in the +2-year Treasury series. Holdings drift; the book trades back to target when +the signal flips or E is more than 25% off; 10 bp per unit traded. Fills: +"lag" trades at the close after the decision (pessimistic), "close" at the +decision close (optimistic). + +Moreira & Muir (2017) and Harvey et al. (2018) motivate the vol target; +Gayed & Bilello (2016) the moving-average switch. The grid is small and +round; the walk-forward chain picks, each January, the configuration with +the best Sharpe over the previous ten years. +""" + +from __future__ import annotations + +import itertools +import os +import pickle +import warnings + +import numpy as np +import pandas as pd + +from research.lab.letf_history import BLOCKS, SRC, build + +warnings.filterwarnings("ignore") +OUT = os.path.join("data", "research", "results", "vol_trend.pkl") +COST, BAND = 0.001, 0.25 +GRID = list(itertools.product(("SPY", "QQQ"), (100, 150, 200), (0.15, 0.20, 0.25, 0.30, 0.40), (1.5, 3.0))) +LETF = {"SPY": "UPRO", "QQQ": "TQQQ"} + + +def backtest(syn: pd.DataFrame, index: str, n: int, target: float, cap: float, + lag: bool = True) -> tuple[pd.Series, pd.Series]: + px = syn[index] + r = px.pct_change().fillna(0.0) + on = (px > px.rolling(n).mean()).astype(float) + sigma = np.sqrt((r ** 2).ewm(halflife=20, min_periods=40).mean() * 252) + E = (on * (target / sigma).clip(upper=cap)).fillna(0.0) + E.iloc[:n] = 0.0 + tgt = (E / 3.0).shift(2 if lag else 1).fillna(0.0).to_numpy() # row t = weight held over t + rl = syn[LETF[index]].pct_change().fillna(0.0).to_numpy() + rs = syn["SHY"].pct_change().fillna(0.0).to_numpy() + w, out, held = 0.0, np.empty(len(tgt)), np.empty(len(tgt)) + for t in range(len(tgt)): + if t > 0: + g = w * (1 + rl[t]) + (1 - w) * (1 + rs[t]) + out[t] = g - 1 + w = w * (1 + rl[t]) / g + else: + out[t] = 0.0 + want = tgt[t + 1] if t + 1 < len(tgt) else w + if (want == 0) != (w == 0) or abs(want - w) > BAND * max(w, 1e-9): + out[t] -= abs(want - w) * COST * 2 # sell one leg, buy the other + w = want + held[t] = w * 3 + return pd.Series(out, index=syn.index), pd.Series(held, index=syn.index) + + +def sharpe(r: pd.Series) -> float: + return float(r.mean() / r.std() * np.sqrt(252)) if r.std() > 0 else 0.0 + + +def block_line(r: pd.Series) -> str: + cells = [] + for a, b in BLOCKS: + x = r.loc[a:b] + e = (1 + x).cumprod() + cells.append(f"{(e.iloc[-1] ** (252 / len(x)) - 1)*100:6.1f}%/{sharpe(x):4.2f}/{((e / e.cummax()) - 1).min()*100:4.0f}%") + return " ".join(c.rjust(17) for c in cells) + + +def walk_forward(R: pd.DataFrame, lookback_years: int = 10, start: int = 2000) -> tuple[pd.Series, dict]: + parts, picks = [], {} + for y in range(start, R.index[-1].year + 1): + hist = R.loc[f"{y - lookback_years}-01-01":f"{y - 1}-12-31"] + best = hist.apply(sharpe).idxmax() + picks[y] = best + parts.append(R.loc[f"{y}-01-01":f"{y}-12-31", best]) + return pd.concat(parts), picks + + +def main() -> None: + with open(SRC, "rb") as fh: + syn, _ = build(pickle.load(fh)) + syn = syn.loc["1989-01-01":] + R, H, params = {}, {}, {} + for idx, n, tv, cap in GRID: + key = f"{idx} ma{n} vt{int(tv*100)} cap{cap:g}" + params[key] = (idx, n, tv, cap) + R[key], H[key] = backtest(syn, idx, n, tv, cap) + R = pd.DataFrame(R).loc["1990-01-01":] + head = " ".join(f"{a[:4]}-{b[2:4]:>2s}".rjust(17) for a, b in BLOCKS) + full = R.apply(sharpe).sort_values(ascending=False) + print(f"CAGR/Sharpe/MaxDD by block, lag fills, synthetic 3x ETFs\n{'':28s}{head}") + for k in list(full.index[:8]) + list(full.index[len(full) // 2:len(full) // 2 + 2]) + list(full.index[-2:]): + print(f"{k:28s}{block_line(R[k])} avg E {H[k].loc['1990':].mean():.2f}") + print(f"\nfull-period Sharpe over the {len(full)}-config grid: median {full.median():.2f}, " + f"range {full.min():.2f}-{full.max():.2f}") + wf, picks = walk_forward(R) + print(f"\nwalk-forward (best trailing-10y Sharpe each January), 2000-2026:\n{'':28s}{head}") + print(f"{'walk-forward chain':28s}{block_line(pd.concat([R.iloc[:, 0].loc[:'1999'] * np.nan, wf]).fillna(0.0))}") + print("picks:", {y: p for y, p in picks.items() if y in (2000, 2005, 2010, 2015, 2020, 2025, 2026)}) + close_mode = {k: backtest(syn, *params[k], lag=False)[0].loc["1990-01-01":] for k in full.index[:3]} + print("\nsame top configs with close fills (optimistic):") + for k, r in close_mode.items(): + print(f"{k:28s}{block_line(r)}") + os.makedirs(os.path.dirname(OUT), exist_ok=True) + with open(OUT, "wb") as fh: + pickle.dump({"returns": R, "exposure": pd.DataFrame(H), "walk_forward": wf, "picks": picks}, fh) + + +if __name__ == "__main__": + main() diff --git a/strategies/qqq_trend.py b/strategies/qqq_trend.py new file mode 100644 index 0000000..3bb44ef --- /dev/null +++ b/strategies/qqq_trend.py @@ -0,0 +1,56 @@ +""" +QQQ trend with a volatility target, held through TQQQ and short Treasuries. + +Exposure to the Nasdaq-100 on each day, from QQQ closes through the day before: + + E = 1{QQQ > its 200-day average} x min(cap, target / EWMA vol of QQQ) + +E / 3 of the book is held in TQQQ (3x daily QQQ) and the rest in SHY, so +E = 1 is QQQ-like exposure with two thirds of the capital in Treasuries. +The 200-day switch is the classic trend rule (Faber 2007; Gayed & Bilello +2016 with leveraged ETFs); the volatility target scales exposure down when +QQQ is turbulent (Moreira & Muir 2017), which is when drawdowns cluster. + +Research (research/lab/vol_trend.py, synthetic 3x ETFs, one-day-late fills): +1990-2026 CAGR 16.0% vs QQQ 14.9%, MaxDD -37% vs -83%, alpha vs QQQ 6%/yr +(t 2.1). It lags QQQ in long calm bull markets (1990s, 2010-16) and earns its +keep in 2000-02, 2008 and 2022. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from strategies.base import Strategy + + +class QQQVolTrend(Strategy): + def __init__(self, ma: int = 200, target_vol: float = 0.30, cap: float = 1.5, + halflife: int = 20, signal: str = "QQQ", levered: str = "TQQQ", + leverage: float = 3.0, park: str = "SHY") -> None: + self.ma, self.target_vol, self.cap, self.halflife = ma, target_vol, cap, halflife + self.signal, self.levered, self.leverage, self.park = signal, levered, leverage, park + + def exposure(self, data: pd.DataFrame) -> pd.Series: + """Index exposure decided on each row's close (not yet shifted).""" + px = data[self.signal].ffill() + r = px.pct_change(fill_method=None) + sigma = np.sqrt((r ** 2).ewm(halflife=self.halflife, min_periods=40).mean() * 252) + on = (px > px.rolling(self.ma).mean()).astype(float) + e = (on * (self.target_vol / sigma).clip(upper=self.cap)).fillna(0.0) + e.iloc[:self.ma] = 0.0 + return e + + def generate_signals(self, data: pd.DataFrame) -> pd.DataFrame: + e = self.exposure(data) + w = pd.DataFrame(0.0, index=data.index, columns=data.columns) + lev = (e / self.leverage).where(data[self.levered].notna(), 0.0) + w[self.levered] = lev + if self.park in w.columns: + w[self.park] = (1.0 - lev).where(data[self.park].notna(), 0.0) + w.iloc[:self.ma] = 0.0 + return w.shift(1).fillna(0.0) + + +__all__ = ["QQQVolTrend"] diff --git a/strategies/sleeve_portfolio.py b/strategies/sleeve_portfolio.py index b26b5b4..295bcc2 100644 --- a/strategies/sleeve_portfolio.py +++ b/strategies/sleeve_portfolio.py @@ -14,6 +14,10 @@ Their monthly return correlations are 0.2-0.4, which is where the portfolio's higher Sharpe and smaller drawdown come from (2011-2026: Sharpe 1.21 vs 0.88-1.03 for the sleeves alone). +``us_momentum_qqq_trend`` is half overnight_recovery_top10, half the QQQ +trend of strategies/qqq_trend.py, from the alpha-vs-QQQ research +(research/lab/beat_qqq.py). + Each sleeve keeps its own strategy's targets. The stock sleeves see only the stock columns and the trend sleeve only its ETFs. A sleeve's share of the account drifts with its own return, like a separate account would, and all @@ -124,6 +128,29 @@ class SleevePortfolioStrategy(BarStrategy): return combined.shift(1).fillna(0.0) +MOMENTUM_QQQ_TREND_ETFS = ["QQQ", "SHY", "TQQQ"] + + +def us_momentum_qqq_trend() -> SleevePortfolioStrategy: + """Half overnight_recovery_top10, half QQQ trend (strategies/qqq_trend.py). + + From the 2026-09-29 alpha-vs-QQQ research (research/lab/beat_qqq.py): + momentum top 10 is the only stock selection whose alpha against QQQ held + in 2017-2026; the QQQ trend sleeve cuts the drawdowns. 2006-2026 alpha vs + QQQ +7.0%/yr (t 2.2), CAGR 17.4% vs 12.0% (2006-16) and 27.9% vs 21.5% + (2017-26), MaxDD -45%.""" + from strategies.qqq_trend import QQQVolTrend + from strategies.robust_momentum import RobustMomentumStrategy + return SleevePortfolioStrategy( + [ + ("momentum", RobustMomentumStrategy("overnight_recovery", top_n=10, rebal_freq=10, + buffer=2.0), 1.0, "stocks"), + ("qqq_trend", QQQVolTrend(), 1.0, "etfs"), + ], + etfs=MOMENTUM_QQQ_TREND_ETFS, + ) + + def us_three_sleeve() -> SleevePortfolioStrategy: from strategies.permanent import ETF_UNIVERSE from strategies.robust_momentum import RobustMomentumStrategy diff --git a/tests/test_sleeve_portfolio.py b/tests/test_sleeve_portfolio.py index 557d3a5..1d752f6 100644 --- a/tests/test_sleeve_portfolio.py +++ b/tests/test_sleeve_portfolio.py @@ -10,8 +10,9 @@ import pandas as pd import trader from strategies.permanent import ETF_UNIVERSE from strategies.robust_momentum import RobustMomentumStrategy +from strategies.qqq_trend import QQQVolTrend from strategies.sleeve_portfolio import (SleevePortfolioStrategy, decisions, drifting_shares, - us_three_sleeve) + us_momentum_qqq_trend, us_three_sleeve) def _bars(n_days: int = 520, n_stocks: int = 30, seed: int = 1) -> dict[str, pd.DataFrame]: @@ -89,5 +90,51 @@ class SleevePortfolioTests(unittest.TestCase): self.assertEqual(trader.strategy_fields("us_three_sleeve"), ("close", "open")) +def _qqq_bars(n_days: int = 520, drift: float = 0.0005, seed: int = 2) -> dict[str, pd.DataFrame]: + """Stocks plus the combo's ETFs only (its live universe: S&P 500 + QQQ, TQQQ, SHY).""" + other = [c for c in ETF_UNIVERSE if c not in ("TQQQ", "SHY")] + bars = {f: b.drop(columns=other) for f, b in _bars(n_days, seed=seed).items()} + rng = np.random.default_rng(seed + 1) + idx = bars["close"].index + q = pd.Series(100 * np.exp(np.cumsum(rng.normal(drift, 0.012, n_days))), index=idx) + for f in bars: + bars[f]["QQQ"] = q if f == "close" else q.shift(1).bfill() + return bars + + +class QQQTrendTests(unittest.TestCase): + def test_weights_on_a_day_ignore_that_day_and_later_closes(self): + close = _qqq_bars()["close"] + base = QQQVolTrend().generate_signals(close) + shocked = close.copy() + cut = close.index[400] + shocked.loc[cut:, "QQQ"] *= 0.5 + pd.testing.assert_frame_equal(base.loc[:cut], QQQVolTrend().generate_signals(shocked).loc[:cut]) + + def test_out_below_the_average_and_capped_above_it(self): + up = _qqq_bars(drift=0.002)["close"] + down = _qqq_bars(drift=-0.002)["close"] + w_up = QQQVolTrend(cap=1.5).generate_signals(up).iloc[260:] + w_down = QQQVolTrend().generate_signals(down).iloc[260:] + self.assertLessEqual(w_up["TQQQ"].max(), 0.5 + 1e-12) + self.assertGreater(w_up["TQQQ"].mean(), 0.3) + np.testing.assert_allclose(w_up[["TQQQ", "SHY"]].sum(axis=1), 1.0) + self.assertLess(w_down["TQQQ"].mean(), 0.05) + + def test_momentum_qqq_trend_is_a_us_mixed_strategy(self): + self.assertIsInstance(trader.STRATEGY_REGISTRY["us_momentum_qqq_trend"](), SleevePortfolioStrategy) + self.assertEqual(set(trader.strategy_extra_tickers("us_momentum_qqq_trend")), {"QQQ", "SHY", "TQQQ"}) + self.assertFalse(trader.strategy_supports_market("us_momentum_qqq_trend", "cn")) + bars = _qqq_bars() + strat = us_momentum_qqq_trend() + w = strat.generate_signals(bars) + self.assertLessEqual(w.sum(axis=1).max(), 1.0 + 1e-9) + stocks = [c for c in w.columns if c[1:].isdigit()] + self.assertLessEqual(int((w[stocks].iloc[-1] > 0).sum()), 10) + self.assertGreater(w[["TQQQ", "SHY"]].iloc[-1].sum(), 0.3) + days = bars["close"].index[-2:] + pd.testing.assert_frame_equal(trader.bar_session_weights(strat, bars, days), w.loc[days]) + + if __name__ == "__main__": unittest.main() diff --git a/trader.py b/trader.py index 3985b71..d4dd59e 100644 --- a/trader.py +++ b/trader.py @@ -58,7 +58,8 @@ from strategies.permanent import ( ) from strategies.recovery_momentum import RecoveryMomentumStrategy from strategies.robust_momentum import RobustMomentumStrategy -from strategies.sleeve_portfolio import us_three_sleeve +from strategies.sleeve_portfolio import (MOMENTUM_QQQ_TREND_ETFS, us_momentum_qqq_trend, + us_three_sleeve) from strategies.trend_following import TrendFollowingStrategy from strategies.trend_rider_v5 import TrendRiderV5 from strategies.trend_rider_voltgt import ( @@ -298,6 +299,8 @@ STRATEGY_REGISTRY = { # 1/3 each of overnight_recovery_top10, lowrisk_top20 and trend_rider_v7_vt24 # in one account (strategies/sleeve_portfolio.py) "us_three_sleeve": lambda **kw: us_three_sleeve(), + # half overnight_recovery_top10, half QQQ 200-day trend with a 30% vol target (TQQQ + SHY) + "us_momentum_qqq_trend": lambda **kw: us_momentum_qqq_trend(), # --- TrendRiderV6: stock-picking + V5 regime engine --- "trend_rider_v6": lambda **kw: TrendRiderV6(), "trend_rider_v6_top10": lambda **kw: TrendRiderV6(top_n=10), @@ -328,6 +331,7 @@ MIXED_STRATEGY_EXTRA_TICKERS = { "trend_rider_v6": sorted(set(ETF_UNIVERSE)), "trend_rider_v6_top10": sorted(set(ETF_UNIVERSE)), "us_three_sleeve": sorted(set(ETF_UNIVERSE)), + "us_momentum_qqq_trend": MOMENTUM_QQQ_TREND_ETFS, } # Registered strategies the live monitor does not run, with the reason.