backtest/engine.py

"""金仓 v2.1 回测引擎：信号复现 + 决策表状态机（纯标准库）。

规则与口径见《金仓v2-全历史回测计划书.md》v2.1（预注册，参数冻结）：
  σ    RV-EWMA6 月度递推，月度口径（逐行复刻 app/model.py sigma_rv / rv_sq）
  q    exp-03 W3 walk-forward 分数序列（exp/results/exp03_regime.json）的
       逐点「前序排序分位」（≤ 比较，不含自身）；2021-04 起可用，此前 q=None
  pace WGC 季度 cbTonnes；季度结束后第 31 个自然日（end+30d）起可见

决策表优先级（计划书 §2.3）：
  1 逻辑退出(月度, pace<102 → 卖 1/3, 残值≤max($10,1%V) 清光) > 2 退出解除(本月不买)
  > 3 q≥90 减仓(只卖不买) > 4 带再平衡([0.8w,1.2w]×V) > 5/6/7 补仓 S1/S2/S3(日频)
执行价 = 信号次一交易日 LBMA PM；成本单边 0.05%；杠杆上限 100%；风险日不买入。
"""
import bisect
import csv
import json
import math
from calendar import monthrange
from datetime import date, timedelta
from pathlib import Path

ROOT = Path(__file__).resolve().parents[1]
DATA = ROOT / "data"
HERE = Path(__file__).resolve().parent

COST = 0.0005          # 单边 0.05%
EXIT_LINE = 102.0      # 央行购金退出线（吨/季），预注册冻结


# ---------------------------------------------------------------- 数据加载

def load_pm():
    """LBMA PM（USD/oz）：[(iso_date, px)]，升序。"""
    rows = json.loads((DATA / "gold-lbma-pm.json").read_text(encoding="utf-8"))
    out = []
    for r in rows:
        v = r.get("v") or []
        if v and v[0]:
            out.append((r["d"], float(v[0])))
    out.sort(key=lambda x: x[0])
    return out


def load_daily_gold():
    lab = json.loads((DATA / "lab-data.json").read_text(encoding="utf-8"))
    return lab["daily"]["dates"], lab["daily"]["gold"]


def rv_sq_map(dates, gold):
    """复刻 model.py:86-95：月内日对数收益平方和（收益计入后一日所在月）；≥15 个收益日的月有效。"""
    acc, cnt = {}, {}
    prev = None
    for d, p in zip(dates, gold):
        if prev is not None:
            mm = d[:7]
            r = math.log(p) - math.log(prev)
            acc[mm] = acc.get(mm, 0.0) + r * r
            cnt[mm] = cnt.get(mm, 0) + 1
        prev = p
    return {k: acc[k] for k in acc if cnt[k] >= 15}


def sigma_series(month_axis, rv):
    """复刻 model.py sigma_rv:158-171：σ(月 t) 只用 ≤ t-1 月，权重 0.5^(age/6)，age>240 截断。"""
    out = {}
    for t, ym in enumerate(month_axis):
        sw = swv = 0.0
        i = t - 1
        while i >= 0:
            r2 = rv.get(month_axis[i])
            if r2 is not None:
                age = t - 1 - i
                w = 0.5 ** (age / 6)
                sw += w
                swv += w * r2
                if age > 240:
                    break
            i -= 1
        if sw > 0:
            out[ym] = math.sqrt(swv / sw)
    return out


def load_q_rows(path=None):
    """exp-03 W3 walk-forward 序列 → [(anchor_date, p, q_prior_pct)]。q=在**此前**已预测锚点中的分位。"""
    p = Path(path) if path else ROOT / "exp" / "results" / "exp03_regime.json"
    d = json.loads(p.read_text(encoding="utf-8"))
    rows = sorted(d["results"]["W3"]["rows"], key=lambda r: r["d"])
    out, past = [], []
    for r in rows:
        q = None
        if past:
            q = 100.0 * bisect.bisect_right(past, r["p"]) / len(past)
        out.append((r["d"], r["p"], q))
        bisect.insort(past, r["p"])
    return out


