claude_dev
404330f29a
旧实现 λI 直接加在被 WLS 权重放大的原始正规阵上: 真实 √市值 量级 (w≈6e4, diag(A)~1e8)下相对收缩 ~1e-9, ridge 解≡OLS 解, 对照通道在 唯一触发场景(VIF>10)下必然 ridge_flip 恒空静默放行——spec §4.7-3 「标准化设计阵固定档」名不副实。 修法: 惩罚随各列方差缩放 A+λ·diag(A)(数学上等价 D^½(R+λI)D^½ 的 标准化解)——权重整体缩放 c 时解不变(尺度不变性), 等权/√市值双通道 同效; 系数统一收缩 1/(1+λ), 截距不罚。当前 quant12 vif≈1.4 不进 触发区, 已产归因件不受影响(通道仅 VIF>10 才跑)。 测试 427→428: ①尺度不变+真实量级生效双性质钉死(旧 no-op 必红); ②flip 测试从 isinstance-only 强化为多路共线构造下非空断言 (src4=三源等权组合: 两两 ρ≈0.58 不并族但对前三源 R²≈1, VIF 全线 34~143; OLS 把 premia 为正的 src1 窗累计拧负, ridge 翻正捕获)。 因子档 §7.5 ridge 语义同步修订(档随码走)。 Co-Authored-By: Claude Code <noreply@anthropic.com>
567 lines
24 KiB
Python
567 lines
24 KiB
Python
"""因子贡献分解(P3 第三刀, spec §4.7 五条, 2026-09-25 拍板).
|
||
|
||
日度 Fama-MacBeth: 个股日收益 ~ 标准化因子暴露 + 控制集(市场截距 +
|
||
log 市值 + 行业哑变量)逐日截面 WLS 回归 → 因子日收益序列 f_k(t);
|
||
组合贡献 = 市场 + Σ_k 暴露·f_k + 行业 + 特质, 恒等式可加和验证
|
||
(特质单列——占比高=解释力弱, 结论降权)。
|
||
|
||
共线 = 决议 K 0.7 族压缩为默认回归元(quant12 纯量价 12 源零 size/行业,
|
||
不加控制集则小票潮汐系统性泼给 vol/量能代表源); VIF 披露; 代表后仍
|
||
VIF>10 触发 ridge 对照(固定 λ=1.0 标准化阵, 月度结论翻转才报警)。
|
||
不借 PCA/lasso。加权 = sqrt(市值) WLS 默认 + 等权敏感性固定双跑。
|
||
|
||
持仓 = 截面确定性函数重建(模型组合: TopN 等权硬切; ST/涨跌停/次新等
|
||
执行级过滤不重建——模型→实盘差距归偏差日报半边, 归因两半天然对齐)。
|
||
跑位 = CLI 独立, 与月度批评首班零耦合。
|
||
"""
|
||
from __future__ import annotations
|
||
|
||
import sqlite3
|
||
|
||
import numpy as np
|
||
import pandas as pd
|
||
|
||
|
||
def load_sections(sections_dir: str, names: list[str]) -> dict[str, pd.DataFrame]:
|
||
"""装载截面宽表(parquet: index=date × columns=vt_symbol). 缺文件 fail-loud."""
|
||
import os
|
||
out: dict[str, pd.DataFrame] = {}
|
||
for n in names:
|
||
p = os.path.join(sections_dir, f"{n}.parquet")
|
||
if not os.path.exists(p):
|
||
raise FileNotFoundError(f"截面缺失: {p}(先跑源级导出批)")
|
||
df = pd.read_parquet(p)
|
||
df.index = pd.to_datetime(df.index)
|
||
out[n] = df
|
||
return out
|
||
|
||
|
||
def load_closes(vnpy_db: str, symbols: list[str],
|
||
start: str, end: str) -> pd.DataFrame:
|
||
"""日线收盘 pivot 宽表(index='YYYY-MM-DD' str). 逐符号索引查, 秒级."""
|
||
from sanguo_data.datareader import guess_exchange
|
||
rows: list[tuple[str, str, float]] = []
|
||
conn = sqlite3.connect(vnpy_db, timeout=60)
|
||
conn.execute("PRAGMA busy_timeout=30000")
|
||
try:
|
||
for vt in symbols:
|
||
sym = vt.split(".", 1)[0]
|
||
ex = guess_exchange(sym).value
|
||
cur = conn.execute(
|
||
"SELECT datetime, close_price FROM dbbardata "
|
||
"WHERE symbol=? AND exchange=? AND interval='d' "
|
||
"AND datetime>=? AND datetime<=? ORDER BY datetime",
|
||
(sym, ex, start, f"{end} 23:59:59"))
|
||
for dt, close in cur.fetchall():
|
||
rows.append((str(dt)[:10], vt, float(close)))
|
||
finally:
|
||
conn.close()
|
||
df = pd.DataFrame(rows, columns=["date", "vt_symbol", "close"])
|
||
if df.empty:
|
||
return pd.DataFrame()
|
||
return df.pivot(index="date", columns="vt_symbol", values="close").sort_index()
|
||
|
||
|
||
def forward_returns(closes: pd.DataFrame) -> pd.DataFrame:
|
||
"""实现日收益: index t 的值 = t-1 收盘→t 收盘(pct_change)."""
|
||
return closes.pct_change()
|
||
|
||
|
||
def zscore_cross(df: pd.DataFrame) -> pd.DataFrame:
|
||
"""逐日截面 z 分数(均值 0 方差 1; 单股日 std=0 → NaN).
|
||
|
||
inf 暴露=缺失(首跑实证: alpha83 部分日含 ±inf, pandas std 不跳 inf
|
||
→ 整行 z 陪葬 NaN)——先替换再标准化, 其余股正常.
|
||
"""
|
||
df = df.replace([np.inf, -np.inf], np.nan)
|
||
return df.sub(df.mean(axis=1), axis=0).div(df.std(axis=1), axis=0)
|
||
|
||
|
||
def _mean_daily_rank_corr(exposures: dict[str, pd.DataFrame]) -> pd.DataFrame:
|
||
"""日均 Spearman: 逐日截面秩相关再对日取均值(精确式, 非全池近似)."""
|
||
names = list(exposures)
|
||
acc = pd.DataFrame(0.0, index=names, columns=names)
|
||
cnt = pd.DataFrame(0, index=names, columns=names)
|
||
for d in exposures[names[0]].index:
|
||
cols = {}
|
||
for n in names:
|
||
s = exposures[n].loc[d].dropna()
|
||
if len(s) >= 3:
|
||
cols[n] = s.rank()
|
||
if len(cols) < 2:
|
||
continue
|
||
c = pd.DataFrame(cols).corr() # 秩的 Pearson = Spearman
|
||
acc = acc.add(c.fillna(0.0), fill_value=0.0)
|
||
cnt = cnt.add(c.notna().astype(int), fill_value=0)
|
||
return acc / cnt.replace(0, 1)
|
||
|
||
|
||
def family_representatives(exposures: dict[str, pd.DataFrame],
|
||
threshold: float = 0.7) -> dict[str, list[str]]:
|
||
"""|日均秩相关|>=threshold 连通分量分族; 代表=族内中位点(
|
||
与成员平均|corr|最大, 平手按名序). 返回 {代表名: 成员名表}."""
|
||
corr = _mean_daily_rank_corr(exposures)
|
||
names = sorted(exposures)
|
||
parent = {n: n for n in names}
|
||
|
||
def find(x: str) -> str:
|
||
while parent[x] != x:
|
||
parent[x] = parent[parent[x]]
|
||
x = parent[x]
|
||
return x
|
||
|
||
for i, a in enumerate(names):
|
||
for b in names[i + 1:]:
|
||
if abs(corr.loc[a, b]) >= threshold:
|
||
ra, rb = find(a), find(b)
|
||
if ra != rb:
|
||
parent[rb] = ra
|
||
groups: dict[str, list[str]] = {}
|
||
for n in names:
|
||
groups.setdefault(find(n), []).append(n)
|
||
out: dict[str, list[str]] = {}
|
||
for members in groups.values():
|
||
ms = sorted(members)
|
||
if len(ms) == 1:
|
||
out[ms[0]] = ms
|
||
continue
|
||
best, best_score = ms[0], -1.0
|
||
for m in ms:
|
||
score = sum(abs(corr.loc[m, o]) for o in ms if o != m)
|
||
if score > best_score + 1e-12:
|
||
best, best_score = m, score
|
||
out[best] = ms
|
||
return out
|
||
|
||
|
||
def vif(zdf: pd.DataFrame) -> dict[str, float]:
|
||
"""标准化设计阵 VIF = 相关阵逆对角线(奇异→inf)."""
|
||
X = zdf.dropna().to_numpy()
|
||
if X.ndim != 2 or X.shape[0] <= X.shape[1]:
|
||
return {c: float("inf") for c in zdf.columns}
|
||
corr = np.corrcoef(X, rowvar=False)
|
||
try:
|
||
inv = np.linalg.inv(corr)
|
||
except np.linalg.LinAlgError:
|
||
return {c: float("inf") for c in zdf.columns}
|
||
return {c: float(inv[i, i]) for i, c in enumerate(zdf.columns)}
|
||
|
||
|
||
RIDGE_LAMBDA = 1.0 # ridge 对照固定档(相关形设计阵 +λI, 尺度不变; spec
|
||
# §4.7-3 不调参. 09-27 审计 F-2: 旧实现 λI 直接加在
|
||
# 被 WLS 权重放大的原始正规阵上, √市值 量级(对角
|
||
# ~1e8)下相对收缩 ~1e-9=no-op, 对照通道恒空)
|
||
VIF_TRIGGER = 10.0 # 代表后仍 >10 → ridge 对照(spec §4.7-3)
|
||
|
||
|
||
def solve_day(X: np.ndarray, y: np.ndarray, w: np.ndarray,
|
||
ridge: float = 0.0) -> np.ndarray:
|
||
"""单日截面 WLS(X 第 0 列=截距全 1; ridge>0 截距不罚).
|
||
|
||
sqrt(市值) 权重由调用方传入(等权敏感性=传全 1, 固定双跑无旋钮).
|
||
ridge 惩罚随各列方差缩放: A+λ·diag(A) ≡ D^½(R+λI)D^½(R=相关形阵)——
|
||
权重整体缩放 c 时解不变(尺度不变), 系数统一收缩 1/(1+λ); 任意权重
|
||
量级(等权/√市值)下对照通道同等生效.
|
||
"""
|
||
sw = np.sqrt(np.maximum(w, 0.0))
|
||
Xw, yw = X * sw[:, None], y * sw
|
||
A, b = Xw.T @ Xw, Xw.T @ yw
|
||
if ridge > 0.0:
|
||
pen = np.diag(A) * ridge
|
||
pen[0] = 0.0 # 截距不罚
|
||
A = A + np.diag(pen)
|
||
return np.linalg.lstsq(A, b, rcond=None)[0]
|
||
|
||
|
||
def _tradable(vt: str) -> bool:
|
||
"""主板+创业板(剔科创 68/北交 4,8/920).
|
||
|
||
与 sanguo_portfolio/strategies/factor_topn._filter_kcbj_keep_chinext
|
||
同规则的 6 行双胞胎——factor 域不反向依赖 portfolio, 改规则双改.
|
||
"""
|
||
code = vt.split(".", 1)[0]
|
||
if not code:
|
||
return False
|
||
return not (code[0] in ("4", "8") or code[:2] in ("68", "92"))
|
||
|
||
|
||
def holdings_schedule(section: pd.DataFrame, top_n: int = 200,
|
||
reb_days: int = 63) -> list[tuple[pd.Timestamp, list[str]]]:
|
||
"""模型组合持仓日程: 每 reb_days 个交易日按截面值重选 TopN(硬切
|
||
band=1.0; 值大=好, 同 factor_topn 方向约定). 执行级过滤(ST/涨跌停/
|
||
次新)不重建——模型→实盘差距归偏差日报半边."""
|
||
out: list[tuple[pd.Timestamp, list[str]]] = []
|
||
held: list[str] | None = None
|
||
for i, d in enumerate(section.index):
|
||
if i % reb_days == 0 or held is None:
|
||
row = section.loc[d].dropna()
|
||
row = row[[c for c in row.index if _tradable(c)]]
|
||
held = list(row.sort_values(ascending=False).index[:top_n])
|
||
if held:
|
||
out.append((d, list(held)))
|
||
return out
|
||
|
||
|
||
def portfolio_returns(schedule: list[tuple[pd.Timestamp, list[str]]],
|
||
rets: pd.DataFrame) -> tuple[pd.Series, pd.DataFrame]:
|
||
"""组合日收益+逐日权重(实现日口径). 同持期间权重随前日收益漂移
|
||
(买入持有); 停牌 NaN 收益按 0(价格冻结, 口径披露)."""
|
||
r_p = pd.Series(0.0, index=[d for d, _ in schedule])
|
||
whist = pd.DataFrame(0.0, index=r_p.index, columns=rets.columns)
|
||
held: list[str] | None = None
|
||
w = pd.Series(dtype=float)
|
||
prev_ret: pd.Series | None = None
|
||
for d, h in schedule:
|
||
if h != held:
|
||
held = h
|
||
w = pd.Series(0.0, index=rets.columns)
|
||
w[h] = 1.0 / len(h)
|
||
elif prev_ret is not None:
|
||
growth = w * (1.0 + prev_ret.reindex(w.index).fillna(0.0))
|
||
w = growth / growth.sum()
|
||
r_p.loc[d] = float((w * rets.loc[d].reindex(w.index).fillna(0.0)).sum())
|
||
whist.loc[d] = w
|
||
prev_ret = rets.loc[d]
|
||
return r_p, whist
|
||
|
||
|
||
# ---------------- 归因主循环(P3 §4.7-1/2/3) ----------------
|
||
|
||
SIZE_LEAK_THRESHOLD = 0.02 # |size 桶窗口 cum| > 2% = 泄漏报警(拍板值)
|
||
MIN_FLIP_MAGNITUDE = 0.001 # 月度 |贡献| 低于此不判翻转(噪声地板)
|
||
|
||
|
||
def size_leak(res: dict) -> tuple[bool, float]:
|
||
"""sanity 锚点: 组合层已 size 中性化 ⇒ size 桶累计贡献应≈0."""
|
||
cum = float((1.0 + res["contributions"]["size"]).prod() - 1.0)
|
||
return abs(cum) > SIZE_LEAK_THRESHOLD, cum
|
||
|
||
|
||
_SUFFIX_MAP = {".SH": ".SSE", ".SZ": ".SZSE", ".BJ": ".BSE"}
|
||
|
||
|
||
def _norm_symbol(s: str) -> str:
|
||
"""通用 .SH/.SZ/.BJ 形(SW 等官方交付) → vnpy vt_symbol 形."""
|
||
for suf, vt in _SUFFIX_MAP.items():
|
||
if s.endswith(suf):
|
||
return s[: -len(suf)] + vt
|
||
return s
|
||
|
||
|
||
def _industry_dummies(industry_map: pd.DataFrame) -> tuple[dict, list]:
|
||
"""行业码映射 + 哑变量码表(去最大桶为基线; 未知 0 单独一列).
|
||
|
||
表 symbol 先经 _norm_symbol 归一(data 交付 .SH/.SZ 通用形态).
|
||
"""
|
||
ind_code = { _norm_symbol(s): int(c) for s, c in
|
||
industry_map.set_index("symbol")["industry_code"].items() }
|
||
codes = sorted(set(ind_code.values()))
|
||
baseline = (max(codes, key=lambda c: sum(1 for v in ind_code.values()
|
||
if v == c)) if codes else None)
|
||
return ind_code, [c for c in codes if c != baseline]
|
||
|
||
|
||
def _day_design(z_row: dict, zs_row: pd.Series, f_ret_next: pd.DataFrame,
|
||
t, ind_code: dict, dummy_codes: list):
|
||
"""单日设计阵装配(主循环与 ridge 通道共用).
|
||
|
||
z_row/zs_row = 当日已切片的暴露行(调用方负责逐日切片一次).
|
||
返回 (X, y, rows_sym) 或 None(样本不足 max(3×列数, 50))."""
|
||
rows_x, rows_y, rows_sym = [], [], []
|
||
base = f_ret_next.loc[t].dropna()
|
||
for s in base.index:
|
||
xs = [z_row[k].get(s, np.nan) for k in z_row]
|
||
xs.append(zs_row.get(s, np.nan))
|
||
if not all(np.isfinite(v) for v in xs):
|
||
continue
|
||
rows_x.append([1.0] + xs + [1.0 if ind_code[s] == c else 0.0
|
||
for c in dummy_codes])
|
||
rows_y.append(float(base[s]))
|
||
rows_sym.append(s)
|
||
n_cols = 1 + len(z_row) + 1 + len(dummy_codes)
|
||
if len(rows_x) < max(3 * n_cols, 50):
|
||
return None
|
||
return np.array(rows_x), np.array(rows_y), rows_sym
|
||
|
||
|
||
def attribute(section: pd.DataFrame, exposures: dict[str, pd.DataFrame],
|
||
size_log: pd.DataFrame, industry_map: pd.DataFrame,
|
||
rets: pd.DataFrame, top_n: int = 200, reb_days: int = 63,
|
||
weighting: str = "sqrt_mktcap") -> dict:
|
||
"""归因主循环: 逐日 FM 回归 → 因子日收益 → 组合贡献分解.
|
||
|
||
时序: t 日暴露解释 t+1 实现收益; 贡献落 t+1(权重=当日开盘权重,
|
||
暴露=t 日值); 逐股分解 r_i=拟合+残差精确 ⇒ 恒等式对任意持仓成立;
|
||
持仓中暴露 NaN 的个股拟合记 0(收益全入特质, 占比披露).
|
||
contributions 不含 ind_ 明细列(聚合进 industry 桶, 消费端 schema
|
||
干净); 明细保留在 factor_returns.
|
||
"""
|
||
fams = family_representatives(exposures)
|
||
reps = list(fams)
|
||
z = {k: zscore_cross(v) for k, v in exposures.items()}
|
||
z_size = zscore_cross(size_log)
|
||
dates = [d for d in section.index if d in rets.index]
|
||
f_ret_next = rets.shift(-1)
|
||
# 行业表契约(data 09-26): 表内只含 SW 已映射股(.SH/.SZ 通用形态, 入口
|
||
# 归一到 vt 形), join 不中 → 0 码独立桶(显式哑变量列; 0 桶最大时自适
|
||
# 应成基线)——不静默丢样本
|
||
universe = sorted(set(rets.columns) | set(section.columns))
|
||
known = industry_map.set_index("symbol")["industry_code"]
|
||
known.index = [_norm_symbol(s) for s in known.index]
|
||
industry_map = pd.DataFrame({
|
||
"symbol": universe,
|
||
"industry_code": known.reindex(universe).fillna(0).astype(int).to_numpy(),
|
||
})
|
||
ind_code, dummy_codes = _industry_dummies(industry_map)
|
||
|
||
sched = holdings_schedule(section, top_n=top_n, reb_days=reb_days)
|
||
_, whist = portfolio_returns(sched, rets)
|
||
whist = whist.reindex(dates)
|
||
|
||
k_rep = len(reps)
|
||
fr_cols = ["market"] + reps + ["size"] + [f"ind_{c}" for c in dummy_codes]
|
||
contributions = pd.DataFrame(
|
||
0.0, index=dates, columns=["market"] + reps + ["size", "industry", "specific"])
|
||
factor_returns = pd.DataFrame(np.nan, index=dates, columns=fr_cols)
|
||
port_ret = pd.Series(0.0, index=dates)
|
||
expo_sum: dict[str, list[float]] = {k: [] for k in list(reps) + ["size"]}
|
||
r2s: list[float] = []
|
||
nan_hold_counts = hold_counts = 0
|
||
unknown = mapped = 0
|
||
mkt_raw = np.exp(size_log) # sqrt(市值) WLS 权的市值
|
||
|
||
for i, t in enumerate(dates[:-1]): # 最后一日无 t+1 收益
|
||
t1 = dates[i + 1]
|
||
w_day = whist.loc[t1]
|
||
held = [s for s in w_day.index[w_day > 0]]
|
||
if not held:
|
||
continue
|
||
hold_counts += len(held)
|
||
z_row = {k: z[k].loc[t] for k in reps}
|
||
zs_row = z_size.loc[t]
|
||
nan_hold_counts += sum(
|
||
1 for s in held
|
||
if not all(np.isfinite(z_row[k].get(s, np.nan)) for k in reps)
|
||
or not np.isfinite(zs_row.get(s, np.nan)))
|
||
got = _day_design(z_row, zs_row, f_ret_next, t, ind_code, dummy_codes)
|
||
if got is None:
|
||
continue
|
||
X, y, rows_sym = got
|
||
unknown += sum(1 for s in rows_sym if ind_code[s] == 0)
|
||
mapped += len(rows_sym)
|
||
if weighting == "sqrt_mktcap":
|
||
w = np.sqrt(mkt_raw.loc[t].reindex(rows_sym).fillna(0.0)
|
||
.clip(lower=1e-8).to_numpy())
|
||
else:
|
||
w = np.ones(len(rows_sym))
|
||
beta = solve_day(X, y, w)
|
||
sw = np.sqrt(w)
|
||
resid = y - X @ beta
|
||
ss_res = float(np.sum(sw ** 2 * resid ** 2))
|
||
ybar = float(np.average(y, weights=w))
|
||
ss_tot = float(np.sum(sw ** 2 * (y - ybar) ** 2))
|
||
r2s.append(1.0 - ss_res / ss_tot if ss_tot > 0 else np.nan)
|
||
factor_returns.loc[t] = beta
|
||
|
||
# 贡献(实现日 t1): E_k = Σ w_i(t1)·x_i(t); ind_ 聚合进 industry 桶
|
||
w_held = w_day.reindex(rows_sym).fillna(0.0)
|
||
expo = (X * w_held.to_numpy()[:, None]).sum(axis=0)
|
||
contrib_vals = expo * beta
|
||
row = {"market": float(contrib_vals[0])}
|
||
for j, k in enumerate(reps):
|
||
row[k] = float(contrib_vals[1 + j])
|
||
expo_sum[k].append(float(expo[1 + j]))
|
||
row["size"] = float(contrib_vals[1 + k_rep])
|
||
expo_sum["size"].append(float(expo[1 + k_rep]))
|
||
row["industry"] = float(contrib_vals[1 + k_rep + 1:].sum())
|
||
# specific: 持仓全体 r_i − 拟合_i(样本外股拟合=0 → 全入特质);
|
||
# NaN 收益=价格冻结口径 0(显式 isfinite, 勿用 or 惯用语——NaN truthy)
|
||
fitted = X @ beta
|
||
r_t1 = rets.loc[t1]
|
||
spec = 0.0
|
||
rows_sym_set = set(rows_sym)
|
||
for j, s in enumerate(rows_sym):
|
||
wgt = float(w_day.get(s, 0.0))
|
||
if wgt > 0:
|
||
rv = r_t1.get(s, np.nan)
|
||
spec += wgt * ((0.0 if not np.isfinite(rv) else float(rv))
|
||
- fitted[j])
|
||
for s in held:
|
||
if s not in rows_sym_set:
|
||
rv = r_t1.get(s, np.nan)
|
||
spec += float(w_day[s]) * (0.0 if not np.isfinite(rv)
|
||
else float(rv))
|
||
row["specific"] = spec
|
||
contributions.loc[t1] = row
|
||
port_ret.loc[t1] = float((w_day * r_t1.reindex(w_day.index)
|
||
.fillna(0.0)).sum())
|
||
|
||
# VIF(逐源 z 分数 stack 须 inner 对齐——各源 NaN 域不同, 错位=静默垃圾)
|
||
stacked = pd.concat({k: z[k].stack() for k in reps}, axis=1, join="inner")
|
||
v = vif(stacked)
|
||
# ridge 对照: 仅当代表 VIF 超阈(spec §4.7-3), 窗口累计符号翻转才报
|
||
ridge_flip: list[str] = []
|
||
if any(vv > VIF_TRIGGER for vv in v.values()):
|
||
ols_m = contributions[reps].groupby(
|
||
contributions.index.to_period("M")).sum()
|
||
ridge_m = _ridge_monthly(z, reps, z_size, f_ret_next, dates,
|
||
ind_code, dummy_codes, mkt_raw, weighting)
|
||
for k in reps:
|
||
a, b = float(ols_m[k].sum()), float(ridge_m[k].sum())
|
||
if (abs(a) > MIN_FLIP_MAGNITUDE and abs(b) > MIN_FLIP_MAGNITUDE
|
||
and np.sign(a) != np.sign(b)):
|
||
ridge_flip.append(k)
|
||
return {
|
||
"factor_returns": factor_returns, "contributions": contributions,
|
||
"portfolio_ret": port_ret, "families": fams,
|
||
"exposure_avg": {k: float(np.nanmean(vv)) for k, vv in expo_sum.items()},
|
||
"r2_mean": float(np.nanmean(r2s)) if r2s else float("nan"),
|
||
"vif": v, "ridge_flip": ridge_flip,
|
||
"nan_exposure_holding_share": nan_hold_counts / max(hold_counts, 1),
|
||
"unknown_industry_share": unknown / max(mapped, 1),
|
||
}
|
||
|
||
|
||
def _ridge_monthly(z, reps, z_size, f_ret_next, dates, ind_code,
|
||
dummy_codes, mkt_raw, weighting) -> pd.DataFrame:
|
||
"""ridge 对照通道(λ=RIDGE_LAMBDA 固定): 复用 _day_design, 只换估计."""
|
||
daily = pd.DataFrame(0.0, index=dates[:-1], columns=reps)
|
||
for t in dates[:-1]:
|
||
z_row = {k: z[k].loc[t] for k in reps}
|
||
got = _day_design(z_row, z_size.loc[t], f_ret_next, t,
|
||
ind_code, dummy_codes)
|
||
if got is None:
|
||
continue
|
||
X, y, rows_sym = got
|
||
if weighting == "sqrt_mktcap":
|
||
w = np.sqrt(mkt_raw.loc[t].reindex(rows_sym).fillna(0.0)
|
||
.clip(lower=1e-8).to_numpy())
|
||
else:
|
||
w = np.ones(len(rows_sym))
|
||
beta = solve_day(X, y, w, ridge=RIDGE_LAMBDA)
|
||
daily.loc[t] = beta[1:1 + len(reps)]
|
||
return daily.groupby(daily.index.to_period("M")).sum()
|
||
|
||
|
||
# ---------------- 月度汇总 / 报告 / CLI(P3 §4.7-4/5) ----------------
|
||
|
||
def monthly_rows(res: dict, families: dict[str, list[str]]) -> list[dict]:
|
||
"""月度行(schema verbatim): 代表行(factor=代表名)→market→size→
|
||
industry→specific(三桶殿后). cum_ret=Π(1+日贡献)−1 全窗累计."""
|
||
contrib = res["contributions"]
|
||
order = [k for k in contrib.columns
|
||
if k not in ("market", "size", "industry", "specific")] \
|
||
+ ["market", "size", "industry", "specific"]
|
||
months = contrib.index.to_period("M")
|
||
rows: list[dict] = []
|
||
for k in order:
|
||
cum = float((1.0 + contrib[k].fillna(0.0)).prod() - 1.0)
|
||
for m, v in contrib[k].groupby(months).sum().items():
|
||
rows.append({"factor": k, "month": str(m),
|
||
"contribution_ret": float(v), "cum_ret": cum,
|
||
"exposure_avg": float(res["exposure_avg"].get(k, 0.0))})
|
||
return rows
|
||
|
||
|
||
def rolling_ir(series: pd.Series, months: int) -> float | None:
|
||
"""滚动窗年化 IR(mean/std×√252), 样本不足或 std=0 → None."""
|
||
s = series.dropna()
|
||
if len(s) < months * 15: # 月≈21 交易日, 15=下限容忍
|
||
return None
|
||
win = s.iloc[-months * 21:]
|
||
sd = float(win.std())
|
||
if not np.isfinite(sd) or sd == 0.0:
|
||
return None
|
||
return float(win.mean() / sd * np.sqrt(252.0))
|
||
|
||
|
||
def build_report(res: dict, families: dict[str, list[str]], as_of: str,
|
||
start: str, end: str) -> dict:
|
||
from datetime import datetime
|
||
leak, size_cum = size_leak(res)
|
||
return {
|
||
"as_of": as_of,
|
||
"generated_at": datetime.now().isoformat(timespec="seconds"),
|
||
"window": {"start": start, "end": end},
|
||
"factorContribution": monthly_rows(res, families),
|
||
"families": families,
|
||
"rolling": {k: {"ir_12m": rolling_ir(res["contributions"][k], 12),
|
||
"ir_24m": rolling_ir(res["contributions"][k], 24)}
|
||
for k in res["contributions"].columns},
|
||
"diagnostics": {
|
||
"r2_mean": res["r2_mean"], "vif": res["vif"],
|
||
"vif_max": max(res["vif"].values()) if res["vif"] else None,
|
||
"ridge_flip": res["ridge_flip"], "size_leak": leak,
|
||
"size_cum": size_cum,
|
||
"nan_exposure_holding_share": res["nan_exposure_holding_share"],
|
||
"unknown_industry_share": res["unknown_industry_share"],
|
||
"weighting": "sqrt_mktcap",
|
||
"industry_note": "行业=最新静态成员回溯(迁移缓慢, 业界同款); "
|
||
"表未覆盖股=0 码独立桶",
|
||
},
|
||
}
|
||
|
||
|
||
def main(argv: list[str] | None = None) -> int:
|
||
"""CLI: python -m sanguo_factor.attribution --sections-dir D --factor F ...
|
||
|
||
产物=<out-dir>/attribution_<end>.json(tmp+os.replace 原子写).
|
||
"""
|
||
import argparse
|
||
import json
|
||
import os
|
||
import tempfile
|
||
|
||
ap = argparse.ArgumentParser(description="因子贡献分解(日度 Fama-MacBeth)")
|
||
ap.add_argument("--sections-dir", required=True)
|
||
ap.add_argument("--factor", required=True, help="组合截面因子名(持仓重建用)")
|
||
ap.add_argument("--sources", default=None,
|
||
help="逗号分隔源截面名(缺省=QUANT_SOURCES)")
|
||
ap.add_argument("--size-name", default="size_log")
|
||
ap.add_argument("--start", required=True)
|
||
ap.add_argument("--end", required=True)
|
||
ap.add_argument("--vnpy-db", required=True)
|
||
ap.add_argument("--industry-path", required=True)
|
||
ap.add_argument("--out-dir", required=True)
|
||
ap.add_argument("--top-n", type=int, default=200)
|
||
ap.add_argument("--reb-days", type=int, default=63)
|
||
ap.add_argument("--weighting", choices=["sqrt_mktcap", "equal"],
|
||
default="sqrt_mktcap")
|
||
args = ap.parse_args(argv)
|
||
|
||
from .composite_library import QUANT_SOURCES
|
||
sources = (args.sources.split(",") if args.sources
|
||
else [n for n, _ in QUANT_SOURCES])
|
||
names = sources + [args.size_name, args.factor]
|
||
secs = load_sections(args.sections_dir, names)
|
||
section = secs[args.factor]
|
||
section = section.loc[(section.index >= pd.Timestamp(args.start))
|
||
& (section.index <= pd.Timestamp(args.end))]
|
||
symbols = sorted(set().union(*[set(df.columns) for df in secs.values()]))
|
||
closes = load_closes(args.vnpy_db, symbols, args.start, args.end)
|
||
closes.index = pd.to_datetime(closes.index)
|
||
rets = forward_returns(closes)
|
||
industry_map = pd.read_parquet(args.industry_path)
|
||
res = attribute(section, {k: secs[k] for k in sources},
|
||
secs[args.size_name], industry_map, rets,
|
||
top_n=args.top_n, reb_days=args.reb_days,
|
||
weighting=args.weighting)
|
||
fams = {k: [s for s in v if s in sources] for k, v in res["families"].items()}
|
||
fams = {k: v for k, v in fams.items() if v}
|
||
report = build_report(res, fams, args.end, args.start, args.end)
|
||
os.makedirs(args.out_dir, exist_ok=True)
|
||
final = os.path.join(args.out_dir, f"attribution_{args.end}.json")
|
||
fd, tmp = tempfile.mkstemp(suffix=".json", dir=args.out_dir)
|
||
os.close(fd)
|
||
with open(tmp, "w", encoding="utf-8") as f:
|
||
json.dump(report, f, ensure_ascii=False, indent=2)
|
||
os.replace(tmp, final)
|
||
print(f"[attribution] {final} rows={len(report['factorContribution'])} "
|
||
f"r2={report['diagnostics']['r2_mean']:.3f} "
|
||
f"size_leak={report['diagnostics']['size_leak']}")
|
||
return 0
|
||
|
||
|
||
if __name__ == "__main__":
|
||
raise SystemExit(main())
|