core/analysis.py

"""
analysis.py — NLG报告、行业影响、评估指标、消融实验
====================================================
从 v2_championship.py 拆出的分析与评估模块。
"""

import pandas as pd
import numpy as np
from sklearn.linear_model import QuantileRegressor
from sklearn.preprocessing import StandardScaler

from config import FACTOR_GROUPS, INDUSTRIES, INDUSTRY_ZH


# ═══════════════════════════════════════════════════════════
# INDUSTRY IMPACT RULE ENGINE
# ═══════════════════════════════════════════════════════════

def apply_industry_rules(rec):
    """基于风险等级、偏置、因子等推断各行业的风险和建议。"""
    risk_level = rec.get('risk_level', 'Low')
    risk_bias = rec.get('risk_bias', 'Balanced')
    vol_ratio = rec.get('vol_ratio', 1.0)
    top_fac = rec.get('top_factor', 'Unknown')

    is_high = risk_level == 'High'
    is_medium = risk_level == 'Medium'
    is_upward = risk_bias == 'Upward'
    is_downward = risk_bias == 'Downward'
    wide = vol_ratio > 1.3

    rules = {}
    for industry in INDUSTRIES:
        ind_risk = 'Low'
        ind_action = 'Routine monitoring'

        if industry == 'Aviation':
            if is_high and is_upward:
                ind_risk, ind_action = 'High', 'Increase hedging coverage; review fuel cost budget'
            elif is_high:
                ind_risk, ind_action = 'High', 'Elevated volatility; prepare contingency liquidity'
            elif is_upward:
                ind_risk, ind_action = 'Medium-High', 'Monitor fuel cost exposure; consider forward contracts'
            elif is_medium:
                ind_risk, ind_action = 'Medium', 'Review quarterly fuel hedging strategy'
            if top_fac == 'Demand': ind_action += '; demand-driven → cost pressure may persist'
            elif top_fac == 'Supply': ind_action += '; supply-driven → watch OPEC decisions'
            elif top_fac == 'Risk_Geo': ind_action += '; geopolitical risk → event monitoring'

        elif industry == 'Logistics':
            if is_high and is_upward:
                ind_risk, ind_action = 'Medium-High', 'Review transport cost pass-through; working capital buffer'
            elif is_high:
                ind_risk, ind_action = 'Medium', 'Monitor diesel/freight cost exposure'
            elif is_medium and is_upward:
                ind_risk, ind_action = 'Medium', 'Review fuel surcharge mechanisms'

        elif industry == 'Chemicals':
            if is_upward and top_fac == 'Supply':
                ind_risk, ind_action = 'High', 'Feedstock cost pressure; margin compression likely'
            elif is_high:
                ind_risk, ind_action = 'Medium-High', 'Monitor naphtha/ethylene spread; review procurement'
            elif wide:
                ind_risk, ind_action = 'Medium', 'Profit uncertainty elevated; scenario planning advised'

        elif industry == 'Manufacturing':
            if is_high and is_upward:
                ind_risk, ind_action = 'High', 'Energy cost surge risk; review energy hedging'
            elif is_high:
                ind_risk, ind_action = 'Medium', 'Elevated input cost volatility'
            elif is_medium:
                ind_risk, ind_action = 'Medium', 'Monitor energy procurement costs'

        elif industry == 'Upstream_OG':
            if is_downward and is_high:
                ind_risk, ind_action = 'High', 'Revenue decline risk; review covenant compliance & liquidity'
            elif is_downward:
                ind_risk, ind_action = 'Medium-High', 'Downside tail expanding; monitor cash flow coverage'
            elif is_upward:
                ind_risk, ind_action = 'Low', 'Revenue tailwind; capex commitment review advised'
            else:
                ind_risk, ind_action = 'Low-Medium', 'Balanced outlook'

        rules[f'{industry}_risk'] = ind_risk
        rules[f'{industry}_action'] = ind_action
    return rules


# ═══════════════════════════════════════════════════════════
# REGIME ECONOMIC NARRATIVES
# ═══════════════════════════════════════════════════════════

