diff --git a/.gitignore b/.gitignore index c465801..8127e83 100644 --- a/.gitignore +++ b/.gitignore @@ -59,3 +59,6 @@ memory/ # Web UI state and precomputed backtests (web/src/server/paths.ts) data/web/ data/precomputed/ + +# Research panels and results — fetched/rebuilt by research/lab +data/research/ diff --git a/research/lab/__init__.py b/research/lab/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/research/lab/baselines.py b/research/lab/baselines.py new file mode 100644 index 0000000..db970ec --- /dev/null +++ b/research/lab/baselines.py @@ -0,0 +1,47 @@ +""" +Re-run the existing factor_combo champions under the research engine, with +the same rules every new candidate gets: next-open execution, PIT universe, +costs, and averaging over rebalance phases. + + PYTHONPATH=. uv run python -m research.lab.baselines --market us +""" + +import argparse +import warnings + +import pandas as pd + +from research.lab.factor_scan import load, phase_avg +from research.lab.report import fmt +from strategies.factor_combo import SIGNAL_REGISTRY + +warnings.filterwarnings("ignore") + +CHAMPIONS = { + "us": ["rec_mfilt+deep_upvol", "ma200+mom7m+rec126", "rec_mfilt+ma200", "mom7m+rec126"], + "cn": ["up_cap+quality_mom", "down_resil+qual_mom", "rec63+mom_gap", "up_cap+mom_gap"], +} + + +def baseline_scores(p, market: str) -> dict[str, pd.DataFrame]: + mk = "us" if market == "us" else "cn" + # factor_combo signals were built on prices masked to index membership + px = p.close.where(p.member).ffill(limit=5) + return {name: SIGNAL_REGISTRY[name](px) for name in CHAMPIONS[mk]} + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + args = ap.parse_args() + p, periods = load(args.market) + scores = baseline_scores(p, args.market) + for name, s in scores.items(): + for every, label in [(5, "weekly"), (21, "monthly")]: + for per, (a, b) in periods.items(): + st = phase_avg(p, s, a, b, 10, every) + print(f"{name:22s} top10 {label:8s} {per:8s} {fmt(st)}", flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/construction.py b/research/lab/construction.py new file mode 100644 index 0000000..bf1901f --- /dev/null +++ b/research/lab/construction.py @@ -0,0 +1,61 @@ +""" +Construction sweep: holdings count, buffer (hysteresis), rebalance interval and +weighting for a fixed score, reported per ~5-year block. + + PYTHONPATH=. uv run python -m research.lab.construction --market us --combo overnight_mom+rec_126 +""" + +import argparse +import itertools +import warnings + +import numpy as np +import pandas as pd + +from research.lab.baselines import baseline_scores +from research.lab.factor_scan import load, phase_avg +from research.lab.factors import FACTORS, zrank +from research.lab.stepwise import SUBPERIODS + +warnings.filterwarnings("ignore") + + +def combo_score(p, combo: str, market: str) -> pd.DataFrame: + if combo.startswith("baseline:"): + return baseline_scores(p, market)[combo.split(":", 1)[1]] + parts = combo.split("+") + return sum(zrank(FACTORS[n](p), p.member) for n in parts) / len(parts) + + +def block_stats(p, score, subs, **kw) -> dict: + per = [phase_avg(p, score, a, b, **kw) for a, b in subs] + return {"sharpe": np.mean([x["sharpe"] for x in per]), "min_sh": min(x["sharpe"] for x in per), + "cagr": np.mean([x["cagr"] for x in per]), "mdd": min(x["mdd"] for x in per), + "turn": np.mean([x["turnover"] for x in per]), + "blocks": " ".join(f"{x['sharpe']:.2f}/{x['cagr']*100:.0f}%" for x in per)} + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + ap.add_argument("--combo", nargs="+", required=True) + ap.add_argument("--tops", nargs="*", type=int, default=[10, 15, 20, 30]) + ap.add_argument("--buffers", nargs="*", type=float, default=[0, 0.5, 1.0, 2.0]) + ap.add_argument("--everys", nargs="*", type=int, default=[5, 10]) + ap.add_argument("--weightings", nargs="*", default=["equal", "invvol"]) + args = ap.parse_args() + mk = "us" if args.market == "us" else "cn" + p, _ = load(args.market) + subs = SUBPERIODS[mk] + for combo in args.combo: + s = combo_score(p, combo, args.market) + print(f"=== {combo}") + for top, buf, every, wt in itertools.product(args.tops, args.buffers, args.everys, args.weightings): + st = block_stats(p, s, subs, top_n=top, every=every, buffer=buf, weighting=wt) + print(f" top{top:<3d} buf{buf:<4} every{every:<3d} {wt:7s} sh {st['sharpe']:.2f} min {st['min_sh']:.2f} " + f"cagr {st['cagr']*100:5.1f}% worstMDD {st['mdd']*100:5.1f}% turn {st['turn']:4.1f}x [{st['blocks']}]", + flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/engine.py b/research/lab/engine.py new file mode 100644 index 0000000..6a8891d --- /dev/null +++ b/research/lab/engine.py @@ -0,0 +1,252 @@ +""" +Event-driven long-only portfolio simulator for factor research. + +Differences from ``main.backtest`` (the engine behind the web UI numbers): + +* Execution lag. A score computed from the close of day s is traded at the + OPEN of day s+1 (``exec_at="open"``) or the CLOSE of day s+1 + (``exec_at="close"``). ``main.backtest`` fills at the same close the signal + saw, which the live trader cannot do. +* Holdings drift between rebalances; costs are charged on actual trades. +* Tradability. A stock that is suspended, or opens at limit-up (buy) or + limit-down (sell) on the execution day keeps its current position. Stocks + outside the point-in-time index on the signal date are never bought. +* Costs are asymmetric (CN stamp tax on sells) and include an optional fixed + per-order fee expressed against a reference capital. + +Conventions: every DataFrame shares one date index and one column set +(``Panel.close.columns``). ``score`` is "higher is better"; NaN means +"not eligible". +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd + + +@dataclass +class Costs: + buy: float = 0.0010 # fraction of traded value + sell: float = 0.0010 + sell_tax: pd.Series | None = None # extra sell-side tax by date (CN stamp duty) + fixed_fee: float = 0.0 # currency per order + capital: float = 10_000.0 # reference capital for the fixed fee + + +@dataclass +class Panel: + market: str + open: pd.DataFrame + high: pd.DataFrame + low: pd.DataFrame + close: pd.DataFrame + volume: pd.DataFrame + member: pd.DataFrame # bool: in PIT index on that date + can_buy: pd.DataFrame # bool: buyable at that day's open + can_sell: pd.DataFrame # bool: sellable at that day's open + cash_ret: pd.Series # daily return of the idle-cash asset + benchmark: pd.Series # index level (cap-weighted) + costs: Costs = field(default_factory=Costs) + extra: dict = field(default_factory=dict) # market-specific fields (turn, pe, ...) + + @property + def dates(self) -> pd.DatetimeIndex: + return self.close.index + + def slice(self, start=None, end=None) -> "Panel": + idx = self.dates + m = np.ones(len(idx), bool) + if start is not None: + m &= idx >= pd.Timestamp(start) + if end is not None: + m &= idx <= pd.Timestamp(end) + def sub(x): + if isinstance(x, (pd.Series, pd.DataFrame)) and x.index.equals(idx): + return x.loc[m] + return x + return Panel(self.market, sub(self.open), sub(self.high), sub(self.low), sub(self.close), + sub(self.volume), sub(self.member), sub(self.can_buy), sub(self.can_sell), + sub(self.cash_ret), sub(self.benchmark), self.costs, + {k: sub(v) for k, v in self.extra.items() if not k.startswith("_")}) + + +@dataclass +class Result: + equity: pd.Series + returns: pd.Series + turnover: pd.Series # one-way traded value / equity, per day + holdings: pd.DataFrame # weights after each rebalance (rows = exec dates) + exposure: pd.Series # invested fraction at each close + name: str = "" + + +def rebalance_dates(dates: pd.DatetimeIndex, every: int | str, offset: int = 0, + start: int = 0) -> np.ndarray: + """Boolean mask of signal dates. + + ``every`` is a trading-day interval (5 = weekly) or "W"/"M" for the last + trading day of each calendar week/month. + """ + mask = np.zeros(len(dates), bool) + if isinstance(every, str): + key = dates.to_period(every[0]) + last = pd.Series(np.arange(len(dates)), index=dates).groupby(key).max().values + mask[last] = True + mask[:start] = False + return mask + mask[start + offset::every] = True + return mask + + +def select_top(score_row: np.ndarray, top_n: int, held: np.ndarray | None = None, + buffer: float = 0.0) -> np.ndarray: + """Indices of the chosen names. With a buffer, current holdings stay while + their rank is within ``top_n * (1 + buffer)``.""" + valid = np.where(~np.isnan(score_row))[0] + if len(valid) == 0: + return valid + order = valid[np.argsort(-score_row[valid], kind="stable")] + if buffer <= 0 or held is None or not held.any(): + return order[:top_n] + rank = np.empty(len(score_row)); rank[:] = np.inf + rank[order] = np.arange(len(order)) + keep = [i for i in np.where(held)[0] if rank[i] < top_n * (1 + buffer)] + keep = sorted(keep, key=lambda i: rank[i])[:top_n] + chosen = list(keep) + for i in order: + if len(chosen) >= top_n: + break + if i not in keep: + chosen.append(i) + return np.array(chosen, int) + + +def simulate(panel: Panel, score: pd.DataFrame | None = None, top_n: int = 10, + every: int | str = 5, offset: int = 0, warmup: int = 0, + weighting: str = "equal", buffer: float = 0.0, + exposure: pd.Series | None = None, exec_at: str = "open", + target_weights: pd.DataFrame | None = None, min_trade: float = 0.002, + name: str = "") -> Result: + """Run a long-only top-N strategy. + + score : higher = better, evaluated at each signal date's close. + target_weights : alternative to score — explicit weights per signal date + (rows need only exist on signal dates; NaN rows skip). + weighting : "equal" | "invvol" (60-day vol) | "score" (score-proportional, clipped >0) + exposure : fraction of equity to invest, read on the signal date + (e.g. a market-trend or vol-target overlay); rest earns cash_ret. + exec_at : "open" → trade at next day's open; "close" → next day's close. + min_trade : resizing trades smaller than this fraction of equity are + skipped (the live trader skips tiny orders too); full + entries and exits always trade. + """ + dates = panel.dates + n_days, n_assets = panel.close.shape + C = panel.close.ffill().values + O = panel.open.values + O = np.where(np.isnan(O), C, O) # missing open → fall back to close + prevC = np.vstack([C[:1], C[:-1]]) + member = panel.member.values + can_buy = panel.can_buy.values + can_sell = panel.can_sell.values + cash_ret = panel.cash_ret.reindex(dates).fillna(0.0).values + sell_tax = (panel.costs.sell_tax.reindex(dates).fillna(0.0).values + if panel.costs.sell_tax is not None else np.zeros(n_days)) + S = score.reindex(index=dates, columns=panel.close.columns).values if score is not None else None + TW = (target_weights.reindex(index=dates, columns=panel.close.columns).values + if target_weights is not None else None) + expo = exposure.reindex(dates).ffill().fillna(1.0).values if exposure is not None else np.ones(n_days) + if weighting == "invvol": + rets = panel.close.pct_change(fill_method=None) + vol = rets.rolling(60, min_periods=40).std().values + sig_mask = rebalance_dates(dates, every, offset, start=warmup) + + h = np.zeros(n_assets) # value held per asset (at last mark) + cash = 1.0 + eq = np.empty(n_days); turn = np.zeros(n_days); inv = np.zeros(n_days) + hold_rows, hold_idx = [], [] + pending = None # target weight vector awaiting execution + + for d in range(n_days): + # mark overnight / full-day moves + with np.errstate(invalid="ignore", divide="ignore"): + gap = np.where(prevC[d] > 0, O[d] / prevC[d], 1.0) + intraday = np.where(O[d] > 0, C[d] / O[d], 1.0) + gap = np.nan_to_num(gap, nan=1.0, posinf=1.0) + intraday = np.nan_to_num(intraday, nan=1.0, posinf=1.0) + cash *= 1 + cash_ret[d] + + if pending is not None: + if exec_at == "open": + h *= gap + else: + h *= gap * intraday + V = h.sum() + cash + tgt = pending * V + # untradable names keep their current value + frozen_buy = (tgt > h) & ~can_buy[d] + frozen_sell = (tgt < h) & ~can_sell[d] + frozen = frozen_buy | frozen_sell + if frozen.any(): + tgt[frozen] = h[frozen] + free = ~frozen & (pending > 0) + room = V - tgt[frozen].sum() + if free.any() and pending[free].sum() > 0: + budget = min(room, pending[free].sum() * V) + tgt[free] = pending[free] / pending[free].sum() * max(budget, 0.0) + if min_trade > 0: + tiny = (np.abs(tgt - h) < min_trade * V) & (tgt > 0) & (h > 0) + tgt[tiny] = h[tiny] + trade = tgt - h + buys = trade[trade > 0].sum() + sells = -trade[trade < 0].sum() + n_orders = int((np.abs(trade) > 1e-4 * V).sum()) + cost = buys * panel.costs.buy + sells * (panel.costs.sell + sell_tax[d]) + # Equity is normalized to 1.0 == costs.capital, so a fixed fee per + # order is fixed_fee / capital in these units at any book size. + if panel.costs.fixed_fee > 0: + cost += n_orders * panel.costs.fixed_fee / panel.costs.capital + cash += h.sum() - tgt.sum() - cost + h = tgt + turn[d] = (buys + sells) / 2 / V if V > 0 else 0.0 + hold_rows.append(h / max(h.sum() + cash, 1e-12)); hold_idx.append(dates[d]) + if exec_at == "open": + h *= intraday + pending = None + else: + h *= gap * intraday + + V = h.sum() + cash + eq[d] = V + inv[d] = h.sum() / V if V > 0 else 0.0 + + if sig_mask[d] and d + 1 < n_days: + if TW is not None: + w = np.nan_to_num(TW[d], nan=0.0) + if np.isnan(TW[d]).all(): + continue + else: + row = np.where(member[d], S[d], np.nan) + held = h > 1e-12 + idx = select_top(row, top_n, held, buffer) + w = np.zeros(n_assets) + if len(idx) >= min(max(1, top_n // 2), 20): + if weighting == "invvol": + iv = 1.0 / np.where(vol[d, idx] > 0, vol[d, idx], np.nan) + iv = np.nan_to_num(iv, nan=np.nanmean(iv) if np.isfinite(iv).any() else 1.0) + w[idx] = iv / iv.sum() + elif weighting == "score": + s = np.clip(row[idx] - np.nanmin(row[idx]) + 1e-6, 0, None) + w[idx] = s / s.sum() + else: + w[idx] = 1.0 / len(idx) + pending = w * float(np.clip(expo[d], 0.0, 1.0)) + + equity = pd.Series(eq, index=dates) + rets = equity.pct_change().fillna(0.0) + holdings = pd.DataFrame(hold_rows, index=hold_idx, columns=panel.close.columns) + return Result(equity, rets, pd.Series(turn, index=dates), holdings, + pd.Series(inv, index=dates), name) diff --git a/research/lab/evaluate.py b/research/lab/evaluate.py new file mode 100644 index 0000000..512c6b2 --- /dev/null +++ b/research/lab/evaluate.py @@ -0,0 +1,143 @@ +""" +Final evaluation of candidate strategies against the existing champions, +under identical rules (PIT universe, next-open execution, costs, phase +averaging). Writes data/research/results/final_{market}.json. + + PYTHONPATH=. uv run python -m research.lab.evaluate --market us +""" + +import argparse +import json +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab import overlays as ov +from research.lab.construction import combo_score +from research.lab.engine import Costs, simulate +from research.lab.factor_scan import RESULTS, load +from research.lab.report import deflated_sharpe, stats, yearly +from research.lab.stepwise import SUBPERIODS + +warnings.filterwarnings("ignore") + +# name -> (combo, construction). "baseline:" combos are the current factor_combo champions. +CANDIDATES = { + "us": { + "champion rec_mfilt+deep_upvol monthly": ("baseline:rec_mfilt+deep_upvol", dict(top_n=10, every=21)), + "champion mom7m+rec126 monthly": ("baseline:mom7m+rec126", dict(top_n=10, every=21)), + "champion mom7m+rec126 weekly": ("baseline:mom7m+rec126", dict(top_n=10, every=5)), + "EW S&P 500 (PIT)": ("ew", dict(top_n=10000, every=5)), + "new overnight+rec": ("overnight_mom+rec_126", dict(top_n=10, every=10, buffer=2.0)), + "new overnight6+rec": ("overnight_mom_6+rec_126", dict(top_n=10, every=10, buffer=2.0)), + "new smooth-mom4": ("overnight_mom+overnight_mom_6+amp_cut_mom+rec_126", dict(top_n=10, every=10, buffer=2.0)), + "new overnight+ampcut+rec": ("overnight_mom+amp_cut_mom+rec_126", dict(top_n=10, every=10, buffer=2.0)), + }, + "cn": { + "champion up_cap+quality_mom weekly": ("baseline:up_cap+quality_mom", dict(top_n=10, every=5)), + "champion up_cap+mom_gap weekly": ("baseline:up_cap+mom_gap", dict(top_n=10, every=5)), + "champion down_resil+qual_mom biweekly": ("baseline:down_resil+qual_mom", dict(top_n=10, every=10)), + "EW CSI 300 (PIT)": ("ew", dict(top_n=10000, every=5)), + }, +} +LIVE_FEES = {"us": (2.0, 10_000.0), "cn": (5.0, 70_000.0)} +WINDOW = ("2016-09-27", "2026-09-24") # the web UI's 10-year window + + +def score_for(p, combo, market): + if combo == "ew": + return p.close.notna().astype(float).where(p.member) + return combo_score(p, combo, market) + + +def phase_runs(p, score, a, b, exposure=None, **kw): + every = kw.get("every", 5) + w0 = p.dates.searchsorted(pd.Timestamp(a)) + sp = p.slice(p.dates[max(w0 - 5, 0)], b) + runs = [] + for off in range(every if kw.get("top_n", 10) < 1000 else 1): + r = simulate(sp, score, offset=off, warmup=5, exposure=exposure, **kw) + runs.append(r) + return runs + + +def summarize(runs, a, b) -> dict: + sts = [stats(r.equity.loc[a:b], r.turnover) for r in runs] + df = pd.DataFrame(sts) + rets = pd.concat([r.returns.loc[a:b] for r in runs], axis=1).mean(axis=1) + eq = (1 + rets).cumprod() + out = {k: float(df[k].mean()) for k in df.columns} + out["cagr_phase_min"] = float(df["cagr"].min()) + out["cagr_phase_max"] = float(df["cagr"].max()) + out["yearly"] = {str(k.year): float(v) for k, v in yearly(eq).items()} + out["_rets"] = rets + return out + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + ap.add_argument("--add", nargs="*", default=[], help="extra 'name=combo|top|every|buffer'") + ap.add_argument("--n-trials", type=int, default=0, help="override trial count for DSR") + args = ap.parse_args() + mk = "us" if args.market == "us" else "cn" + p, _ = load(args.market) + cands = dict(CANDIDATES[mk]) + for spec in args.add: + name, rest = spec.split("=", 1) + combo, top, every, buf = rest.split("|") + cands[name] = (combo, dict(top_n=int(top), every=int(every), buffer=float(buf))) + trials_path = os.path.join(RESULTS, f"trials_{args.market}.csv") + n_trials = args.n_trials or (len(pd.read_csv(trials_path)) if os.path.exists(trials_path) else 100) + full = (SUBPERIODS[mk][0][0], WINDOW[1]) + out = {} + fee, cap = LIVE_FEES[mk] + for name, (combo, kw) in cands.items(): + s = score_for(p, combo, args.market) + res = {"combo": combo, "construction": kw} + base_runs = phase_runs(p, s, *WINDOW, **kw) + res["window"] = summarize(base_runs, *WINDOW) + res["full"] = summarize(phase_runs(p, s, *full, **kw), *full) + res["blocks"] = [summarize(phase_runs(p, s, a, b, **kw), a, b) for a, b in SUBPERIODS[mk]] + # live-book fixed fee + pf = p.slice(); pf.costs = Costs(p.costs.buy, p.costs.sell, p.costs.sell_tax, fee, cap) + res["window_livefee"] = summarize(phase_runs(pf, s, *WINDOW, **kw), *WINDOW) + # double proportional cost + p2 = p.slice(); p2.costs = Costs(p.costs.buy * 2, p.costs.sell * 2, p.costs.sell_tax, 0.0, cap) + res["window_2xcost"] = summarize(phase_runs(p2, s, *WINDOW, **kw), *WINDOW) + res["window_exec_close"] = summarize(phase_runs(p, s, *WINDOW, exec_at="close", **kw), *WINDOW) + # drawdown-control overlay: slow/fast market trend blend on the EW index + ew = ov.ew_index(p) + slow = (ew / ew.shift(252) > 1).astype(float) + fast = (ew / ew.shift(21) > 1).astype(float) + expo = (0.5 + 0.25 * slow + 0.25 * fast).where(slow + fast < 2, 1.0) + res["window_trend_overlay"] = summarize(phase_runs(p, s, *WINDOW, exposure=expo, **kw), *WINDOW) + dsr = deflated_sharpe(res["full"]["_rets"], n_trials) + res["dsr_full"] = dsr + out[name] = res + w, f = res["window"], res["full"] + print(f"{name:42s} 2016-26: CAGR {w['cagr']*100:5.1f}% [{w['cagr_phase_min']*100:4.1f},{w['cagr_phase_max']*100:4.1f}] " + f"Sh {w['sharpe']:.2f} MDD {w['mdd']*100:5.1f}% | full: CAGR {f['cagr']*100:5.1f}% Sh {f['sharpe']:.2f} " + f"MDD {f['mdd']*100:5.1f}% | blocks Sh {' '.join(f'{b_['sharpe']:.2f}' for b_ in res['blocks'])} " + f"| livefee {res['window_livefee']['cagr']*100:5.1f}% 2x {res['window_2xcost']['cagr']*100:5.1f}% " + f"close {res['window_exec_close']['cagr']*100:5.1f}% trend {res['window_trend_overlay']['cagr']*100:5.1f}%/" + f"{res['window_trend_overlay']['sharpe']:.2f}/{res['window_trend_overlay']['mdd']*100:.0f}% " + f"| DSR {dsr['dsr']:.2f} (N={n_trials})", flush=True) + + def clean(o): + if isinstance(o, dict): + return {k: clean(v) for k, v in o.items() if not k.startswith("_")} + if isinstance(o, list): + return [clean(v) for v in o] + if isinstance(o, (np.floating, np.integer)): + return float(o) + return o + os.makedirs(RESULTS, exist_ok=True) + with open(os.path.join(RESULTS, f"final_{args.market}.json"), "w") as fh: + json.dump(clean(out), fh, indent=1) + + +if __name__ == "__main__": + main() diff --git a/research/lab/export_curves.py b/research/lab/export_curves.py new file mode 100644 index 0000000..2561e05 --- /dev/null +++ b/research/lab/export_curves.py @@ -0,0 +1,78 @@ +""" +Export phase-averaged equity curves of the final candidates and baselines for +the report (weekly samples), plus the US/CN correlation of the new strategies. + + PYTHONPATH=. uv run python -m research.lab.export_curves +""" + +import json +import os +import warnings + +import numpy as np +import pandas as pd + +from research.lab.evaluate import phase_runs, score_for +from research.lab.factor_scan import RESULTS, load +from research.lab.report import stats + +warnings.filterwarnings("ignore") + +SETS = { + "us": ("2006-09-27", "2026-09-24", { + "新: overnight+recovery": ("overnight_mom+rec_126", dict(top_n=10, every=10, buffer=2.0)), + "现: rec_mfilt+deep_upvol 月频": ("baseline:rec_mfilt+deep_upvol", dict(top_n=10, every=21)), + "现: mom7m+rec126 月频": ("baseline:mom7m+rec126", dict(top_n=10, every=21)), + "等权 S&P 500 (PIT)": ("ew", dict(top_n=10000, every=5)), + }), + "cn_hs300": ("2010-01-04", "2026-09-24", { + "新: 低风险动量4因子": ("amp_cut_mom+low_beta+pv_corr+sharpe_mom", dict(top_n=15, every=10, buffer=2.0)), + "新: low_beta+sharpe_mom": ("low_beta+sharpe_mom", dict(top_n=20, every=10, buffer=1.0)), + "现: down_resil+qual_mom 双周": ("baseline:down_resil+qual_mom", dict(top_n=10, every=10)), + "现: up_cap+quality_mom 周频": ("baseline:up_cap+quality_mom", dict(top_n=10, every=5)), + "等权 CSI 300 (PIT)": ("ew", dict(top_n=10000, every=5)), + }), +} + + +def main() -> None: + out, daily = {}, {} + for market, (a, b, cands) in SETS.items(): + p, _ = load(market) + curves = {} + for name, (combo, kw) in cands.items(): + s = score_for(p, combo, market) + runs = phase_runs(p, s, a, b, **kw) + rets = pd.concat([r.returns.loc[a:b] for r in runs], axis=1).mean(axis=1) + daily[(market, name)] = rets + eq = (1 + rets).cumprod() + wk = eq.resample("W-FRI").last().dropna() + dd = eq / eq.cummax() - 1 + curves[name] = {"dates": [d.strftime("%Y-%m-%d") for d in wk.index], + "equity": [round(float(v), 4) for v in wk.values], + "drawdown": [round(float(v), 4) for v in dd.resample("W-FRI").min().dropna().values], + "stats": {k: float(v) for k, v in stats(eq).items()}} + print(market, name, {k: round(v, 3) for k, v in curves[name]["stats"].items()}, flush=True) + out[market] = curves + us = daily[("us", "新: overnight+recovery")] + cn = daily[("cn_hs300", "新: 低风险动量4因子")] + # compound each market on its own calendar before joining: the exchanges + # have different holidays, and an inner join on days would drop returns + m = pd.concat([us, cn], axis=1, keys=["us", "cn"]).dropna().loc["2016-09-27":] + mon = pd.concat({k: (1 + v.loc["2016-09-27":]).resample("ME").prod() - 1 + for k, v in (("us", us), ("cn", cn))}, axis=1).dropna() + cross = {"daily_corr": float(m.corr().iloc[0, 1]), "monthly_corr": float(mon.corr().iloc[0, 1])} + for col in ["us", "cn", "5050"]: + r = mon.mean(axis=1) if col == "5050" else mon[col] + eq = (1 + r).cumprod() + cross[col] = {"cagr": float(eq.iloc[-1] ** (12 / len(r)) - 1), + "sharpe_monthly": float(r.mean() / r.std() * np.sqrt(12)), + "mdd_monthly": float((eq / eq.cummax() - 1).min())} + out["cross"] = cross + print("cross", out["cross"]) + with open(os.path.join(RESULTS, "curves.json"), "w") as fh: + json.dump(out, fh, ensure_ascii=False, default=float) + + +if __name__ == "__main__": + main() diff --git a/research/lab/factor_scan.py b/research/lab/factor_scan.py new file mode 100644 index 0000000..e6568ae --- /dev/null +++ b/research/lab/factor_scan.py @@ -0,0 +1,95 @@ +""" +Single-factor scan: rank IC and long-only top-N backtests per period. + + PYTHONPATH=. uv run python -m research.lab.factor_scan --market us + PYTHONPATH=. uv run python -m research.lab.factor_scan --market cn_hs300 + +Writes data/research/results/factor_scan_{market}.csv. +""" + +import argparse +import os +import time +import warnings + +import pandas as pd + +from research.lab import panel as pnl +from research.lab.engine import simulate +from research.lab.factors import FACTORS +from research.lab.report import (forward_returns, ic_series, ic_summary, stats, + top_minus_universe) + +warnings.filterwarnings("ignore") + +PERIODS = { + "us": {"design": ("2006-01-01", "2016-09-26"), "holdout": ("2016-09-27", "2026-09-24")}, + "cn": {"design": ("2010-01-01", "2016-09-26"), "holdout": ("2016-09-27", "2026-09-24")}, +} +RESULTS = os.path.join("data", "research", "results") + + +def load(market: str): + if market == "us": + return pnl.load_us(), PERIODS["us"] + return pnl.load_cn(market.split("_", 1)[1]), PERIODS["cn"] + + +def phase_avg(p, score, start, end, top_n, every, **kw) -> dict: + """Average stats over every rebalance phase (offset 0..every-1).""" + w0 = p.dates.searchsorted(pd.Timestamp(start)) + sp = p.slice(p.dates[max(w0 - 5, 0)], end) + rows = [] + for off in range(every): + r = simulate(sp, score, top_n=top_n, every=every, offset=off, warmup=5, **kw) + rows.append(stats(r.equity.loc[start:end], r.turnover)) + return pd.DataFrame(rows).mean().to_dict() + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + ap.add_argument("--top-n", type=int, default=20) + ap.add_argument("--every", type=int, default=5) + ap.add_argument("--only", nargs="*") + args = ap.parse_args() + mk = "us" if args.market == "us" else "cn" + p, periods = load(args.market) + fwd5 = forward_returns(p, 5) + fwd20 = forward_returns(p, 20) + rows = [] + names = args.only or [n for n, f in FACTORS.items() if mk in f.markets] + for name in names: + t = time.time() + try: + s = FACTORS[name](p) + except Exception as exc: # noqa: BLE001 + print(f"{name}: failed {exc}") + continue + row = {"factor": name} + for per, (a, b) in periods.items(): + sl = slice(a, b) + ic5 = ic_summary(ic_series(s.loc[sl], fwd5.loc[sl], p.member.loc[sl], 5)) + ic20 = ic_summary(ic_series(s.loc[sl], fwd20.loc[sl], p.member.loc[sl], 20), 12) + tmu = top_minus_universe(s.loc[sl], fwd5.loc[sl], p.member.loc[sl], args.top_n, 5) + st = phase_avg(p, s, a, b, args.top_n, args.every) + row.update({f"{per}_ic5": ic5["ic"], f"{per}_ic5_t": ic5["t"], + f"{per}_ic20": ic20["ic"], f"{per}_ic20_t": ic20["t"], + f"{per}_top_excess_ann": tmu.mean() * 52, + f"{per}_cagr": st["cagr"], f"{per}_sharpe": st["sharpe"], + f"{per}_mdd": st["mdd"], f"{per}_turn": st["turnover"]}) + rows.append(row) + print(f"{name:16s} " + f"D: ic5 {row['design_ic5']:+.3f} (t{row['design_ic5_t']:+.1f}) ic20 {row['design_ic20']:+.3f} " + f"xs {row['design_top_excess_ann']*100:+5.1f}% cagr {row['design_cagr']*100:5.1f}% sh {row['design_sharpe']:.2f} | " + f"H: ic5 {row['holdout_ic5']:+.3f} (t{row['holdout_ic5_t']:+.1f}) ic20 {row['holdout_ic20']:+.3f} " + f"xs {row['holdout_top_excess_ann']*100:+5.1f}% cagr {row['holdout_cagr']*100:5.1f}% sh {row['holdout_sharpe']:.2f} " + f"turn {row['holdout_turn']:.0f}x [{time.time()-t:.0f}s]", flush=True) + # equal-weight universe reference + os.makedirs(RESULTS, exist_ok=True) + out = pd.DataFrame(rows) + out.to_csv(os.path.join(RESULTS, f"factor_scan_{args.market}.csv"), index=False) + + +if __name__ == "__main__": + main() diff --git a/research/lab/factors.py b/research/lab/factors.py new file mode 100644 index 0000000..2f10ab3 --- /dev/null +++ b/research/lab/factors.py @@ -0,0 +1,529 @@ +""" +Factor library. Every factor maps a Panel to a DataFrame (dates x assets) +where HIGHER = MORE ATTRACTIVE and the value on date t uses data up to the +close of t only. Sign conventions follow the literature: e.g. ``ivol`` returns +NEGATIVE idiosyncratic volatility because low-ivol stocks earn more. + +Windows are in trading days: 21 ≈ 1 month, 63 ≈ 1 quarter, 252 ≈ 1 year. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from research.lab.engine import Panel + +FACTORS: dict[str, callable] = {} + + +def factor(name: str, *, markets: tuple[str, ...] = ("us", "cn")): + def deco(fn): + fn.markets = markets + FACTORS[name] = fn + return fn + return deco + + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + +def _ret(p: Panel) -> pd.DataFrame: + key = "_ret" + if key not in p.extra: + p.extra[key] = p.close.pct_change(fill_method=None) + return p.extra[key] + + +def _mkt(p: Panel) -> pd.Series: + """Equal-weight return of PIT members (the investable universe).""" + key = "_mkt" + if key not in p.extra: + p.extra[key] = _ret(p).where(p.member).mean(axis=1).fillna(0.0) + return p.extra[key] + + +def _rolling_beta(r: pd.DataFrame, m: pd.Series, win: int, minp: int) -> pd.DataFrame: + mr = r.mul(m, axis=0) + e_xy = mr.rolling(win, min_periods=minp).mean() + e_x = r.rolling(win, min_periods=minp).mean() + e_y = m.rolling(win, min_periods=minp).mean() + var_y = (m ** 2).rolling(win, min_periods=minp).mean() - e_y ** 2 + return (e_xy - e_x.mul(e_y, axis=0)).div(var_y, axis=0) + + +def resid_returns(p: Panel, win: int = 252) -> pd.DataFrame: + key = f"_resid{win}" + if key not in p.extra: + r = _ret(p) + m = _mkt(p) + beta = _rolling_beta(r, m, win, win // 2).shift(1) + p.extra[key] = r - beta.mul(m, axis=0) + return p.extra[key] + + +def sector_map(p: Panel) -> pd.Series | None: + return p.extra.get("sector") + + +def demean_by_sector(df: pd.DataFrame, p: Panel) -> pd.DataFrame: + sec = sector_map(p) + if sec is None: + return df.sub(df.where(p.member).mean(axis=1), axis=0) + out = df.copy() + m = df.where(p.member) + for s in sec.dropna().unique(): + cols = [c for c in sec.index[sec == s] if c in df.columns] + if cols: + out[cols] = df[cols].sub(m[cols].mean(axis=1), axis=0) + return out + + +# --------------------------------------------------------------------------- +# momentum family +# --------------------------------------------------------------------------- + +@factor("mom_12_1") +def mom_12_1(p): + c = p.close.ffill() + return c.shift(21) / c.shift(252) - 1 + + +@factor("mom_6_1") +def mom_6_1(p): + c = p.close.ffill() + return c.shift(21) / c.shift(126) - 1 + + +@factor("mom_12_7") # Novy-Marx (2012) intermediate momentum +def mom_12_7(p): + c = p.close.ffill() + return c.shift(147) / c.shift(252) - 1 + + +@factor("mom_3") +def mom_3(p): + c = p.close.ffill() + return c / c.shift(63) - 1 + + +@factor("resid_mom") # Blitz, Huij & Martens (2011): residual 12-1, vol-scaled +def resid_mom(p): + e = resid_returns(p) + s = e.shift(21).rolling(231, min_periods=180) + return s.sum() / s.std() + + +@factor("resid_mom_6") +def resid_mom_6(p): + e = resid_returns(p) + s = e.shift(21).rolling(105, min_periods=80) + return s.sum() / s.std() + + +@factor("sharpe_mom") # 12-1 return / 12-month vol +def sharpe_mom(p): + r = _ret(p) + return mom_12_1(p) / r.rolling(252, min_periods=180).std() + + +@factor("high_52w") # George & Hwang (2004) +def high_52w(p): + c = p.close.ffill() + return c / c.rolling(252, min_periods=200).max() + + +@factor("fip") # Da, Gurun & Warachka (2014) frog-in-the-pan (continuous info) +def fip(p): + r = _ret(p).shift(21) + pos = (r > 0).astype(float).where(r.notna()).rolling(231, min_periods=180).mean() + neg = (r < 0).astype(float).where(r.notna()).rolling(231, min_periods=180).mean() + m = mom_12_1(p) + return -np.sign(m) * (neg - pos) # high = smooth, continuous gains + + +@factor("mom_fip") # momentum restricted to continuous-information winners +def mom_fip(p): + return mom_12_1(p).rank(axis=1, pct=True) + fip(p).rank(axis=1, pct=True) + + +@factor("rec_126") # existing house factor: distance from 126d low +def rec_126(p): + c = p.close.ffill() + return c / c.rolling(126, min_periods=126).min() - 1 + + +@factor("ma200_dist") +def ma200_dist(p): + c = p.close.ffill() + return c / c.rolling(200, min_periods=200).mean() - 1 + + +@factor("ts_trend") # fraction of 10 lookbacks (1..12m) with positive return +def ts_trend(p): + c = p.close.ffill() + looks = [21, 42, 63, 84, 105, 126, 147, 168, 189, 210, 231, 252] + return sum((c / c.shift(k) > 1).astype(float) for k in looks) / len(looks) + + +@factor("sector_mom") # Moskowitz & Grinblatt (1999) industry momentum, 6-1 +def sector_mom(p): + sec = sector_map(p) + m = mom_6_1(p) + if sec is None: + return m * np.nan + out = pd.DataFrame(np.nan, index=m.index, columns=m.columns) + mm = m.where(p.member) + for s in sec.dropna().unique(): + cols = [c for c in sec.index[sec == s] if c in m.columns] + if cols: + avg = mm[cols].mean(axis=1) + out[cols] = np.repeat(avg.values[:, None], len(cols), axis=1) + return out + + +# --------------------------------------------------------------------------- +# reversal family +# --------------------------------------------------------------------------- + +@factor("rev_5") +def rev_5(p): + c = p.close.ffill() + return -(c / c.shift(5) - 1) + + +@factor("rev_21") +def rev_21(p): + c = p.close.ffill() + return -(c / c.shift(21) - 1) + + +@factor("rev_resid_21") # residual (beta-adjusted) 1-month reversal +def rev_resid_21(p): + return -resid_returns(p).rolling(21, min_periods=15).sum() + + +@factor("rev_sector_21") # industry-adjusted reversal (Da, Liu & Schaumburg 2014 style) +def rev_sector_21(p): + return demean_by_sector(rev_21(p), p) + + +@factor("overnight_mom") # Lou, Polk & Skouras (2019): overnight returns persist +def overnight_mom(p): + o = p.open + cprev = p.close.ffill().shift(1) + on = np.log(o / cprev) + return on.rolling(252, min_periods=180).sum().shift(21) + + +@factor("intraday_rev") # intraday component of last month reverses +def intraday_rev(p): + intr = np.log(p.close / p.open) + return -intr.rolling(21, min_periods=15).sum() + + +# --------------------------------------------------------------------------- +# risk / lottery family (low-risk anomaly) +# --------------------------------------------------------------------------- + +@factor("ivol") # Ang, Hodrick, Xing & Zhang (2006): low ivol +def ivol(p): + return -resid_returns(p).rolling(63, min_periods=40).std() + + +@factor("low_vol") +def low_vol(p): + return -_ret(p).rolling(252, min_periods=180).std() + + +@factor("low_beta") # Frazzini & Pedersen (2014) BAB, long leg +def low_beta(p): + return -_rolling_beta(_ret(p), _mkt(p), 252, 180) + + +@factor("max_ret") # Bali, Cakici & Whitelaw (2011): avoid lottery stocks (MAX(1)) +def max_ret(p): + return -_ret(p).rolling(21, min_periods=15).max() + + +@factor("skew") # low skewness earns more +def skew(p): + return -_ret(p).rolling(252, min_periods=180).skew() + + +# --------------------------------------------------------------------------- +# volume / liquidity family +# --------------------------------------------------------------------------- + +def _dollar_volume(p): + if "amount" in p.extra: + return p.extra["amount"] + return p.close * p.volume + + +def _turnover(p): + if "turn" in p.extra: + return p.extra["turn"] + # US: no shares outstanding → dollar volume relative to its own 1y mean + return None + + +@factor("amihud") # Amihud (2002): illiquid earns more +def amihud(p): + dv = _dollar_volume(p).replace(0, np.nan) + return np.log((_ret(p).abs() / dv).rolling(21, min_periods=15).mean()) + + +@factor("abn_volume") # Gervais, Kaniel & Mingelgrin (2001) high-volume return premium +def abn_volume(p): + dv = _dollar_volume(p) + return np.log(dv.rolling(5, min_periods=4).mean() / dv.rolling(252, min_periods=180).mean()) + + +@factor("low_turnover") # CN: low turnover / low attention earns more +def low_turnover(p): + t = _turnover(p) + if t is None: + dv = _dollar_volume(p) + return -np.log(dv.rolling(21, min_periods=15).mean() / dv.rolling(252, min_periods=180).mean()) + return -t.rolling(21, min_periods=15).mean() + + +@factor("abn_turnover") # CN: turnover now vs 1y (negative = attention fading) +def abn_turnover(p): + t = _turnover(p) + if t is None: + return abn_volume(p) * -1 + return -np.log(t.rolling(21, min_periods=15).mean() / t.rolling(252, min_periods=180).mean()) + + +@factor("turn_vol") # CN: volatility of turnover (low is better) +def turn_vol(p): + t = _turnover(p) + if t is None: + t = _dollar_volume(p) + s = t.rolling(21, min_periods=15) + return -(s.std() / s.mean()) + + +@factor("pv_corr") # CN: price-volume correlation (negative corr is better) +def pv_corr(p): + v = _dollar_volume(p) + return -p.close.ffill().rolling(21, min_periods=15).corr(v) + + +@factor("upvol_ratio") # share of volume on up days, 21d +def upvol_ratio(p): + r = _ret(p) + v = _dollar_volume(p) + up = v.where(r > 0, 0.0).rolling(21, min_periods=15).sum() + return up / v.rolling(21, min_periods=15).sum() + + +# --------------------------------------------------------------------------- +# value / fundamentals (CN: baostock daily ratios) +# --------------------------------------------------------------------------- + +def _inv(x: pd.DataFrame) -> pd.DataFrame: + return (1.0 / x).replace([np.inf, -np.inf], np.nan) + + +@factor("ep", markets=("cn",)) +def ep(p): + return _inv(p.extra["pe"]) + + +@factor("bp", markets=("cn",)) +def bp(p): + return _inv(p.extra["pb"]) + + +@factor("sp", markets=("cn",)) +def sp(p): + return _inv(p.extra["ps"]) + + +@factor("cfp", markets=("cn",)) +def cfp(p): + return _inv(p.extra["pcf"]) + + +@factor("roe_proxy", markets=("cn",)) # ROE = PB / PE (TTM earnings / book) +def roe_proxy(p): + roe = p.extra["pb"] / p.extra["pe"] + return roe.where(p.extra["pe"] > 0) + + +@factor("size", markets=("cn",)) # small float cap earns more (within index) +def size(p): + return -np.log(p.extra["float_cap"].rolling(21, min_periods=5).mean()) + + +# --------------------------------------------------------------------------- +# composite utilities +# --------------------------------------------------------------------------- + +def zrank(df: pd.DataFrame, member: pd.DataFrame | None = None) -> pd.DataFrame: + """Cross-sectional percentile rank in [0, 1] among members (NaN elsewhere).""" + if member is not None: + df = df.where(member) + return df.rank(axis=1, pct=True) + + +def composite(p: Panel, weights: dict[str, float], cache: dict | None = None) -> pd.DataFrame: + total = None + wsum = 0.0 + for name, w in weights.items(): + f = cache[name] if cache is not None and name in cache else FACTORS[name](p) + r = zrank(f, p.member) + total = r * w if total is None else total.add(r * w, fill_value=np.nan) + wsum += abs(w) + return total / wsum + + +# --------------------------------------------------------------------------- +# literature shortlist additions +# --------------------------------------------------------------------------- + +@factor("overnight_mom_6") # 6-month overnight return, no skip +def overnight_mom_6(p): + on = np.log(p.open / p.close.ffill().shift(1)) + return on.rolling(126, min_periods=90).sum() + + +@factor("overnight_mom_1") # last month's overnight return +def overnight_mom_1(p): + on = np.log(p.open / p.close.ffill().shift(1)) + return on.rolling(21, min_periods=15).sum() + + +@factor("stmom_hiturn") # Medhat & Schmeling (2022): 1-month momentum among high-turnover +def stmom_hiturn(p): # names, skipping the last 3 days; turnover proxied by + c = p.close.ffill() # abnormal dollar volume when shares outstanding are unknown + r1 = c.shift(3) / c.shift(24) - 1 + t = _turnover(p) + if t is None: + dv = _dollar_volume(p) + t = dv.rolling(21, min_periods=15).mean() / dv.rolling(252, min_periods=180).mean() + else: + t = t.rolling(21, min_periods=15).mean() + hi = t.where(p.member).rank(axis=1, pct=True) >= 0.7 + return r1.where(hi) + + +@factor("seasonality") # Heston & Sadka (2008): same-calendar-month return, years 1-10 +def seasonality(p): + m = p.close.ffill().resample("ME").last() + mret = m.pct_change(fill_method=None) + lags = [mret.shift(12 * k) for k in range(1, 11)] + avg = sum(l.fillna(0) for l in lags) / sum(l.notna().astype(float) for l in lags).replace(0, np.nan) + # value known at month-end t predicts month t+1: shift so the forecast for + # next month (lag 11, 23, ...) is available at the start of that month + fwd = avg.shift(-1) # uses months t+1-12k, all <= t + daily = fwd.reindex(p.close.index, method="ffill") + return daily + + +@factor("trend_hzz") # Han, Zhou & Zhu (2016) trend factor +def trend_hzz(p): + c = p.close.ffill() + lags = [3, 5, 10, 20, 50, 100, 200, 400] + feats = {L: (c.rolling(L, min_periods=L).mean() / c) for L in lags} + # monthly cross-sectional regression of next-month return on MA signals + month_end = c.resample("ME").last().index + me = c.index[c.index.searchsorted(month_end, side="right") - 1] + fwd = c.shift(-21) / c - 1 + betas = [] + for d in me: + y = fwd.loc[d].where(p.member.loc[d]) + X = pd.DataFrame({L: feats[L].loc[d] for L in lags}) + ok = y.notna() & X.notna().all(axis=1) + if ok.sum() < 50: + betas.append(pd.Series(np.nan, index=lags, name=d)); continue + Xo = np.column_stack([np.ones(ok.sum()), X[ok].values]) + b = np.linalg.lstsq(Xo, y[ok].values, rcond=None)[0][1:] + betas.append(pd.Series(b, index=lags, name=d)) + B = pd.DataFrame(betas) + # the regression at month d uses returns through d+21, so it is usable at d+21; + # shift by one month, then average the last 12 known betas + Ebeta = B.shift(1).rolling(12, min_periods=6).mean() + Ed = Ebeta.reindex(c.index, method="ffill") + return sum(feats[L].mul(Ed[L], axis=0) for L in lags) + + +# --------------------------------------------------------------------------- +# A-share additions (CN literature review). Volume ratios stand in for +# turnover ratios: shares outstanding rarely change inside the window. +# --------------------------------------------------------------------------- + +def _limit_hit(p) -> pd.DataFrame: + r = _ret(p) + lim = p.extra.get("limit") + if lim is None: + return r > 0.095 + return r >= lim.reindex_like(r) - 0.005 + + +@factor("mom_ex_limitup", markets=("cn",)) # 国信: 12-1 momentum without limit-up days +def mom_ex_limitup(p): + lr = np.log1p(_ret(p)).where(~_limit_hit(p), 0.0) + return lr.shift(21).rolling(231, min_periods=180).sum() + + +@factor("amp_cut_mom", markets=("us", "cn")) # 开源: momentum on the 70% lowest-amplitude days (160d) +def amp_cut_mom(p): + amp = (p.high / p.low - 1).values + r = _ret(p).values + n, keep = 160, 0.7 + out = np.full(r.shape, np.nan) + from numpy.lib.stride_tricks import sliding_window_view as swv + A = swv(np.nan_to_num(amp, nan=np.inf), n, axis=0) # (T-n+1, N, n) + R = swv(np.nan_to_num(r, nan=0.0), n, axis=0) + for t0 in range(0, A.shape[0], 256): # chunk to bound memory + a, rr = A[t0:t0 + 256], R[t0:t0 + 256] + thr = np.quantile(np.where(np.isinf(a), np.nan, a), keep, axis=2, keepdims=True) + out[t0 + n - 1:t0 + n - 1 + a.shape[0]] = np.where(a <= thr, rr, 0.0).sum(axis=2) + valid = pd.DataFrame(amp, index=p.close.index).notna().rolling(n).sum() >= 0.8 * n + return pd.DataFrame(out, index=p.close.index, columns=p.close.columns).where(valid.values) + + +@factor("ideal_amp", markets=("us", "cn")) # 开源 理想振幅: amp on high-price days minus low-price days +def ideal_amp(p): + n, lam = 20, 0.25 + amp = (p.high / p.low - 1).values + c = p.close.ffill().values + from numpy.lib.stride_tricks import sliding_window_view as swv + A = swv(amp, n, axis=0) + Cw = swv(c, n, axis=0) + out = np.full(amp.shape, np.nan) + hi_q = np.nanquantile(Cw, 1 - lam, axis=2, keepdims=True) + lo_q = np.nanquantile(Cw, lam, axis=2, keepdims=True) + with np.errstate(invalid="ignore"): + v_hi = np.nanmean(np.where(Cw >= hi_q, A, np.nan), axis=2) + v_lo = np.nanmean(np.where(Cw <= lo_q, A, np.nan), axis=2) + out[n - 1:] = v_hi - v_lo + return -pd.DataFrame(out, index=p.close.index, columns=p.close.columns) + + +@factor("abn_vol_ratio", markets=("us", "cn")) # LSY PMO: -mean(vol,20)/mean(vol,250) +def abn_vol_ratio(p): + v = p.volume.where(p.volume > 0) + return -(v.rolling(20, min_periods=15).mean() / v.rolling(250, min_periods=180).mean()) + + +@factor("vol_cv", markets=("us", "cn")) # HQZ turnover volatility, scale-free +def vol_cv(p): + v = p.volume.where(p.volume > 0) + s = v.rolling(20, min_periods=15) + return -(s.std() / s.mean()) + + +@factor("overnight_20", markets=("us", "cn")) # 国信 overnight momentum, 20d +def overnight_20(p): + return (p.open / p.close.ffill().shift(1) - 1).rolling(20, min_periods=15).sum() + + +@factor("ivol_20", markets=("us", "cn")) +def ivol_20(p): + return -resid_returns(p).rolling(20, min_periods=15).std() diff --git a/research/lab/fetch_cn.py b/research/lab/fetch_cn.py new file mode 100644 index 0000000..e18aabb --- /dev/null +++ b/research/lab/fetch_cn.py @@ -0,0 +1,157 @@ +""" +Fetch point-in-time CSI 300 / CSI 500 membership and daily bars from baostock. + +The live CN panel (data/cn.csv) holds today's CSI 300 members only, so every +CN backtest on it carries survivorship bias. baostock serves historical index +membership, delisted stocks, suspensions (tradestatus), ST flags, turnover and +valuation ratios, which lets research backtests trade the index as it was. + +Usage: + uv run --with baostock python -m research.lab.fetch_cn [--start 2007-01-01] [--workers 4] + +Writes data/research/cn_ohlcv.pkl: + fields: {"open","high","low","close","preclose","volume","amount","turn", + "pctchg","tradestatus","isst","pe","pb","ps","pcf"} DataFrames + (dates x codes; prices back-adjusted, 后复权) + members: {"hs300": {date: [codes]}, "zz500": {date: [codes]}} + index: DataFrame of CSI 300 / CSI 500 index closes +""" + +import argparse +import os +import pickle +from datetime import date +from multiprocessing import Pool + +import pandas as pd + +OUT = os.path.join("data", "research", "cn_ohlcv.pkl") +K_FIELDS = ("date,code,open,high,low,close,preclose,volume,amount,turn," + "tradestatus,pctChg,peTTM,pbMRQ,psTTM,pcfNcfTTM,isST") +NUMERIC = ["open", "high", "low", "close", "preclose", "volume", "amount", "turn", + "pctChg", "peTTM", "pbMRQ", "psTTM", "pcfNcfTTM", "tradestatus", "isST"] +RENAME = {"pctChg": "pctchg", "peTTM": "pe", "pbMRQ": "pb", "psTTM": "ps", + "pcfNcfTTM": "pcf", "isST": "isst"} + + +def _rows(rs) -> list[list[str]]: + out = [] + while rs.error_code == "0" and rs.next(): + out.append(rs.get_row_data()) + return out + + +def _login() -> None: + import baostock as bs + bs.login() + + +def fetch_members(start: str) -> dict[str, dict[str, list[str]]]: + path = os.path.join(BAR_DIR, f"members_{start}_{date.today()}.pkl") + if os.path.exists(path): + return pd.read_pickle(path) + members = _fetch_members(start) + pd.to_pickle(members, path) + return members + + +def _fetch_members(start: str) -> dict[str, dict[str, list[str]]]: + import baostock as bs + bs.login() + months = pd.date_range(start, date.today(), freq="MS") + members: dict[str, dict[str, list[str]]] = {"hs300": {}, "zz500": {}} + for d in months: + ds = d.strftime("%Y-%m-%d") + members["hs300"][ds] = [r[1] for r in _rows(bs.query_hs300_stocks(date=ds))] + members["zz500"][ds] = [r[1] for r in _rows(bs.query_zz500_stocks(date=ds))] + bs.logout() + return members + + +BAR_DIR = os.path.join("data", "research", "cn_bars") + + +def fetch_one(args: tuple[str, str, str]) -> pd.DataFrame | None: + """Fetch one stock's bars, cached per code so an interrupted run resumes.""" + code, start, end = args + path = os.path.join(BAR_DIR, f"{code}_{end}.pkl") + if os.path.exists(path): + return pd.read_pickle(path) + df = _fetch_one(code, start, end) + pd.to_pickle(df, path) + return df + + +def _fetch_one(code: str, start: str, end: str) -> pd.DataFrame | None: + import baostock as bs + for _ in range(3): + rs = bs.query_history_k_data_plus(code, K_FIELDS, start_date=start, end_date=end, + frequency="d", adjustflag="1") + if rs.error_code == "0": + rows = _rows(rs) + if not rows: + return None + df = pd.DataFrame(rows, columns=rs.fields) + df["date"] = pd.to_datetime(df["date"]) + for c in NUMERIC: + df[c] = pd.to_numeric(df[c], errors="coerce") + return df + bs.logout() + bs.login() + print(f" failed {code}: {rs.error_msg}", flush=True) + return None + + +def fetch_index(start: str, end: str) -> pd.DataFrame: + import baostock as bs + bs.login() + out = {} + for code in ("sh.000300", "sh.000905"): + rs = bs.query_history_k_data_plus(code, "date,close", start_date=start, end_date=end, + frequency="d") + df = pd.DataFrame(_rows(rs), columns=rs.fields) + out[code] = pd.Series(pd.to_numeric(df["close"]).values, index=pd.to_datetime(df["date"])) + bs.logout() + return pd.DataFrame(out) + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--start", default="2007-01-01") + ap.add_argument("--member-start", default="2008-01-01") + ap.add_argument("--workers", type=int, default=4) + ap.add_argument("--hs300-only", action="store_true") + args = ap.parse_args() + end = date.today().strftime("%Y-%m-%d") + + os.makedirs(BAR_DIR, exist_ok=True) + print("Fetching index membership...", flush=True) + members = fetch_members(args.member_start) + hs300 = sorted({c for lst in members["hs300"].values() for c in lst}) + codes = list(hs300) + if not args.hs300_only: + codes += sorted({c for lst in members["zz500"].values() for c in lst} - set(hs300)) + print(f" {len(codes)} codes ever in the index set", flush=True) + + frames = [] + with Pool(args.workers, initializer=_login) as pool: + for i, df in enumerate(pool.imap(fetch_one, [(c, args.start, end) for c in codes])): + if df is not None: + frames.append(df) + if (i + 1) % 100 == 0: + print(f" bars {i + 1}/{len(codes)}", flush=True) + bars = pd.concat(frames, ignore_index=True) + fields = {} + for c in NUMERIC: + fields[RENAME.get(c, c)] = bars.pivot(index="date", columns="code", values=c).sort_index() + + panel = {"fields": fields, "members": members, "index": fetch_index(args.start, end)} + os.makedirs(os.path.dirname(OUT), exist_ok=True) + with open(OUT, "wb") as fh: + pickle.dump(panel, fh) + c = fields["close"] + print(f"Saved {OUT}: {c.shape[1]} codes, {c.index[0].date()} → {c.index[-1].date()}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_cn_yahoo.py b/research/lab/fetch_cn_yahoo.py new file mode 100644 index 0000000..cfadd47 --- /dev/null +++ b/research/lab/fetch_cn_yahoo.py @@ -0,0 +1,70 @@ +""" +Fetch adjusted daily OHLCV for CSI 300 / CSI 500 ever-members from Yahoo +(through the market-data proxy). Delisted names Yahoo lacks stay missing +(listed in data/research/cn_missing.txt): Eastmoney and Tencent throttle bulk +requests, and Tencent's back-adjusted series damps returns (limit-down days +show -7% instead of -10%), so neither is used. + + uv run python -m research.lab.fetch_cn_yahoo [--hs300-only] + +Writes data/research/cn_ohlcv.pkl (fields open/high/low/close/volume/preclose/ +pctchg/tradestatus + members + index + per-code "source"). +""" + +import argparse +import glob +import os +import pickle + +import pandas as pd + +from research.lab.fetch_us import _proxy_env, fetch + +OUT = os.path.join("data", "research", "cn_ohlcv.pkl") + + +def to_yahoo(code: str) -> str: + ex, num = code.split(".") + return f"{num}.{'SS' if ex == 'sh' else 'SZ'}" + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--hs300-only", action="store_true") + ap.add_argument("--start", default="2007-01-01") + args = ap.parse_args() + member_file = sorted(glob.glob(os.path.join("data", "research", "cn_bars", "members_*.pkl")))[-1] + members = pd.read_pickle(member_file) + hs300 = sorted({c for lst in members["hs300"].values() for c in lst}) + codes = list(hs300) + if not args.hs300_only: + codes += sorted({c for lst in members["zz500"].values() for c in lst} - set(hs300)) + ymap = {to_yahoo(c): c for c in codes} + _proxy_env() + raw = fetch(list(ymap) + ["000300.SS", "000905.SS"], args.start) + fields = {k: raw[k].rename(columns=ymap) for k in ("open", "high", "low", "close", "volume")} + source = {c: "yahoo" for c in fields["close"].columns if c in codes} + + idx = pd.DataFrame({"sh.000300": raw["close"].get("000300.SS"), + "sh.000905": raw["close"].get("000905.SS")}) + keep = [c for c in fields["close"].columns if c in source] + fields = {k: v[keep].sort_index() for k, v in fields.items()} + # Yahoo prints flat bars (volume 0, O=H=L=C) on suspended days: treat as no trade + close = fields["close"] + traded = (fields["volume"] > 0) & close.notna() + fields["tradestatus"] = traded.astype(float) + fields["preclose"] = close.where(traded).ffill().shift(1) + fields["pctchg"] = (close / fields["preclose"] - 1) * 100 + missing = sorted(set(codes) - set(source)) + panel = {"fields": fields, "members": members, "index": idx, "source": source, + "missing": missing} + with open(OUT, "wb") as fh: + pickle.dump(panel, fh) + print(f"Saved {OUT}: {len(keep)} codes, {len(missing)} missing, " + f"{close.index[0].date()} → {close.index[-1].date()}") + with open(os.path.join("data", "research", "cn_missing.txt"), "w") as fh: + fh.write("\n".join(missing)) + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_us.py b/research/lab/fetch_us.py new file mode 100644 index 0000000..0af114a --- /dev/null +++ b/research/lab/fetch_us.py @@ -0,0 +1,91 @@ +""" +Fetch a long OHLCV history for every ticker that was ever an S&P 500 member. + +The live panels in data/us.csv hold 10 years of close/open for tickers with +>50% coverage. Factor research needs more: volume and high/low, 20+ years of +history for an out-of-sample design split, and recently listed members (CEG, +GEV, ...) that the coverage rule drops. Delisted names that Yahoo no longer +serves stay missing; ``coverage_report`` in panel.py quantifies that gap. + +Usage: + uv run python -m research.lab.fetch_us [--start 2004-01-01] + +Writes data/research/us_ohlcv.pkl: {"open","high","low","close","volume"} +DataFrames (dates x tickers), split/dividend adjusted. +""" + +import argparse +import os +import pickle +import time + +import pandas as pd +import yfinance as yf + +import universe_history as uh + +OUT = os.path.join("data", "research", "us_ohlcv.pkl") +FIELDS = ["Open", "High", "Low", "Close", "Volume"] +# Benchmarks, cash proxies and sector ETFs used by overlays and attribution. +EXTRA = ["SPY", "RSP", "QQQ", "SHY", "BIL", "IEF", "TLT", "GLD", + "XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY"] + + +def _proxy_env() -> None: + proxy = os.environ.get("QUANT_MARKET_PROXY") + if not proxy and os.path.exists(".env"): + for line in open(".env"): + if line.startswith("QUANT_MARKET_PROXY="): + proxy = line.split("=", 1)[1].strip().strip("'\"") + if proxy: + for k in ("http_proxy", "https_proxy", "HTTP_PROXY", "HTTPS_PROXY"): + os.environ[k] = proxy + + +def fetch(tickers: list[str], start: str, chunk: int = 60) -> dict[str, pd.DataFrame]: + parts: dict[str, list[pd.DataFrame]] = {f: [] for f in FIELDS} + for i in range(0, len(tickers), chunk): + batch = tickers[i:i + chunk] + for attempt in range(4): + try: + raw = yf.download(batch, start=start, auto_adjust=True, progress=False, + group_by="column", threads=True) + break + except Exception as exc: # noqa: BLE001 — retry transient Yahoo errors + print(f" batch {i}: {exc}; retry {attempt + 1}") + time.sleep(10 * (attempt + 1)) + else: + continue + for f in FIELDS: + if f in raw.columns.get_level_values(0): + parts[f].append(raw[f]) + print(f" fetched {min(i + chunk, len(tickers))}/{len(tickers)}", flush=True) + time.sleep(2) + out = {} + for f in FIELDS: + if not parts[f]: + raise RuntimeError(f"no {f} data downloaded") + df = pd.concat(parts[f], axis=1) + df = df.loc[:, ~df.columns.duplicated()].sort_index() + out[f.lower()] = df.dropna(axis=1, how="all") + return out + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--start", default="2004-01-01") + args = ap.parse_args() + _proxy_env() + intervals = uh.load_sp500_history() + tickers = sorted(set(uh.all_tickers_ever(intervals)) | set(EXTRA)) + print(f"Fetching {len(tickers)} tickers from {args.start}") + panel = fetch(tickers, args.start) + os.makedirs(os.path.dirname(OUT), exist_ok=True) + with open(OUT, "wb") as fh: + pickle.dump(panel, fh) + c = panel["close"] + print(f"Saved {OUT}: {c.shape[1]} tickers, {c.index[0].date()} → {c.index[-1].date()}") + + +if __name__ == "__main__": + main() diff --git a/research/lab/fetch_us_sectors.py b/research/lab/fetch_us_sectors.py new file mode 100644 index 0000000..ba2e189 --- /dev/null +++ b/research/lab/fetch_us_sectors.py @@ -0,0 +1,43 @@ +"""Fetch GICS-like sectors from Yahoo for S&P 500 ever-members missing from +data/us_sectors.csv. Writes data/research/us_sectors_extra.csv. + + uv run python -m research.lab.fetch_us_sectors +""" + +import os +import pickle + +import pandas as pd +import yfinance as yf + +from research.lab.fetch_us import _proxy_env + +# Yahoo sector names → GICS sector names used in data/us_sectors.csv +YAHOO_TO_GICS = { + "Technology": "Information Technology", "Healthcare": "Health Care", + "Financial Services": "Financials", "Consumer Cyclical": "Consumer Discretionary", + "Consumer Defensive": "Consumer Staples", "Basic Materials": "Materials", + "Communication Services": "Communication Services", "Industrials": "Industrials", + "Energy": "Energy", "Utilities": "Utilities", "Real Estate": "Real Estate", +} + + +def main() -> None: + _proxy_env() + with open(os.path.join("data", "research", "us_ohlcv.pkl"), "rb") as fh: + tickers = list(pickle.load(fh)["close"].columns) + known = set(pd.read_csv(os.path.join("data", "us_sectors.csv"))["ticker"]) + rows = [] + for t in sorted(set(tickers) - known): + try: + s = yf.Ticker(t).info.get("sector") + except Exception: # noqa: BLE001 + s = None + rows.append({"ticker": t, "sector": YAHOO_TO_GICS.get(s, s)}) + out = pd.DataFrame(rows) + out.to_csv(os.path.join("data", "research", "us_sectors_extra.csv"), index=False) + print(out["sector"].value_counts(dropna=False)) + + +if __name__ == "__main__": + main() diff --git a/research/lab/ml.py b/research/lab/ml.py new file mode 100644 index 0000000..329b851 --- /dev/null +++ b/research/lab/ml.py @@ -0,0 +1,95 @@ +""" +Walk-forward signal combination. + +Features are cross-sectional percentile ranks of every factor among PIT +members on each date; the target is the cross-sectional rank of the forward +open-to-open return over ``horizon`` days. Each calendar year is predicted by +a model trained only on samples whose labels ended before that year starts +(a ``horizon``-day purge), so every score is out-of-sample. + + from research.lab import ml + feats = ml.feature_ranks(panel, names) + score = ml.walk_forward(panel, feats, model="lgbm", horizon=20) +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from research.lab.factors import FACTORS +from research.lab.report import forward_returns + + +def feature_ranks(p, names: list[str]) -> dict[str, pd.DataFrame]: + out = {} + for n in names: + f = FACTORS[n](p).reindex(index=p.dates, columns=p.close.columns) + out[n] = f.where(p.member).rank(axis=1, pct=True) + return out + + +def _stack(feats: dict[str, pd.DataFrame], rows: pd.DatetimeIndex) -> pd.DataFrame: + cols = {n: f.loc[rows].stack(future_stack=True) for n, f in feats.items()} + return pd.DataFrame(cols) + + +def walk_forward(p, feats: dict[str, pd.DataFrame], model: str = "lgbm", horizon: int = 20, + sample_every: int = 5, first_test_year: int = 2010, min_train_years: int = 4, + params: dict | None = None, seed: int = 0) -> pd.DataFrame: + fwd = forward_returns(p, horizon).where(p.member) + target = fwd.rank(axis=1, pct=True) + dates = p.dates + sample_rows = dates[::sample_every] + X_all = _stack(feats, sample_rows) + y_all = target.loc[sample_rows].stack(future_stack=True).reindex(X_all.index) + ok = y_all.notna() & (X_all.notna().sum(axis=1) >= len(feats) // 2) + X_all, y_all = X_all[ok], y_all[ok] + sample_dates = X_all.index.get_level_values(0) + + score = pd.DataFrame(np.nan, index=dates, columns=p.close.columns) + years = sorted({d.year for d in dates if d.year >= first_test_year}) + for yr in years: + test_start = pd.Timestamp(f"{yr}-01-01") + pos = dates.searchsorted(test_start) + if pos >= len(dates): + break + # purge: the label of a sample at d uses opens through d + horizon + 1 + cutoff = dates[max(pos - horizon - 2, 0)] + tr = sample_dates <= cutoff + if (sample_dates[tr].max() - sample_dates[tr].min()).days < 365 * min_train_years: + continue + Xtr, ytr = X_all[tr], y_all[tr] + test_rows = dates[(dates >= test_start) & (dates < pd.Timestamp(f"{yr + 1}-01-01"))] + Xte = _stack(feats, test_rows) + pred = _fit_predict(model, Xtr, ytr, Xte, params, seed) + s = pd.Series(pred, index=Xte.index).unstack() + score.loc[s.index, s.columns] = s.values + return score + + +def _fit_predict(model, Xtr, ytr, Xte, params, seed): + if model == "lgbm": + import lightgbm as lgb + prm = dict(objective="regression", learning_rate=0.03, num_leaves=31, + min_data_in_leaf=1000, bagging_fraction=0.7, bagging_freq=1, + feature_fraction=0.7, lambda_l2=1.0, seed=seed, verbose=-1, num_threads=8) + prm.update(params or {}) + rounds = prm.pop("n_estimators", 300) + booster = lgb.train(prm, lgb.Dataset(Xtr.values, ytr.values), num_boost_round=rounds) + return booster.predict(Xte.values) + if model == "ridge": + Xf = Xtr.fillna(0.5).values - 0.5 + lam = (params or {}).get("lam", 10.0) * len(Xf) + beta = np.linalg.solve(Xf.T @ Xf + lam * np.eye(Xf.shape[1]), Xf.T @ (ytr.values - 0.5)) + return (Xte.fillna(0.5).values - 0.5) @ beta + if model == "ic_weight": + # sign-and-size each feature by its training-period mean rank IC + d = Xtr.index.get_level_values(0) + ics = {} + for c in Xtr.columns: + df = pd.DataFrame({"x": Xtr[c], "y": ytr, "d": d}).dropna() + ics[c] = df.groupby("d").apply(lambda g: g["x"].corr(g["y"])).mean() + w = pd.Series(ics).fillna(0.0) + return (Xte.fillna(0.5) - 0.5).values @ w.values + raise ValueError(model) diff --git a/research/lab/overlays.py b/research/lab/overlays.py new file mode 100644 index 0000000..669b316 --- /dev/null +++ b/research/lab/overlays.py @@ -0,0 +1,55 @@ +""" +Exposure overlays: series in [0, 1] read on signal dates by engine.simulate +(``exposure=``). The un-invested part earns the panel's cash return. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from research.lab.engine import simulate + + +def ew_index(p) -> pd.Series: + """Equal-weight index of PIT members (daily rebalanced).""" + r = p.close.pct_change(fill_method=None).where(p.member.shift(1, fill_value=False)) + return (1 + r.mean(axis=1).fillna(0)).cumprod() + + +def trend(level: pd.Series, ma: int = 200, floor: float = 0.0) -> pd.Series: + """1 when the index is above its moving average, else ``floor``.""" + above = level > level.rolling(ma, min_periods=ma).mean() + return above.astype(float).where(above, floor) + + +def breadth(p, ma: int = 200, lo: float = 0.4, hi: float = 0.6) -> pd.Series: + """Share of members above their own MA mapped linearly from lo→0 to hi→1.""" + c = p.close.ffill() + above = (c > c.rolling(ma, min_periods=ma).mean()).where(p.member) + share = above.mean(axis=1) + return ((share - lo) / (hi - lo)).clip(0, 1) + + +def drawdown_brake(level: pd.Series, dd: float = 0.15, ma: int = 50) -> pd.Series: + """Cut to zero when the index is > dd below its 1y high AND below its MA.""" + high = level.rolling(252, min_periods=60).max() + bad = (level / high - 1 < -dd) & (level < level.rolling(ma, min_periods=ma).mean()) + return (~bad).astype(float) + + +def vol_target(p, score, target: float, lookback: int = 60, cap: float = 1.0, + base_kwargs: dict | None = None, exposure: pd.Series | None = None) -> pd.Series: + """Scale exposure by target / realized vol of the unscaled strategy. + + Two passes: run the strategy fully invested, measure its trailing + daily-return vol, then size the next rebalance; any given ``exposure`` + (e.g. a trend filter) multiplies the result. + """ + kw = dict(base_kwargs or {}) + base = simulate(p, score, **kw) + rv = base.returns.rolling(lookback, min_periods=lookback // 2).std() * np.sqrt(252) + scale = (target / rv).clip(upper=cap).fillna(cap) + if exposure is not None: + scale = scale * exposure.reindex(scale.index).ffill().fillna(1.0) + return scale diff --git a/research/lab/panel.py b/research/lab/panel.py new file mode 100644 index 0000000..5754378 --- /dev/null +++ b/research/lab/panel.py @@ -0,0 +1,196 @@ +""" +Build research Panels (see engine.Panel) from the fetched raw data. + + from research.lab import panel + us = panel.load_us() # S&P 500 PIT, 2004 → + cn = panel.load_cn("hs300") # CSI 300 PIT, 2007 → ("zz800" = CSI 300 + 500) + +Panels are cached in data/research/*.panel.pkl; delete them after refetching. +""" + +from __future__ import annotations + +import json +import os +import pickle + +import numpy as np +import pandas as pd + +from research.lab.engine import Costs, Panel + +RAW_DIR = os.path.join("data", "research") +US_ETFS = {"SPY", "RSP", "QQQ", "SHY", "BIL", "IEF", "TLT", "GLD", + "XLB", "XLE", "XLF", "XLI", "XLK", "XLP", "XLU", "XLV", "XLY"} + + +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) + return p + + +# --------------------------------------------------------------------------- +# US +# --------------------------------------------------------------------------- + +def _us_membership(dates: pd.DatetimeIndex, tickers: list[str]) -> pd.DataFrame: + with open(os.path.join("data", "sp500_history.json")) as fh: + intervals = json.load(fh)["intervals"] + mask = pd.DataFrame(False, index=dates, columns=tickers) + for t in tickers: + for start, end in intervals.get(t, []): + s = pd.Timestamp(start) if start else dates[0] + e = pd.Timestamp(end) if end else dates[-1] + pd.Timedelta(days=1) + mask.loc[(dates >= s) & (dates < e), t] = True + return mask + + +def _build_us() -> Panel: + with open(os.path.join(RAW_DIR, "us_ohlcv.pkl"), "rb") as fh: + raw = pickle.load(fh) + close = raw["close"] + # keep rows where most listed tickers printed (drops stray holiday rows) + alive = close.notna().sum(axis=1) + close = close[alive >= 0.5 * alive.rolling(20, min_periods=1).max()] + dates = close.index + etf = {k: raw[k].reindex(dates)[[c for c in raw[k].columns if c in US_ETFS]] + for k in ("close",)}["close"] + stocks = sorted(c for c in close.columns if c not in US_ETFS) + f = {k: raw[k].reindex(index=dates, columns=stocks) for k in ("open", "high", "low", "close", "volume")} + # Yahoo sometimes prints zero/negative adjusted opens for old bars + for k in ("open", "high", "low"): + f[k] = f[k].where(f[k] > 0) + has_px = f["close"].notna() + member = _us_membership(dates, stocks) & has_px + cash_px = etf["BIL"].combine_first(etf["SHY"]) if "BIL" in etf else etf["SHY"] + 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) + del cash_px + return Panel( + market="us", 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=etf["SPY"], + costs=Costs(buy=0.001, sell=0.001, fixed_fee=0.0, capital=10_000.0), + extra={"etf": etf, "sector": _us_sectors(stocks)}, + ) + + +def _us_sectors(stocks: list[str]) -> pd.Series: + """GICS sector per ticker; current members from data/us_sectors.csv, the + rest from data/research/us_sectors_extra.csv when it has been fetched.""" + sec = pd.read_csv(os.path.join("data", "us_sectors.csv")).set_index("ticker")["GICS Sector"] + extra = os.path.join(RAW_DIR, "us_sectors_extra.csv") + if os.path.exists(extra): + sec = sec.combine_first(pd.read_csv(extra).set_index("ticker")["sector"]) + return sec.reindex(stocks) + + +def load_us() -> Panel: + return _cached("us", _build_us) + + +# --------------------------------------------------------------------------- +# CN +# --------------------------------------------------------------------------- + +def _limit_pct(codes: list[str], dates: pd.DatetimeIndex, isst: pd.DataFrame) -> pd.DataFrame: + """Daily price-limit fraction per stock (10% main board, 20% ChiNext after + 2020-08-24 and STAR, 5% ST).""" + lim = pd.DataFrame(0.10, index=dates, columns=codes) + star = [c for c in codes if c.split(".")[1].startswith("688")] + chinext = [c for c in codes if c.split(".")[1].startswith("30")] + lim[star] = 0.20 + lim.loc[dates >= "2020-08-24", chinext] = 0.20 + st = isst.reindex(index=dates, columns=codes).fillna(0).astype(bool) + lim = lim.mask(st & (lim < 0.2), 0.05) + return lim + + +def _cn_membership(members: dict[str, list[str]], dates: pd.DatetimeIndex, + codes: list[str]) -> pd.DataFrame: + snap = sorted(members.items()) + mask = pd.DataFrame(False, index=dates, columns=codes) + col = {c: i for i, c in enumerate(codes)} + arr = mask.values + for i, (d, lst) in enumerate(snap): + s = pd.Timestamp(d) + e = pd.Timestamp(snap[i + 1][0]) if i + 1 < len(snap) else dates[-1] + pd.Timedelta(days=1) + rows = (dates >= s) & (dates < e) + cols = [col[c] for c in lst if c in col] + arr[np.ix_(rows, cols)] = True + return pd.DataFrame(arr, index=dates, columns=codes) + + +def _build_cn(index_set: str) -> Panel: + with open(os.path.join(RAW_DIR, "cn_ohlcv.pkl"), "rb") as fh: + raw = pickle.load(fh) + F = raw["fields"] + # Codes an earlier fetcher filled from Tencent's 后复权 series are dropped: + # it damps returns for long dividend histories (limit-down shows -7%). + bad = {c for c, src in raw.get("source", {}).items() if src == "tencent"} + close = F["close"] + dates = close.index + codes = [c for c in close.columns if c not in bad] + f = {k: v.reindex(index=dates, columns=codes) for k, v in F.items()} + # Yahoo bars lack ST flags and valuation ratios; index members + # are rarely ST (the index drops them at the next review). + if "isst" not in f: + f["isst"] = pd.DataFrame(0.0, index=dates, columns=codes) + traded = (f["tradestatus"] == 1) & (f["volume"] > 0) & f["close"].notna() + if index_set == "hs300": + member = _cn_membership(raw["members"]["hs300"], dates, codes) + elif index_set == "zz500": + member = _cn_membership(raw["members"]["zz500"], dates, codes) + else: # zz800 + member = (_cn_membership(raw["members"]["hs300"], dates, codes) + | _cn_membership(raw["members"]["zz500"], dates, codes)) + st = f["isst"].fillna(0).astype(bool) + member = member & traded.rolling(5, min_periods=1).max().astype(bool) & ~st + lim = _limit_pct(codes, dates, f["isst"]) + open_ret = f["open"] / f["preclose"] - 1 + can_buy = traded & ~(open_ret >= lim - 0.002) + can_sell = traded & ~(open_ret <= -lim + 0.002) + # halted days carry the last price, so returns across a halt are zero + for k in ("open", "high", "low", "close"): + f[k] = f[k].where(traded) + stamp = pd.Series(0.001, index=dates) + stamp[dates >= "2023-08-28"] = 0.0005 + idx = raw["index"].reindex(dates).ffill() + return Panel( + market=f"cn_{index_set}", open=f["open"], high=f["high"], low=f["low"], close=f["close"], + volume=f["volume"].where(traded), member=member, can_buy=can_buy, can_sell=can_sell, + cash_ret=pd.Series(0.02 / 252, index=dates), benchmark=idx["sh.000300"], + costs=Costs(buy=0.0008, sell=0.0008, sell_tax=stamp, fixed_fee=0.0, capital=70_000.0), + extra=_cn_extra(f, traded, lim, idx), + ) + + +def _cn_extra(f: dict, traded: pd.DataFrame, lim: pd.DataFrame, idx: pd.DataFrame) -> dict: + extra = {"limit": lim, "index": idx, "source": "baostock" if "pe" in f else "yahoo"} + if "turn" in f: + extra["turn"] = f["turn"].where(traded) + if "amount" in f: + extra["amount"] = f["amount"].where(traded) + if "turn" in f: + extra["float_cap"] = (f["amount"] / (f["turn"] / 100)).where(traded & (f["turn"] > 0)) + for k in ("pe", "pb", "ps", "pcf"): + if k in f: + extra[k] = f[k] + return extra + + +def load_cn(index_set: str = "hs300") -> Panel: + return _cached(f"cn_{index_set}", lambda: _build_cn(index_set)) + + +def coverage_report(p: Panel) -> pd.DataFrame: + """Per year: PIT members with price data vs expected index size.""" + m = p.member.resample("YE").last() + return pd.DataFrame({"members_with_data": m.sum(axis=1)}) diff --git a/research/lab/report.py b/research/lab/report.py new file mode 100644 index 0000000..e00d4f3 --- /dev/null +++ b/research/lab/report.py @@ -0,0 +1,118 @@ +"""Performance statistics, factor IC analysis and overfitting controls.""" + +from __future__ import annotations + +import math + +import numpy as np +import pandas as pd +from scipy import stats as sps + +import metrics + + +def stats(equity: pd.Series, turnover: pd.Series | None = None) -> dict: + eq = equity.dropna() + r = eq.pct_change().dropna() + out = { + "cagr": metrics.annualized_return(eq), + "vol": metrics.annualized_volatility(r), + "sharpe": metrics.sharpe_ratio(r), + "sortino": metrics.sortino_ratio(r), + "mdd": metrics.max_drawdown(eq), + } + out["calmar"] = out["cagr"] / abs(out["mdd"]) if out["mdd"] < 0 else np.inf + if turnover is not None: + t = turnover.loc[eq.index[0]:eq.index[-1]] + years = max((eq.index[-1] - eq.index[0]).days / 365.25, 1e-9) + out["turnover"] = t.sum() / years # one-way, x per year + return out + + +def fmt(s: dict) -> str: + base = (f"CAGR {s['cagr']*100:6.1f}% vol {s['vol']*100:5.1f}% Sharpe {s['sharpe']:5.2f} " + f"Sortino {s['sortino']:5.2f} MDD {s['mdd']*100:6.1f}% Calmar {s['calmar']:5.2f}") + if "turnover" in s: + base += f" turn {s['turnover']:5.1f}x" + return base + + +def yearly(equity: pd.Series) -> pd.Series: + e = equity.dropna() + ye = e.resample("YE").last() + first = e.iloc[0] + return ye.pct_change().fillna(ye.iloc[0] / first - 1) + + +# --------------------------------------------------------------------------- +# Factor analytics +# --------------------------------------------------------------------------- + +def forward_returns(panel, horizon: int = 5) -> pd.DataFrame: + """Return from the open of t+1 to the open of t+1+horizon — the return a + signal formed at the close of t can actually earn.""" + o = panel.open.where(panel.open > 0).ffill() + return o.shift(-(1 + horizon)) / o.shift(-1) - 1 + + +def ic_series(score: pd.DataFrame, fwd: pd.DataFrame, member: pd.DataFrame, + every: int = 5) -> pd.Series: + """Cross-sectional Spearman rank IC on every ``every``-th date.""" + idx = score.index[::every] + out = {} + for d in idx: + s = score.loc[d].where(member.loc[d]) + f = fwd.loc[d] + ok = s.notna() & f.notna() + if ok.sum() < 30: + continue + out[d] = sps.spearmanr(s[ok], f[ok]).statistic + return pd.Series(out) + + +def top_minus_universe(score: pd.DataFrame, fwd: pd.DataFrame, member: pd.DataFrame, + top_n: int = 20, every: int = 5) -> pd.Series: + """Mean forward return of the top-N names minus the member average.""" + idx = score.index[::every] + out = {} + for d in idx: + s = score.loc[d].where(member.loc[d]).dropna() + f = fwd.loc[d] + s = s[f.reindex(s.index).notna()] + if len(s) < top_n * 3: + continue + top = s.nlargest(top_n).index + out[d] = f[top].mean() - f[s.index].mean() + return pd.Series(out) + + +def ic_summary(ic: pd.Series, periods_per_year: float = 52) -> dict: + ic = ic.dropna() + if len(ic) < 5: + return {"ic": np.nan, "icir": np.nan, "t": np.nan, "hit": np.nan, "n": len(ic)} + return {"ic": ic.mean(), "icir": ic.mean() / ic.std() * math.sqrt(periods_per_year), + "t": ic.mean() / ic.std() * math.sqrt(len(ic)), "hit": (ic > 0).mean(), "n": len(ic)} + + +# --------------------------------------------------------------------------- +# Overfitting controls (Bailey & Lopez de Prado 2014) +# --------------------------------------------------------------------------- + +def deflated_sharpe(returns: pd.Series, n_trials: int, sr_trials_var: float | None = None, + periods: int = 252) -> dict: + """Probability that the true Sharpe exceeds the best Sharpe expected from + ``n_trials`` unskilled strategies. Uses per-period Sharpe internally.""" + r = returns.dropna() + T = len(r) + sr = r.mean() / r.std() + g3 = sps.skew(r) + g4 = sps.kurtosis(r, fisher=False) + if sr_trials_var is None: + sr_trials_var = 1.0 / T # variance of SR estimates under H0 + emc = 0.5772156649 + z1 = sps.norm.ppf(1 - 1.0 / n_trials) if n_trials > 1 else 0.0 + z2 = sps.norm.ppf(1 - 1.0 / (n_trials * math.e)) if n_trials > 1 else 0.0 + sr0 = math.sqrt(sr_trials_var) * ((1 - emc) * z1 + emc * z2) + denom = math.sqrt(max(1 - g3 * sr + (g4 - 1) / 4 * sr ** 2, 1e-12)) + dsr = sps.norm.cdf((sr - sr0) * math.sqrt(T - 1) / denom) + return {"sr_annual": sr * math.sqrt(periods), "sr0_annual": sr0 * math.sqrt(periods), "dsr": dsr} diff --git a/research/lab/run_ml.py b/research/lab/run_ml.py new file mode 100644 index 0000000..2e1686c --- /dev/null +++ b/research/lab/run_ml.py @@ -0,0 +1,52 @@ +""" +Walk-forward ML combination experiment. + + PYTHONPATH=. uv run --with lightgbm python -m research.lab.run_ml --market us +""" + +import argparse +import os +import pickle +import time +import warnings + + +from research.lab import ml +from research.lab.factor_scan import load, phase_avg +from research.lab.factors import FACTORS +from research.lab.report import fmt + +warnings.filterwarnings("ignore") +OUT = os.path.join("data", "research", "results") + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + ap.add_argument("--models", nargs="*", default=["lgbm", "ridge", "ic_weight"]) + ap.add_argument("--horizons", nargs="*", type=int, default=[5, 20]) + ap.add_argument("--first-test-year", type=int, default=2010) + args = ap.parse_args() + mk = "us" if args.market == "us" else "cn" + p, periods = load(args.market) + names = [n for n, f in FACTORS.items() if mk in f.markets] + t = time.time() + feats = ml.feature_ranks(p, names) + print(f"{len(names)} features in {time.time() - t:.0f}s", flush=True) + evalp = {"design": (f"{args.first_test_year}-01-01", periods["design"][1]), "holdout": periods["holdout"]} + os.makedirs(OUT, exist_ok=True) + for model in args.models: + for h in args.horizons: + t = time.time() + score = ml.walk_forward(p, feats, model=model, horizon=h, first_test_year=args.first_test_year) + with open(os.path.join(OUT, f"ml_{args.market}_{model}_h{h}.pkl"), "wb") as fh: + pickle.dump(score, fh) + for top_n in (10, 20): + parts = [] + for per, (a, b) in evalp.items(): + parts.append(f"{per}: {fmt(phase_avg(p, score, a, b, top_n, 5))}") + print(f"{model:9s} h{h:<3d} top{top_n:<3d}", " | ".join(parts), f"[{time.time() - t:.0f}s]", flush=True) + + +if __name__ == "__main__": + main() diff --git a/research/lab/stepwise.py b/research/lab/stepwise.py new file mode 100644 index 0000000..b06898c --- /dev/null +++ b/research/lab/stepwise.py @@ -0,0 +1,116 @@ +""" +Forward stepwise selection of equal-weight rank composites, scored on the +DESIGN period only; the holdout is evaluated once for the final picks. + + PYTHONPATH=. uv run python -m research.lab.stepwise --market us --max-k 4 + +Every evaluated composite is appended to the trial ledger +(data/research/results/trials_{market}.csv) so deflated Sharpe can count them. +""" + +import argparse +import os +import time +import warnings + +import numpy as np +import pandas as pd + +from research.lab.factor_scan import RESULTS, load, phase_avg +from research.lab.factors import FACTORS, zrank + +warnings.filterwarnings("ignore") + +# Four ~5-year blocks; "maximin" scores a composite by its worst block Sharpe. +SUBPERIODS = { + "us": [("2006-09-27", "2011-09-26"), ("2011-09-27", "2016-09-26"), + ("2016-09-27", "2021-09-26"), ("2021-09-27", "2026-09-24")], + "cn": [("2010-01-01", "2013-06-30"), ("2013-07-01", "2016-09-26"), + ("2016-09-27", "2021-09-26"), ("2021-09-27", "2026-09-24")], +} + + +def main() -> None: + ap = argparse.ArgumentParser() + ap.add_argument("--market", default="us") + ap.add_argument("--max-k", type=int, default=4) + ap.add_argument("--top-n", type=int, default=20) + ap.add_argument("--every", type=int, default=5) + ap.add_argument("--beam", type=int, default=3) + ap.add_argument("--objective", default="sharpe") + ap.add_argument("--exclude", nargs="*", default=[]) + ap.add_argument("--buffer", type=float, default=0.0) + args = ap.parse_args() + mk = "us" if args.market == "us" else "cn" + p, periods = load(args.market) + names = [n for n, f in FACTORS.items() if mk in f.markets and n not in args.exclude] + t = time.time() + ranks = {} + for n in names: + try: + ranks[n] = zrank(FACTORS[n](p), p.member) + except Exception as exc: # noqa: BLE001 + print(f"skip {n}: {exc}") + names = list(ranks) + print(f"{len(names)} factor ranks in {time.time() - t:.0f}s", flush=True) + a, b = periods["design"] + ledger = [] + subs = SUBPERIODS[mk] + + def evaluate(combo: tuple[str, ...]) -> dict: + s = sum(ranks[n] for n in combo) / len(combo) + if args.objective == "maximin": + per = [phase_avg(p, s, sa, sb, args.top_n, args.every, buffer=args.buffer) for sa, sb in subs] + st = {k: float(np.mean([x[k] for x in per])) for k in per[0]} + st["maximin"] = min(x["sharpe"] for x in per) + st["sub_sharpes"] = " ".join(f"{x['sharpe']:.2f}" for x in per) + st["sub_cagrs"] = " ".join(f"{x['cagr']*100:.0f}" for x in per) + else: + st = phase_avg(p, s, a, b, args.top_n, args.every, buffer=args.buffer) + ledger.append({"combo": "+".join(combo), "objective": args.objective, **st}) + return st + + beam = [()] + best_overall = [] + for k in range(1, args.max_k + 1): + cands = {} + for base in beam: + for n in names: + if n in base: + continue + combo = tuple(sorted(base + (n,))) + if combo in cands: + continue + cands[combo] = evaluate(combo) + ranked = sorted(cands.items(), key=lambda kv: -kv[1][args.objective]) + beam = [c for c, _ in ranked[:args.beam]] + print(f"--- k={k}: {len(cands)} combos", flush=True) + for c, st in ranked[:6]: + extra = (f" min {st['maximin']:.2f} subSh [{st['sub_sharpes']}] subCAGR [{st['sub_cagrs']}]" + if "maximin" in st else "") + print(f" {'+'.join(c):55s} sh {st['sharpe']:.2f} cagr {st['cagr']*100:5.1f}% " + f"mdd {st['mdd']*100:5.1f}% turn {st['turnover']:.0f}x{extra}", flush=True) + best_overall += ranked[:args.beam] + + if args.objective == "maximin": + _save(ledger, args.market) + return + # holdout, once, for the beam winners of each k + print("=== holdout of design winners ===") + ha, hb = periods["holdout"] + for c, st in best_overall: + s = sum(ranks[n] for n in c) / len(c) + hs = phase_avg(p, s, ha, hb, args.top_n, args.every) + print(f" {'+'.join(c):55s} D sh {st['sharpe']:.2f} cagr {st['cagr']*100:5.1f}% | " + f"H sh {hs['sharpe']:.2f} cagr {hs['cagr']*100:5.1f}% mdd {hs['mdd']*100:5.1f}%", flush=True) + _save(ledger, args.market) + + +def _save(ledger: list[dict], market: str) -> None: + os.makedirs(RESULTS, exist_ok=True) + path = os.path.join(RESULTS, f"trials_{market}.csv") + pd.DataFrame(ledger).to_csv(path, mode="a", header=not os.path.exists(path), index=False) + + +if __name__ == "__main__": + main()