def load_cb():
    """WGC 季度央行购金：[(qlabel, quarter_end_date, tonnes)]，按季度末升序。"""
    lines = (DATA / "wgc" / "gdt-quarterly.csv").read_text(encoding="utf-8-sig").strip().splitlines()
    hdr = lines[0].split(",")
    out = []
    for ln in lines[1:]:
        r = dict(zip(hdr, ln.split(",")))
        qlab = r["quarter"].strip()
        y, q = int(qlab[:4]), int(qlab[5:6])
        nxt = date(y + (1 if q == 4 else 0), 1 if q == 4 else q * 3 + 1, 1)
        out.append((qlab, nxt - timedelta(days=1), float(r["cbTonnes"])))
    out.sort(key=lambda x: x[1])
    return out


def pace_at(cb, d):
    """t 日可见 pace：end+30 天 ≤ d 的最近季度值。"""
    best = None
    for _, end, t in cb:
        if end + timedelta(days=30) <= d:
            best = t
        else:
            break
    return best


def pace_last_two(cb, d):
    vals = [t for _, end, t in cb if end + timedelta(days=30) <= d]
    return vals[-2:] if len(vals) >= 2 else []


def prev_month_last(ym):
    y, m = int(ym[:4]), int(ym[5:7])
    py, pm = (y - 1, 12) if m == 1 else (y, m - 1)
    return date(py, pm, monthrange(py, pm)[1]).isoformat()


def build_signals(q_file=None):
    """共用信号集（主回测与对照任务 §9 共用，防止实现漂移）：
    σ 映射（RV-EWMA6 月度递推）、D14 回撤深度代理分位、exp-03 walk-forward q 序列、WGC pace、PM 序列。"""
    dates_all, gold_all = load_daily_gold()
    rv = rv_sq_map(dates_all, gold_all)
    month_axis = sorted({d[:7] for d in dates_all})
    sig = sigma_series(month_axis, rv)
    month_last_idx = {}
    for i, d in enumerate(dates_all):
        month_last_idx[d[:7]] = i
    depth63 = {}
    for ym, i in month_last_idx.items():
        lo = max(0, i - 62)
        seg = gold_all[lo:i + 1]
        peak = seg[0]
        worst = 0.0
        for p in seg:
            if p > peak:
                peak = p
            dd = p / peak - 1.0
            if dd < worst:
                worst = dd
        depth63[ym] = -worst                      # 深度取正值，越深分位越高
    yms_sorted = sorted(depth63)
    dvals = [depth63[ym] for ym in yms_sorted]
    q_proxy = {}
    for t, ym in enumerate(yms_sorted):
        window = dvals[max(0, t - 119):t + 1]
        if len(window) < 24:
            continue
        q_proxy[ym] = 100.0 * sum(1 for v in window if v <= dvals[t]) / len(window)
    return {"dates_all": dates_all, "gold_all": gold_all, "month_axis": month_axis,
            "sig": sig, "q_proxy": q_proxy,
            "qrows": load_q_rows(q_file), "cb": load_cb(), "pm": load_pm()}


def month_q_eff(S, qrows, ym):
    """有效 q（D13/D14）：真实 exp-03 walk-forward 分位（≤上月末最新锚点）优先，否则回撤深度代理；均缺 → None。"""
    lim = prev_month_last(ym)
    q = None
    for ad, _p, qq in qrows:
        if ad <= lim:
            q = qq
        else:
            break
    return q if q is not None else S["q_proxy"].get(ym)


# ---------------------------------------------------------------- 引擎