REGIME_NARRATIVES = {
    '2008 金融危机': {
        'history': '2008年全球金融危机期间，油价从$147暴跌至$32，降幅78%。需求端崩塌是主因——全球GDP收缩、贸易锐减、制造业PMI普遍跌破40。',
        'implication': '需求崩塌格局下，下游成本端企业短期受益于低油价，但总需求萎缩拖累整体营收。上游企业面临最大冲击。',
        'hedge_advice': '上游企业应立即锁定远期销售价格；下游企业可逢低建仓，锁定低价原料。',
    },
    '2014 页岩油冲击': {
        'history': '2014-16年美国页岩油产量爆发（+400万桶/日），叠加OPEC拒绝减产，油价从$110跌至$26。供给过剩主导。',
        'implication': '供给过剩格局持续时间长（18个月以上），下游企业成本优势可持续；上游企业需要重组债务、缩减资本开支。',
        'hedge_advice': '重点关注远期曲线结构（contango加深），利用期货锁价窗口。下游企业延长采购合约期限。',
    },
    '2020 COVID': {
        'history': '2020年COVID-19导致全球需求暴减2000万桶/日，WTI期货历史性跌至负值。需求冲击+仓储危机双重打击。',
        'implication': '极端需求冲击下，航空业客运量降90%+，物流链中断。但复苏速度可能超预期——V型反弹是历史常态。',
        'hedge_advice': '短期：保持现金流弹性，避免过度套保。中期：关注OPEC+协调减产信号，逢低建立多头头寸。',
    },
    '2022 俄乌冲突': {
        'history': '2022年俄乌冲突导致俄罗斯原油出口受制裁，供给缺口+地缘溢价推升布伦特至$130+。地缘政治+供给双重冲击。',
        'implication': '地缘驱动的价格飙升通常突然但短暂（3-6个月），随后制裁适应和替代供给逐步消化溢价。',
        'hedge_advice': '事件驱动行情中，期权策略优于期货——买入看涨期权锁定上限成本，保留价格回落的收益空间。',
    },
    '2023 OPEC减产': {
        'history': '2023年OPEC+主动减产200万桶/日，托底油价在$70-90区间。供给管理型市场，价格波动率较低。',
        'implication': 'OPEC减产格局下，价格区间可预测性较高，但下行风险来自减产执行率下滑和非OPEC增产。',
        'hedge_advice': '低波动环境适合使用零成本领（collar）策略，锁定窄价格带。',
    },
    '常态/低波动': {
        'history': '油价处于常态波动区间，无明显单一因子主导。市场处于供需基本平衡状态。',
        'implication': '常态下关注结构性变化信号——OPEC会议决策、美国钻井数趋势、中国PMI走向。',
        'hedge_advice': '常态下对冲比例可适当降低（25-40%），使用低成本期货锁价即可。',
    },
}


# ═══════════════════════════════════════════════════════════
# NLG REPORT GENERATION (ENHANCED)
# ═══════════════════════════════════════════════════════════

