PythagorasGoal/scripts/research/xsec/runs/XD03.py

"""XD03 — Coskewness with Market

Mechanism: For each asset, compute rolling coskewness of asset returns
with the equal-weight market return. Assets with LOW coskewness (they do
not co-skew positively with the market) tend to earn a premium because
investors disfavor assets with negative coskewness (they hurt in crashes
when skewness matters most). Classic Harvey & Siddique (2000) anomaly.

Coskew(i, M) = E[(r_i - mu_i)(r_M - mu_M)^2] / (sigma_i * sigma_M^2)
Causally computed. LOWER coskew = LONG signal.

Grid: 5 backtests varying (win, H, k, universe, long_short).
"""

import sys
sys.path.insert(0, "/opt/docker/PythagorasGoal/scripts/research/xsec")
import xslib as xs
import numpy as np
import pandas as pd


def coskew_score(ret: np.ndarray, win: int = 60) -> np.ndarray:
    """Rolling coskewness of each asset with the equal-weight market.

    coskew(i, M) = E[(r_i - mu_i)(r_M - mu_M)^2] / (std_i * std_M^2)

    Returns (n_days x n_assets). LOWER = should be LONG (earns premium).
    So for long-low, negate: score = -coskew
    """
    n, A = ret.shape
    mkt = xs.market_ret(ret)  # (n,)

    out = np.full((n, A), np.nan)

    # Use pandas rolling for causality
    mkt_s = pd.Series(mkt)

    for a in range(A):
        asset_s = pd.Series(ret[:, a])

        # Rolling window stats
        mu_a = asset_s.rolling(win, min_periods=max(10, win // 3)).mean()
        mu_m = mkt_s.rolling(win, min_periods=max(10, win // 3)).mean()
        std_a = asset_s.rolling(win, min_periods=max(10, win // 3)).std()
        std_m = mkt_s.rolling(win, min_periods=max(10, win // 3)).std()

        # Centered series (element-wise)
        da = asset_s - mu_a
        dm = mkt_s - mu_m

        # coskew numerator = mean(da * dm^2)
        coskew_num = (da * dm ** 2).rolling(win, min_periods=max(10, win // 3)).mean()

        # Normalize by std_a * std_m^2
        denom = std_a * std_m ** 2
        denom = denom.replace(0, np.nan)

        coskew = coskew_num / denom
        out[:, a] = coskew.values

    return out


def score_fn_60(P):
    """Long low-coskew: negate so that lower coskew = higher score."""
    return -coskew_score(P.ret, win=60)


def score_fn_90(P):
    """Longer lookback for coskewness."""
    return -coskew_score(P.ret, win=90)


def score_fn_30(P):
    """Shorter lookback — more reactive."""
    return -coskew_score(P.ret, win=30)


if __name__ == "__main__":
    print("=== XD03: Coskewness with Market ===\n")

    results = []

    # Run 1: baseline config (win=60, all, H=10, k=5, LS)
    print("Run 1/5: win=60, universe=all, H=10, k=5, long_short=True")
    r1 = xs.study_xs("XD03-w60-H10-k5-LS", score_fn_60,
                     universe="all", H=10, k=5, long_short=True)
    print(xs.fmt(r1))
    print("JSON:", xs.as_json(r1))
    results.append(r1)

    # Run 2: vary rebalance period (H=20, looser)
    print("\nRun 2/5: win=60, universe=all, H=20, k=5, long_short=True")
    r2 = xs.study_xs("XD03-w60-H20-k5-LS", score_fn_60,
                     universe="all", H=20, k=5, long_short=True)
    print(xs.fmt(r2))
    print("JSON:", xs.as_json(r2))
    results.append(r2)

    # Run 3: longer win=90 (more stable coskewness estimate)
    print("\nRun 3/5: win=90, universe=all, H=10, k=5, long_short=True")
    r3 = xs.study_xs("XD03-w90-H10-k5-LS", score_fn_90,
                     universe="all", H=10, k=5, long_short=True)
    print(xs.fmt(r3))
    print("JSON:", xs.as_json(r3))
    results.append(r3)

    # Run 4: majors only (19 assets, cleaner signal)
    print("\nRun 4/5: win=60, universe=majors, H=10, k=5, long_short=True")
    r4 = xs.study_xs("XD03-w60-H10-k5-LS-maj", score_fn_60,
                     universe="majors", H=10, k=5, long_short=True)
    print(xs.fmt(r4))
    print("JSON:", xs.as_json(r4))
    results.append(r4)

    # Run 5: long-only on majors (captures risk-premium differently)
    print("\nRun 5/5: win=60, universe=majors, H=10, k=5, long_only")
    r5 = xs.study_xs("XD03-w60-H10-k5-LO-maj", score_fn_60,
                     universe="majors", H=10, k=5, long_short=False)
    print(xs.fmt(r5))
    print("JSON:", xs.as_json(r5))
    results.append(r5)

    # Summary: pick best by (earns_slot, then hold-out sharpe, then full sharpe)
    def rank_key(r):
        es = 1 if r["earns_slot"] else 0
        hs = r["holdout"].get("sharpe", -99)
        fs = r["full"]["sharpe"]
        corr_ok = (r.get("corr_xs01") or 1.0) < 0.6
        return (es, int(corr_ok), hs, fs)

    best = max(results, key=rank_key)

    print("\n\n=== BEST CONFIG ===")
    print(xs.fmt(best))
    print("JSON:", xs.as_json(best))