def run(start="2013-01-01", end="2025-12-31", target_vol_ann=12.0, qrows=None,
        V0=10000.0, tag="run", overlay=None, init_state=None):
    """跑一段回测，写 results/<tag>/，返回 (summary, final_state)。

    overlay（v3 研究 D17，默认全关）：
      cash_yield  闲置现金年化收益率 %（日计息）
      cap_half    dict ym->bool：True 的月份目标仓位 ×0.5（F2 慢趋势/气候地板；该月暂停买入档）
      breaker     True 启用日频闪断：RV5 > 1.5×当月σ（日频口径）→ 目标 ×0.5，翻转次日卖出至目标，
                  期间暂停买入档，RV5 回落自动解除（F3）
      cb_line     dict ym->CB*：逻辑退出线替换固定 102（F5；缺月或 γ≈0 回落 102）
    """
    sigma_star = (target_vol_ann / 100.0) / math.sqrt(12.0)
    ov = overlay or {}
    cash_yield_daily = (ov.get("cash_yield", 0.0) or 0.0) / 100.0 / 252.0
    cys = ov.get("cash_yield_series")                      # D19：[(date, 年利率%)] 逐日真实利率
    cys_d = [x[0] for x in cys] if cys else None
    cap_half = ov.get("cap_half")
    breaker_on = bool(ov.get("breaker", False))
    cb_line = ov.get("cb_line")

    S = build_signals()
    dates_all, gold_all = S["dates_all"], S["gold_all"]
    sig = S["sig"]
    q_proxy = S["q_proxy"]
    qrows = qrows if qrows is not None else S["qrows"]
    cb = S["cb"]
    pm = S["pm"]

    days = [(d, px) for d, px in pm if start <= d <= end]
    if not days:
        assert init_state is not None, f"窗口内无交易日: {start}~{end}（数据截至 {pm[-1][0]}）"
        days = [(pm[-1][0], pm[-1][1])]                      # 锚定今天：数据末日单点（只出状态与价格地图）
    # F3 日频熔断的 RV5（日频已实现波动，近 5 日对数收益 RMS）
    rv5 = [None] * len(days)
    _r = []
    for i, (d, px) in enumerate(days):
        if i > 0:
            _r.append(math.log(px / days[i - 1][1]))
        if len(_r) > 5:
            _r.pop(0)
        if len(_r) == 5:
            rv5[i] = math.sqrt(sum(x * x for x in _r) / 5.0)
    month_start_dates = set()                  # 日历意义上每月首个交易日（全序列判定）
    _prev_ym = None
    for d, _px in pm:
        if d[:7] != _prev_ym:
            month_start_dates.add(d)
            _prev_ym = d[:7]

    def w_target(s, q):
        if not s:
            return None
        base = min(1.0, sigma_star / s)
        rf = 1.0
        if q is not None:
            if q >= 90:
                rf = 0.5
            elif q >= 75:
                rf = 0.7
            elif q >= 60:
                rf = 0.85
        return base * rf

    # ---- 状态 ----
    cash, units = V0, 0.0
    mode = "FLAT"                      # FLAT / IN / EXITING
    entry_price, r0 = None, 0.0
    cost_basis = 0.0                   # 当前持仓的买入成本（卖出按比例结转）
    used = {"S1": False, "S2": False, "S3": False}
    routine_order = None               # (exec_j, side, amount, reason, w_sig or None)
    pending_add = None                 # (exec_j, amount, reason, [stage,...])
    pending_forced = None              # (exec_j, "sell", amount, reason) 熔断减仓（F3）
    brk = False                        # 熔断状态
    brk_flip = False
    fees = 0.0
    trades, curve = [], []
    bh_units = V0 * (1 - COST) / days[0][1]
    exit_spans = []                    # [(start_d, end_d)] EXITING 区间（供画图）

    # 现实锚定（D22）：以用户真实持仓为初始状态，从窗口起点继续执行规则
    if init_state:
        cash = float(init_state.get("cash", 0.0) or 0.0)
        units = float(init_state.get("units", 0.0) or 0.0)
        cost_basis = float(init_state.get("cost_basis", 0.0) or 0.0)
        entry_price = float(init_state["entry_price"]) if init_state.get("entry_price") else None
        mode = "IN" if units > 1e-12 else "FLAT"
        r0 = cash if mode == "IN" else 0.0
        for k in used:
            used[k] = bool((init_state.get("stages_used") or {}).get(k, False))
        bh_units = (cash + units * days[0][1]) * (1 - COST) / days[0][1]

    def do_buy(j, d, px, amount, reason, w_sig=None, is_entry=False):
        nonlocal cash, units, fees, mode, entry_price, r0, cost_basis
        amount = min(amount, cash / (1 + COST))
        if amount <= 1e-9:
            return None
        fee = amount * COST
        units += amount / px
        cost_basis += amount
        cash -= amount + fee
        fees += fee
        if is_entry:
            mode = "IN"
            entry_price = px
            r0 = (cash + units * px) * (1 - w_sig) if w_sig else 0.0
            for k in used:
                used[k] = False
        trades.append({"date": d, "action": "BUY", "reason": reason,
                       "amount_usd": round(amount, 2), "price": px,
                       "V_after": round(cash + units * px, 2)})
        return amount

    def do_sell(j, d, px, amount, reason):
        nonlocal cash, units, fees, mode, cost_basis
        amount = min(amount, units * px)
        if amount <= 1e-9:
            return None
        cost_basis *= (1 - amount / (units * px))       # 卖出部分按比例结转成本
        fee = amount * COST
        units -= amount / px
        cash += amount - fee
        fees += fee
        trades.append({"date": d, "action": "SELL", "reason": reason,
                       "amount_usd": round(amount, 2), "price": px,
                       "V_after": round(cash + units * px, 2)})
        return amount

    for j, (d, px) in enumerate(days):
        ym = d[:7]
        s = sig.get(ym)
        q = month_q_eff(S, qrows, ym)                           # 有效 q：真实优先，否则回撤深度代理（D13/D14）
        w = w_target(s, q)
        dd = date.fromisoformat(d)
        pace = pace_at(cb, dd)

        # ---- 0) 覆盖信号（D17：F2 地板 / F3 熔断 / F5 退出线）----
        line = EXIT_LINE
        if cb_line:
            _cl = cb_line.get(ym)
            if _cl is not None:
                line = _cl
        half_m = bool(cap_half and cap_half.get(ym))
        brk_now = bool(breaker_on and s and rv5[j] is not None
                       and rv5[j] > 1.5 * s / math.sqrt(21.0))
        if brk_now and not brk:
            brk, brk_flip = True, True
        elif not brk_now and brk:
            brk = False
        w_day = w * (0.5 if (half_m or brk) else 1.0) if w is not None else None

        # ---- 1) 执行到期订单 ----
        sold_today = False
        if routine_order and routine_order[0] == j:
            _, side, amount, reason, w_sig = routine_order
            routine_order = None
            if side == "buy":
                got = do_buy(j, d, px, amount, reason, w_sig=w_sig, is_entry=(reason == "ENTRY"))
            else:
                got = do_sell(j, d, px, amount, reason)
                sold_today = got is not None
                if units * px <= 1e-9 or units * px <= max(10.0, 0.01 * (cash + units * px)):
                    if mode == "EXITING" and units > 0 and units * px <= max(10.0, 0.01 * (cash + units * px)):
                        got = do_sell(j, d, px, units * px, reason + "+final")
                        sold_today = sold_today or (got is not None)
                    if units * px <= 1e-6:
                        units = 0.0
                        if mode == "EXITING":
                            mode = "FLAT"
        if pending_forced and pending_forced[0] == j:
            _, _side, amount, reason = pending_forced
            pending_forced = None
            if not sold_today and mode == "IN":
                got = do_sell(j, d, px, amount, reason)
                sold_today = sold_today or (got is not None)
        if pending_add and pending_add[0] == j:
            _, amount, reason, stages = pending_add
            pending_add = None
            if not sold_today and mode == "IN":
                got = do_buy(j, d, px, amount, reason)
                if got is not None:
                    for st in stages:
                        used[st] = True
            elif sold_today:
                trades.append({"date": d, "action": "SKIP", "reason": reason,
                               "amount_usd": 0, "price": px,
                               "V_after": round(cash + units * px, 2)})

        # ---- 2) 月度例行信号（当月首个交易日收盘评估，次日执行）----
        if d in month_start_dates and j + 1 < len(days):
            band_skip = False
            if mode in ("IN", "EXITING"):
                if pace is not None and pace < line:
                    hv = units * px
                    routine_order = (j + 1, "sell", hv / 3.0, f"EXIT pace<{line:.0f}", None)
                    mode = "EXITING"
                else:
                    if mode == "EXITING":
                        mode = "IN"
                        band_skip = True          # 退出解除当月不买（计划书 D2）
                    if mode == "IN" and w_day is not None:
                        hv = units * px
                        V = cash + hv
                        frac = hv / V if V > 0 else 0.0
                        if q is not None and q >= 90:
                            routine_order = (j + 1, "sell",
                                             max(hv / 3.0, hv - 1.2 * w_day * V), "RISK q>=90", None)
                        elif not band_skip:
                            if frac < 0.8 * w_day:
                                routine_order = (j + 1, "buy", w_day * V - hv, "BAND buy", None)
                            elif frac > 1.2 * w_day:
                                routine_order = (j + 1, "sell", hv - w_day * V, "BAND sell", None)
            elif mode == "FLAT" and w_day is not None:
                last2 = pace_last_two(cb, dd)
                if len(last2) == 2 and min(last2) >= line:
                    routine_order = (j + 1, "buy", w_day * cash, "ENTRY", w_day)

        # ---- 2b) 熔断翻转：次日卖出至减半目标（F3）----
        if brk_flip:
            brk_flip = False
            if j + 1 < len(days) and mode == "IN" and w_day is not None:
                hv = units * px
                V = cash + hv
                tgt = w_day * V                       # w_day 已含 ×0.5
                if hv > tgt:
                    pending_forced = (j + 1, "sell", hv - tgt, "BREAKER")

        # ---- 3) 日频补仓触发（收盘评估，次日执行；地板/熔断月暂停买入档）----
        if (mode == "IN" and s and w_day is not None and entry_price is not None
                and j + 1 < len(days) and not brk and not half_m):
            q_ok = (q is None) or (q < 90)
            trig = []
            if q_ok and not used["S1"] and px <= entry_price * (1 - s):
                trig.append((0.3 * r0, "ADD_S1", "S1"))
            if q_ok and not used["S2"] and px <= entry_price * (1 - 2 * s):
                trig.append((0.3 * r0, "ADD_S2", "S2"))
            s3_cond = s <= 0.06 or q is None or (q is not None and q <= 60)
            if s3_cond and not used["S3"]:
                hv = units * px
                V = cash + hv
                need = w_day * V - hv
                if need > 1e-9:
                    trig.append((min(need, cash), "ADD_S3", "S3"))
            if trig:
                pending_add = (j + 1, sum(a for a, _, _ in trig),
                               "+".join(rs for _, rs, _ in trig),
                               [st for _, _, st in trig])

        # ---- 4) 计息（当日利息计入当日净值）+ 记录 ----
        if cys:
            k = bisect.bisect_right(cys_d, d) - 1
            cash *= (1.0 + (cys[k][1] if k >= 0 else 0.0) / 100.0 / 252.0)
        else:
            cash *= (1.0 + cash_yield_daily)                # F1：闲置现金日计息
        hv = units * px
        V = cash + hv
        curve.append({"date": d, "V": V, "bh_V": bh_units * px, "frac": hv / V if V else 0.0,
                      "mode": mode, "sigma": s, "q": q, "pace": pace})
        if mode == "EXITING" and (not exit_spans or exit_spans[-1][1] is not None):
            exit_spans.append([d, None])
        elif mode != "EXITING" and exit_spans and exit_spans[-1][1] is None:
            exit_spans[-1][1] = d

    end_d = days[-1][0]
    for sp in exit_spans:
        if sp[1] is None:
            sp[1] = end_d

    final = curve[-1]
    state = {
        "as_of": end_d, "mode": mode, "V": final["V"], "cash": cash,
        "units": units, "px_last": days[-1][1], "cost_basis": cost_basis,
        "rv5_last": rv5[-1] if rv5 else None,
        "holdings_usd": units * days[-1][1], "frac": final["frac"],
        "entry_price": entry_price, "reserve_r0": r0,
        "stages_used": dict(used), "sigma_last": final["sigma"], "q_last": final["q"],
        "pace_visible": final["pace"], "fees_total": fees,
        "w_last": w, "s_last": s, "q_eff_last": q,
        "pending_routine": routine_order, "pending_add": pending_add,
    }
    summary = {"tag": tag, "start": days[0][0], "end": end_d,
               "target_vol_ann": target_vol_ann, "V0": V0, "V": final["V"],
               "fees": fees, "n_trades": sum(1 for t in trades if t["action"] != "SKIP"),
               "exit_spans": exit_spans}
    return summary, state, {"curve": curve, "trades": trades}


# ---------------------------------------------------------------- 输出

def write_outputs(out_dir, curve, trades):
    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)
    with open(out / "curve.csv", "w", newline="", encoding="utf-8") as f:
        wtr = csv.DictWriter(f, fieldnames=list(curve[0].keys()))
        wtr.writeheader()
        wtr.writerows(curve)
    with open(out / "trades.csv", "w", newline="", encoding="utf-8") as f:
        wtr = csv.DictWriter(f, fieldnames=["date", "action", "reason", "amount_usd",
                                            "price", "V_after"])
        wtr.writeheader()
        wtr.writerows(trades)