def generate_nlg_report(row):
    """生成单月深度风险研判报告，含经济学叙事和对冲建议。"""
    date = pd.Timestamp(row['test_date']).strftime('%Y年%m月')
    rl = row['risk_level']
    rb = row['risk_bias']
    top = row.get('top_factor', 'Unknown')
    q10 = row['pred_q10_1m'] * 100
    q50 = row['pred_q50_1m'] * 100
    q90 = row['pred_q90_1m'] * 100
    vol = row['pred_vol'] * 100

    rl_zh = {'Low': '低', 'Medium': '中等', 'High': '高'}.get(rl, rl)
    rb_zh = {'Upward': '偏上行', 'Downward': '偏下行', 'Balanced': '均衡'}.get(rb, rb)
    factor_zh = {
        'Price': '价格联动', 'Supply': '供给端', 'Demand': '需求端',
        'Risk_Geo': '地缘政治/风险', 'Technical': '技术面',
    }.get(top, top)

    # ── Part 1: 核心判断 ──
    summary = (
        f"【{date}油价风险研判】\n"
        f"■ 核心判断：风险等级{rl_zh}，方向{rb_zh}，由{factor_zh}因子主导。\n"
        f"■ 1M预测区间：[{q10:+.1f}%, {q90:+.1f}%]，中枢{q50:+.1f}%，波动率{vol:.1f}%。\n"
    )

    # 3M
    if pd.notna(row.get('pred_q10_3m')):
        summary += f"■ 3M预测区间：[{row['pred_q10_3m']*100:+.1f}%, {row['pred_q90_3m']*100:+.1f}%]，中枢{row['pred_q50_3m']*100:+.1f}%。\n"

    # CQR
    cqr_lo = row.get('cqr_q10_1m', None)
    if cqr_lo is not None and pd.notna(cqr_lo):
        summary += f"■ CQR校准区间：[{cqr_lo*100:+.1f}%, {row['cqr_q90_1m']*100:+.1f}%]（分布自由覆盖保证）。\n"

    # ── Part 2: Regime 经济学叙事 ──
    regime = row.get('regime_match', '')
    regime_sim = row.get('regime_similarity', 0)
    if regime and regime != 'Unknown':
        summary += f"\n▶ 格局识别：当前最接近「{regime}」（相似度{regime_sim*100:.0f}%）\n"
        narr = REGIME_NARRATIVES.get(regime, {})
        if narr:
            summary += f"  历史参照：{narr['history']}\n"
            summary += f"  当前启示：{narr['implication']}\n"
            summary += f"  对冲建议：{narr['hedge_advice']}\n"

    # ── Part 3: 行业影响 ──
    high_risk = []
    med_risk = []
    for ind in INDUSTRIES:
        risk = str(row.get(f'{ind}_risk', 'Low'))
        if 'High' in risk:
            high_risk.append(INDUSTRY_ZH.get(ind, ind))
        elif 'Medium' in risk:
            med_risk.append(INDUSTRY_ZH.get(ind, ind))
    if high_risk:
        summary += f"\n▶ 高风险行业：{'、'.join(high_risk)}——建议提升套保覆盖率至60-80%。\n"
    if med_risk:
        summary += f"▶ 中风险行业：{'、'.join(med_risk)}——建议维持25-50%套保覆盖。\n"
    if not high_risk and not med_risk:
        summary += f"\n▶ 各行业风险均处于可控水平，建议维持常规套保比例。\n"

    # ── Part 4: 压力测试 ──
    base = row.get('scenario_base', 0) * 100
    vix = row.get('scenario_vix_shock', 0) * 100
    supply = row.get('scenario_supply_cut', 0) * 100
    demand = row.get('scenario_demand_crash', 0) * 100
    worst = min(supply, demand)
    summary += (
        f"\n▶ 压力测试：基准{base:+.1f}% | VIX翻倍{vix:+.1f}% | "
        f"供给中断{supply:+.1f}% | 需求崩塌{demand:+.1f}%\n"
        f"  最大下行风险：{worst:+.1f}%，建议预留相应流动性缓冲。\n"
    )

    return summary


def generate_all_reports(results):
    """为所有月份生成 NLG 报告。"""
    reports = {}
    for _, row in results.iterrows():
        dt = pd.Timestamp(row['test_date']).strftime('%Y-%m')
        reports[dt] = generate_nlg_report(row)
    return reports


# ═══════════════════════════════════════════════════════════
# EVALUATION METRICS
# ═══════════════════════════════════════════════════════════

