# -*- coding: utf-8 -*-
"""商品价格 → 正股:怎么用现货信号炒股赚钱 (可复现回测)。

问题:追踪 钨/锑/稀土 这些商品价格, 到底能不能帮你在 A 股赚钱?
核心是「谁领先谁」——商品先动、正股后跟(有 edge), 还是正股已提前定价(没用)?

无期货的小金属没有历史现货序列, 无法直接回测。所以在**有期货、有长历史**的
商品上验证这条普适规律(铜/铝/锌/钢/煤/油/金银/锂/硅/纯碱/尿素/生猪/PTA),
每个商品配纯正 A 股正股。结论再外推到 /commodities 的无期货篮子(机制相同,
且散户驱动、覆盖更薄 → lead-lag 大概率更强, 但需 tracker 积累后另测)。

数据: 商品期货主力连续 = akshare futures_main_sina; 正股日线 = tushare pro.daily
(pct_chg 复权). 月频。成本每次换仓 0.2% 双边。train<2022 / test>=2022 OOS。
A 股只报长仓口径。
"""
import sys, json, time, warnings
warnings.filterwarnings("ignore")
sys.path.insert(0, "/root/cb-allotment/scripts")
import numpy as np, pandas as pd
import akshare as ak
from tushare_client import get_pro

pro = get_pro()
START, END = "20160101", "20260701"
import os
CACHE = "/tmp/cs_panels.parquet"; CACHE2 = "/tmp/cs_panels_stk.parquet"

# 商品期货主力(sina符号) -> (中文, [纯正正股 ts_code])
PAIRS = {
    "CU0": ("铜",   ["600362.SH", "000630.SZ", "000878.SZ"]),
    "AL0": ("铝",   ["000807.SZ", "601600.SH", "000933.SZ"]),
    "ZN0": ("锌",   ["600497.SH", "600961.SH"]),
    "RB0": ("螺纹钢", ["000932.SZ", "002110.SZ"]),
    "J0":  ("焦炭",  ["601015.SH", "600985.SH"]),
    "JM0": ("焦煤",  ["000983.SZ", "601666.SH"]),
    "SC0": ("原油",  ["601857.SH", "600938.SH"]),
    "AU0": ("黄金",  ["600547.SH", "600489.SH", "601899.SH"]),
    "AG0": ("白银",  ["000603.SZ", "000426.SZ"]),
    "LC0": ("碳酸锂", ["002466.SZ", "002460.SZ"]),
    "SI0": ("工业硅", ["603260.SH"]),
    "FG0": ("玻璃",  ["601636.SH"]),
    "SA0": ("纯碱",  ["600409.SH", "000683.SZ"]),
    "UR0": ("尿素",  ["600426.SH", "000422.SZ"]),
    "LH0": ("生猪",  ["002714.SZ", "300498.SZ"]),
    "TA0": ("PTA",  ["002493.SZ", "601233.SH"]),
    "V0":  ("PVC",  ["002092.SZ"]),
}

def month_key(s):
    return pd.to_datetime(s).dt.to_period("M")

if os.path.exists(CACHE) and os.path.exists(CACHE2) and "--refresh" not in sys.argv:
    COM = pd.read_parquet(CACHE); STK = pd.read_parquet(CACHE2)
    COM.index = pd.PeriodIndex(COM.index, freq="M"); STK.index = pd.PeriodIndex(STK.index, freq="M")
else:
  # ---- 商品月收益 ----
  com_m = {}
  for sym, (nm, _) in PAIRS.items():
    try:
        d = ak.futures_main_sina(symbol=sym, start_date=START, end_date=END)
        d = d.rename(columns={"日期": "date", "收盘价": "close"})[["date", "close"]]
        d["m"] = month_key(d["date"]); mc = d.groupby("m")["close"].last()
        com_m[sym] = mc.pct_change()
        print(f"  期货 {nm}({sym}): {len(mc)} 月 {mc.index.min()}..{mc.index.max()}", file=sys.stderr)
    except Exception as e:
        print(f"  期货 {nm}({sym}) FAIL {str(e)[:60]}", file=sys.stderr)
    time.sleep(0.4)
  COM = pd.DataFrame(com_m)  # index=月, col=期货符号

  # ---- 正股月收益 ----
  stk_m = {}
  allcodes = sorted({c for _, codes in PAIRS.values() for c in codes})
  for c in allcodes:
    try:
        d = pro.daily(ts_code=c, start_date=START, end_date=END)
        if d is None or not len(d): continue
        d = d.sort_values("trade_date")
        d["m"] = month_key(d["trade_date"])
        mr = d.groupby("m").apply(lambda x: (1 + x["pct_chg"] / 100).prod() - 1)
        stk_m[c] = mr
    except Exception as e:
        print(f"  股 {c} FAIL {str(e)[:50]}", file=sys.stderr)
    time.sleep(0.25)
  STK = pd.DataFrame(stk_m)
  COM.index = COM.index.astype(str); STK.index = STK.index.astype(str)
  COM.to_parquet(CACHE); STK.to_parquet(CACHE2)
  COM.index = pd.PeriodIndex(COM.index, freq="M"); STK.index = pd.PeriodIndex(STK.index, freq="M")

allcodes = sorted({c for _, codes in PAIRS.values() for c in codes})
# 对齐月区间
idx = COM.index.intersection(STK.index)
COM, STK = COM.loc[idx], STK.loc[idx]
TRAIN = pd.Period("2022-01", "M")

def ann_stats(rets):
    r = pd.Series(rets).dropna()
    if len(r) < 6: return {}
    cum = (1 + r).prod()
    cagr = cum ** (12 / len(r)) - 1
    vol = r.std() * np.sqrt(12)
    shp = (r.mean() * 12) / vol if vol > 0 else 0
    nav = (1 + r).cumprod()
    mdd = (nav / nav.cummax() - 1).min()
    return {"cagr": round(cagr * 100, 1), "sharpe": round(shp, 2),
            "mdd": round(mdd * 100, 1), "n": int(len(r))}

# ============ 测试 A:lead-lag 回归(商品领先正股吗) ============
# 池化: 正股当月收益 ~ 商品同月收益 + 商品上月收益 (+ 正股上月收益)
rows = []
for sym, (nm, codes) in PAIRS.items():
    if sym not in COM.columns: continue
    c_now = COM[sym]; c_lag = COM[sym].shift(1)
    for code in codes:
        if code not in STK.columns: continue
        s = STK[code]; s_lag = s.shift(1)
        df = pd.concat([s, c_now, c_lag, s_lag], axis=1).dropna()
        df.columns = ["s", "cn", "cl", "sl"]
        for _, r in df.iterrows():
            rows.append(r.to_dict())
P = pd.DataFrame(rows)
# 多元回归 s ~ cn + cl + sl
X = np.column_stack([np.ones(len(P)), P["cn"], P["cl"], P["sl"]])
y = P["s"].values
beta, *_ = np.linalg.lstsq(X, y, rcond=None)
resid = y - X @ beta
se = np.sqrt(np.diag(np.linalg.inv(X.T @ X) * (resid @ resid) / (len(y) - 4)))
tvals = beta / se
leadlag = {"n_obs": int(len(P)),
           "beta_同月商品": round(beta[1], 3), "t_同月": round(tvals[1], 1),
           "beta_上月商品": round(beta[2], 3), "t_上月": round(tvals[2], 1),
           "beta_上月正股": round(beta[3], 3), "t_正股动量": round(tvals[3], 1)}