def evaluate_results(results):
    """计算全面的评估指标并打印。"""
    # 只评估有实际值的月份（排除 live forecast）
    eval_mask = results['actual_ret_1m'].notna()
    results = results[eval_mask].copy()
    ar = results['actual_ret_1m'].values
    av = results['actual_vol'].values
    n = len(results)

    print(f"\n{'='*65}")
    print("V2 CHAMPIONSHIP — EVALUATION")
    print("=" * 65)

    # 1M Interval
    print(f"\n--- 1M INTERVAL ---")
    models = [('QR (vol-adapt)', 'pred'), ('LightGBM', 'lgb')]
    if 'cqr_q10_1m' in results.columns:
        models.append(('Conformal QR', 'cqr'))
    for model_name, prefix in models:
        q10 = results[f'{prefix}_q10_1m'].values
        q90 = results[f'{prefix}_q90_1m'].values
        cov = ((ar >= q10) & (ar <= q90)).mean()
        wis = ((q90-q10) + (2/0.2)*np.maximum(q10-ar, 0) + (2/0.2)*np.maximum(ar-q90, 0)).mean()
        vm = np.median(av)
        hi = av > vm
        cov_hi = ((ar[hi] >= q10[hi]) & (ar[hi] <= q90[hi])).mean() if hi.sum() > 0 else 0
        print(f"  {model_name:<18} Cov={cov:.1%} HiCov={cov_hi:.1%} WIS={wis:.4f}")

    nq10 = np.full(n, np.quantile(ar, 0.10))
    nq90 = np.full(n, np.quantile(ar, 0.90))
    naive_wis = ((nq90-nq10) + (2/0.2)*np.maximum(nq10-ar, 0) + (2/0.2)*np.maximum(ar-nq90, 0)).mean()
    print(f"  {'Naive':<18} WIS={naive_wis:.4f}")

    # 3M
    print(f"\n--- 3M INTERVAL ---")
    m3 = results['actual_ret_3m'].notna()
    if m3.sum() > 10:
        ar3 = results.loc[m3, 'actual_ret_3m'].values
        q10_3 = results.loc[m3, 'pred_q10_3m'].values
        q90_3 = results.loc[m3, 'pred_q90_3m'].values
        cov3 = ((ar3 >= q10_3) & (ar3 <= q90_3)).mean()
        wis3 = ((q90_3-q10_3) + (2/0.2)*np.maximum(q10_3-ar3, 0) + (2/0.2)*np.maximum(ar3-q90_3, 0)).mean()
        print(f"  QR 3M (vol-adapt)  Cov={cov3:.1%} WIS={wis3:.4f} n={m3.sum()}")

    # Vol
    print(f"\n--- VOL SCORE ---")
    pv = results['pred_vol'].values
    bl = results['baseline_ewma'].values
    for nm, prd in [('V2 BL+Resid', pv), ('EWMA', bl)]:
        from sklearn.metrics import mean_squared_error
        rmse = np.sqrt(mean_squared_error(av, prd))
        corr = np.corrcoef(av, prd)[0, 1]
        print(f"  {nm:<18} RMSE={rmse:.4f} corr={corr:+.3f}")

    # Risk levels
    print(f"\n--- RISK LEVELS ---")
    for lvl in ['Low', 'Medium', 'High']:
        mask = results['risk_level'] == lvl
        if mask.sum() > 0:
            print(f"  {lvl:<8}: vol={av[mask].mean():.4f} n={mask.sum()}")

    # Factor frequency
    print(f"\n--- FACTOR FREQ ---")
    if 'top_factor' in results.columns:
        for fac, cnt in results['top_factor'].value_counts().items():
            print(f"  {fac:<12}: {cnt} ({cnt/n:.1%})")

    # SHAP
    print(f"\n--- SHAP (avg) ---")
    for g in FACTOR_GROUPS:
        col = f'shap_{g}'
        if col in results.columns:
            print(f"  {g:<12}: {results[col].abs().mean():.4f}")

    # Scenario
    print(f"\n--- SCENARIO (latest) ---")
    lat = results.iloc[-1]
    print(f"  Base:         {lat['scenario_base']*100:+.1f}%")
    print(f"  VIX shock:    {lat['scenario_vix_shock']*100:+.1f}%")
    print(f"  Supply cut:   {lat['scenario_supply_cut']*100:+.1f}%")
    print(f"  Demand crash: {lat['scenario_demand_crash']*100:+.1f}%")

    # NLG sample
    print(f"\n--- NLG REPORT (latest) ---")
    print(generate_nlg_report(results.iloc[-1]))


# ═══════════════════════════════════════════════════════════
# ABLATION EXPERIMENTS
# ═══════════════════════════════════════════════════════════