# ============ 测试 B:商品趋势择时 (核心赚钱逻辑) ============
# 每只正股: 下月持有 当且仅当 其商品 3 月动量>0, 否则空仓(现金). 对比 买入持有.
def build_strategy(mom_win=3, filt=True):
    port = []  # 月度组合收益
    for m_i in range(1, len(STK)):
        m = STK.index[m_i]
        rr = []
        for sym, (nm, codes) in PAIRS.items():
            if sym not in COM.columns: continue
            # 商品截至上月末的 mom_win 月动量
            past = COM[sym].iloc[max(0, m_i - mom_win):m_i]
            sig = (1 + past.fillna(0)).prod() - 1
            hold = (not filt) or (pd.notna(sig) and sig > 0)
            if not hold: continue
            for code in codes:
                if code in STK.columns and pd.notna(STK[code].iloc[m_i]):
                    rr.append(STK[code].iloc[m_i])
        port.append((m, np.mean(rr) if rr else 0.0, len(rr)))
    return pd.DataFrame(port, columns=["m", "ret", "nheld"]).set_index("m")

bh = build_strategy(filt=False)     # 买入持有全篮子
tf = build_strategy(3, True)        # 商品3月动量>0 才持有
# 成本: 择时版按持仓变化估算换手 (进出各 0.2%*仓位变动)。粗略: 每月 nheld 变化比例 * 0.2%
def apply_cost(df, base_n):
    prev = None; net = []
    for m, row in df.iterrows():
        turn = 1.0 if prev is None else abs(row["nheld"] - prev) / max(base_n, 1)
        net.append(row["ret"] - 0.002 * min(turn, 1.0))
        prev = row["nheld"]
    return pd.Series(net, index=df.index)
base_n = len(allcodes)
tf_net = apply_cost(tf, base_n)

def split_stats(series):
    s = pd.Series(series)
    return {"all": ann_stats(s),
            "train": ann_stats(s[s.index < TRAIN]),
            "test": ann_stats(s[s.index >= TRAIN])}

# ============ 测试 C:横截面商品动量倾斜 (长仓超额, 变现 0.08 lag) ============
# 每月按 商品上月收益(或N月动量) 排序, 只买 上半区商品的正股 vs 等权全篮子基准。
def cross_section(mom_win=1):
    excess = []; eq = []; tilt = []
    for m_i in range(mom_win, len(STK)):
        m = STK.index[m_i]
        # 各商品截至上月的动量
        sig = {}
        for sym in COM.columns:
            past = COM[sym].iloc[m_i - mom_win:m_i]
            sig[sym] = (1 + past.fillna(0)).prod() - 1 if past.notna().any() else np.nan
        ss = pd.Series(sig).dropna()
        if len(ss) < 6: continue
        hi = set(ss[ss >= ss.median()].index)
        hi_r, all_r = [], []
        for sym, (nm, codes) in PAIRS.items():
            if sym not in COM.columns: continue
            for code in codes:
                if code in STK.columns and pd.notna(STK[code].iloc[m_i]):
                    all_r.append(STK[code].iloc[m_i])
                    if sym in hi: hi_r.append(STK[code].iloc[m_i])
        if hi_r and all_r:
            eq.append(np.mean(all_r)); tilt.append(np.mean(hi_r))
            excess.append(np.mean(hi_r) - np.mean(all_r))
    di = STK.index[mom_win:mom_win + len(excess)]
    return pd.Series(tilt, index=di), pd.Series(eq, index=di), pd.Series(excess, index=di)

cs_variants = {}
for w in (1, 3, 6):
    tlt, eq, exc = cross_section(w)
    # 倾斜组扣成本(每月部分换手, 近似 0.2%*0.5)
    tlt_net = tlt - 0.001
    cs_variants[f"mom{w}"] = {
        "tilt_net": split_stats(tlt_net), "benchmark_eq": split_stats(eq),
        "excess_ann": round(exc.mean() * 12 * 100, 1),
        "excess_t": round(exc.mean() / exc.std() * np.sqrt(len(exc)), 1),
        "excess_train": round(exc[exc.index < TRAIN].mean() * 12 * 100, 1),
        "excess_test": round(exc[exc.index >= TRAIN].mean() * 12 * 100, 1),
    }