def _run_qr_eval(panel, feat_list, train_window=120):
    """内部辅助：对给定特征子集跑 walk-forward QR 并返回 (cov, wis, n)。"""
    hits, totals, wis_list = 0, 0, []
    for i in range(train_window, len(panel) - 1):
        train_df = panel.iloc[max(0, i - train_window):i]
        test_df = panel.iloc[i:i + 1]
        avail = [f for f in feat_list if f in train_df.columns and train_df[f].notna().mean() > 0.8]
        if len(avail) < 2:
            continue
        X_tr = train_df[avail].fillna(train_df[avail].median())
        X_te = test_df[avail].fillna(train_df[avail].median())
        sc = StandardScaler()
        X_tr_s = sc.fit_transform(X_tr)
        X_te_s = sc.transform(X_te)
        y = train_df['target_ret_1m'].dropna()
        mask = y.index.isin(X_tr.index)
        y = y[mask]
        X_tr_s = X_tr_s[:len(y)]
        actual = panel['target_ret_1m'].iloc[i]
        if np.isnan(actual):
            continue
        try:
            qr10 = QuantileRegressor(quantile=0.10, alpha=0.1, solver='highs')
            qr90 = QuantileRegressor(quantile=0.90, alpha=0.1, solver='highs')
            qr10.fit(X_tr_s, y)
            qr90.fit(X_tr_s, y)
            p10, p90 = qr10.predict(X_te_s)[0], qr90.predict(X_te_s)[0]
            if p10 > p90: p10, p90 = p90, p10
            hit = 1 if p10 <= actual <= p90 else 0
            hits += hit
            totals += 1
            wis_list.append((p90-p10) + (2/0.2)*max(p10-actual, 0) + (2/0.2)*max(actual-p90, 0))
        except:
            continue
    cov = hits / totals if totals > 0 else 0
    wis = np.mean(wis_list) if wis_list else 999
    return cov, wis, totals


def run_ablation(panel, features):
    """运行消融实验：训练窗口 + 因子组 leave-one-out。"""
    print(f"\n--- ABLATION EXPERIMENTS ---")
    ablation_results = []

    # ── Part 1: Window ablation ──
    print("  [窗口消融]")
    for w in [84, 120, 180]:
        cov, wis, n_ab = _run_qr_eval(panel, features, train_window=w)
        ablation_results.append({
            'type': 'window', 'param': w, 'param_label': f'{w}月',
            'cov': round(cov, 3), 'wis': round(wis, 4), 'n': n_ab
        })
        print(f"    Window={w:>3}: Cov={cov:.1%} WIS={wis:.4f} n={n_ab}")

    # ── Part 2: Factor-group leave-one-out ──
    print("  [因子组消融 — Leave-One-Out]")
    # Baseline: all factors
    base_cov, base_wis, base_n = _run_qr_eval(panel, features)
    ablation_results.append({
        'type': 'factor_group', 'param': 'ALL', 'param_label': '全部因子',
        'cov': round(base_cov, 3), 'wis': round(base_wis, 4), 'n': base_n
    })
    print(f"    ALL (baseline): Cov={base_cov:.1%} WIS={base_wis:.4f}")

    group_zh = {'Price': '价格联动', 'Supply': '供给端', 'Demand': '需求端',
                'Risk_Geo': '地缘/风险', 'Technical': '技术面'}

    for group, members in FACTOR_GROUPS.items():
        # Remove this group's features
        reduced = [f for f in features if f not in members]
        if len(reduced) < 2:
            continue
        cov, wis, n_ab = _run_qr_eval(panel, reduced)
        delta_cov = cov - base_cov
        delta_wis = wis - base_wis
        ablation_results.append({
            'type': 'factor_group', 'param': group,
            'param_label': f'去除{group_zh.get(group, group)}',
            'cov': round(cov, 3), 'wis': round(wis, 4), 'n': n_ab,
            'delta_cov': round(delta_cov, 3), 'delta_wis': round(delta_wis, 4),
        })
        print(f"    -{group:<12}: Cov={cov:.1%} (Δ{delta_cov:+.1%}) "
              f"WIS={wis:.4f} (Δ{delta_wis:+.4f})")

    return ablation_results