result = {
    "universe": {"n_commodities": int(COM.shape[1]), "n_stocks": int(STK.shape[1]),
                 "months": [str(idx.min()), str(idx.max())], "n_months": int(len(idx))},
    "leadlag": leadlag,
    "buy_hold": split_stats(bh["ret"]),
    "trend_filter_gross": split_stats(tf["ret"]),
    "trend_filter_net": split_stats(tf_net),
    "cross_section": cs_variants,
    "avg_held_filtered": round(tf["nheld"].mean(), 1),
    "avg_held_total": base_n,
    "nav_bh": [round(x, 4) for x in (1 + bh["ret"]).cumprod().tolist()],
    "nav_tf": [round(x, 4) for x in (1 + tf_net).cumprod().tolist()],
    "nav_dates": [str(m) for m in bh.index],
}
# ============ 测试 D:商品创新高事件 (regime 切换, 非连续趋势) ============
# 用商品价格水平(把月收益累乘回价格). breakout = 收盘创 win 月新高. 比较正股 forward。
comlvl = {}
for sym in COM.columns:
    comlvl[sym] = (1 + COM[sym].fillna(0)).cumprod()
CL = pd.DataFrame(comlvl)
def breakout_event(win=12, fwd=1):
    hit, base = [], []
    for sym, (nm, codes) in PAIRS.items():
        if sym not in CL.columns: continue
        lvl = CL[sym]
        rollmax = lvl.shift(1).rolling(win).max()
        is_bo = lvl > rollmax  # 本月创 win 月新高
        for code in codes:
            if code not in STK.columns: continue
            f = STK[code]  # 月收益
            for i in range(len(f) - fwd):
                if not bool(is_bo.iloc[i]) and not True: continue
                fwd_cum = (1 + f.iloc[i + 1:i + 1 + fwd]).prod() - 1  # 下月起 fwd 月累计
                if pd.isna(fwd_cum): continue
                (hit if bool(is_bo.iloc[i]) else base).append(fwd_cum)
    hit, base = np.array(hit), np.array(base)
    # 年化换算
    def annz(x): return (np.nanmean(x)) * (12 / fwd) * 100
    return {"n_breakout": int(len(hit)), "n_base": int(len(base)),
            "fwd_months": fwd, "win": win,
            "breakout_fwd_ann": round(annz(hit), 1),
            "base_fwd_ann": round(annz(base), 1),
            "edge_ann": round(annz(hit) - annz(base), 1),
            "hit_winrate": round((hit > 0).mean() * 100, 0)}
result["breakout"] = {f"{w}m_hi_{f}m_fwd": breakout_event(w, f)
                      for w, f in [(12, 1), (12, 3), (6, 1), (24, 3)]}

# 每商品:同月 vs 次月 beta (哪些商品最"迟钝"= 正股跟得慢 = 现货信号最有用)
per = {}
for sym, (nm, codes) in PAIRS.items():
    if sym not in COM.columns: continue
    cn, cl = COM[sym], COM[sym].shift(1)
    for code in codes:
        if code not in STK.columns: continue
        df = pd.concat([STK[code], cn, cl], axis=1).dropna(); df.columns = ["s", "cn", "cl"]
        if len(df) < 24: continue
        b_now = np.polyfit(df["cn"], df["s"], 1)[0]
        b_lag = np.polyfit(df["cl"], df["s"], 1)[0]
        per.setdefault(nm, []).append((code, round(b_now, 2), round(b_lag, 2)))
result["per_commodity_beta"] = {k: v for k, v in per.items()}

json.dump(result, open("/root/cb-allotment/my-app/public/reports/commodity-stocks/backtest_result.json", "w"),
          ensure_ascii=False, indent=1)
print(json.dumps({k: result[k] for k in ["universe", "leadlag", "buy_hold", "trend_filter_net", "cross_section"]},
                 ensure_ascii=False, indent=1))