Blame for scripts/run_sentiment_replication.py

"""Replication of arXiv 2608.26174 -- sentiment + on-chain + market, daily BTC.

That paper is the strongest sentiment-inclusive result I could find that also
evaluates honestly: strict expanding walk-forward, a stated baseline set, and
transaction costs in its profit function. It reports ROC-AUC 0.6316 for its
"ImprovedTCN", against TCN 0.5692, ImprovedTCN_GRU 0.5314, XGBoost 0.5011,
LSTM 0.4586 and Random Forest 0.3627.

The reason its AUC beats anything achievable on plain hourly direction is its
target: not "will price go up" but "will price rise more than 5% within the next
7 days" -- a large-move question, with a positive class near a quarter of days.
Magnitude is forecastable in a way that sign is not, which is the same reason
the volatility-regime model here works while the direction models do not.

Data, and where it differs from the paper:

  market      Binance BTCUSDT hourly, resampled to daily        (as paper)
  sentiment   Crypto Fear & Greed Index, alternative.me         (as paper)
  on-chain    CoinMetrics community: MVRV, hashrate, active
              addresses, tx count, supply, market cap. Realized
              cap is derived as marketcap / MVRV and NUPL as
              1 - 1/MVRV, both exact identities.                (as paper)
  MISSING     Google Trends. The paper uses it; pytrends is rate
              limited and unreproducible, so it is omitted and
              the gap is reported rather than papered over.

Protocol follows the paper exactly: TimeSeriesSplit(n_splits=5, test_size=100),
expanding training window, no look-ahead. Every rolling statistic is causal.

The ablation is the part the paper does not run and the thesis needs: the same
models on market-only, market+on-chain, and market+on-chain+sentiment, so the
contribution of sentiment is measured rather than assumed.

Usage:
    python scripts/run_sentiment_replication.py
    python scripts/run_sentiment_replication.py --threshold 0.05 --horizon 7
"""
from __future__ import annotations

import argparse
import json
import sys
import time
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import requests
import torch
import torch.nn as nn
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import f1_score, precision_score, recall_score, roc_auc_score
from sklearn.model_selection import TimeSeriesSplit
from sklearn.preprocessing import StandardScaler

ROOT = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))
warnings.filterwarnings("ignore")

from lightgbm import LGBMClassifier                            # noqa: E402
from xgboost import XGBClassifier                              # noqa: E402
from scripts.run_direction_sota_experiment import load_candles  # noqa: E402

CACHE = ROOT / "outputs" / "sentiment_replication"
CM_METRICS = ["CapMVRVCur", "HashRate", "AdrActCnt", "TxCnt", "SplyCur", "CapMrktCurUSD"]
SEED = 42


# --------------------------------------------------------------------------- #
# data sources
# --------------------------------------------------------------------------- #
def fear_greed() -> pd.DataFrame:
    """Daily Crypto Fear & Greed Index, 2018-02-01 to now."""
    cache = CACHE / "fear_greed.csv"
    if cache.exists():
        return pd.read_csv(cache, parse_dates=["date"])
    payload = requests.get("https://api.alternative.me/fng/",
                           params={"limit": 0, "format": "json"}, timeout=60).json()["data"]
    df = pd.DataFrame({
        "date": pd.to_datetime([int(r["timestamp"]) for r in payload], unit="s", utc=True).normalize(),
        "fear_greed": [float(r["value"]) for r in payload],
    }).sort_values("date").drop_duplicates("date").reset_index(drop=True)
    cache.parent.mkdir(parents=True, exist_ok=True)
    df.to_csv(cache, index=False)
    return df


def coinmetrics() -> pd.DataFrame:
    """Free daily on-chain metrics, paged through the community endpoint."""
    cache = CACHE / "onchain.csv"
    if cache.exists():
        return pd.read_csv(cache, parse_dates=["date"])
    rows: list = []
    token = None
    while True:
        params = {"assets": "btc", "metrics": ",".join(CM_METRICS), "frequency": "1d",
                  "start_time": "2018-01-01", "page_size": 10000}
        if token:
            params["next_page_token"] = token
        body = requests.get("https://community-api.coinmetrics.io/v4/timeseries/asset-metrics",
                            params=params, timeout=60).json()
        rows += body.get("data", [])
        token = body.get("next_page_token")
        if not token:
            break
        time.sleep(0.2)
    df = pd.DataFrame(rows)
    df["date"] = pd.to_datetime(df["time"], utc=True).dt.normalize()
    for col in CM_METRICS:
        df[col] = pd.to_numeric(df.get(col), errors="coerce")
    df = df[["date"] + CM_METRICS].sort_values("date").drop_duplicates("date").reset_index(drop=True)
    cache.parent.mkdir(parents=True, exist_ok=True)
    df.to_csv(cache, index=False)
    return df


def daily_market() -> pd.DataFrame:
    candles = load_candles("2018-01-01", ROOT / "outputs" / "volatility" / "btc_1h.csv")
    agg = {"open": "first", "high": "max", "low": "min", "close": "last",
           "volume": "sum", "quote_volume": "sum", "num_trades": "sum"}
    have = {k: v for k, v in agg.items() if k in candles.columns}
    daily = candles.set_index("timestamp").resample("1D").agg(have).dropna().reset_index()
    daily["date"] = daily["timestamp"].dt.normalize()
    return daily.drop(columns=["timestamp"])


# --------------------------------------------------------------------------- #
# features -- every window is causal, nothing uses a future bar
# --------------------------------------------------------------------------- #
def build(threshold: float, horizon: int) -> tuple[pd.DataFrame, dict]:
    df = (daily_market()
          .merge(coinmetrics(), on="date", how="left")
          .merge(fear_greed(), on="date", how="left")
          .sort_values("date").reset_index(drop=True))
    # On-chain and sentiment publish once per day; forward-fill a late print
    # rather than dropping the row, but never backfill, which would look ahead.
    df[CM_METRICS + ["fear_greed"]] = df[CM_METRICS + ["fear_greed"]].ffill()

    eps = 1e-12
    x = df.copy()
    x["log_return"] = np.log(x["close"].clip(lower=eps)).diff()
    for w in (7, 14, 30):
        x[f"ret_{w}d"] = np.log(x["close"].clip(lower=eps)).diff(w)
        x[f"vol_{w}d"] = x["log_return"].rolling(w).std()
        x[f"close_sma{w}_log"] = np.log(x["close"].clip(lower=eps)
                                        / x["close"].rolling(w).mean().clip(lower=eps))
    x["range_pct"] = (x["high"] - x["low"]) / x["close"].clip(lower=eps)
    x["volume_logchg"] = np.log1p(x["volume"].clip(lower=0)).diff()
    x["trades_logchg"] = np.log1p(x["num_trades"].clip(lower=0)).diff()
    delta = x["close"].diff()
    up = delta.clip(lower=0).rolling(14).mean()
    down = (-delta.clip(upper=0)).rolling(14).mean()
    x["rsi_14"] = 100 - 100 / (1 + up / down.replace(0, np.nan))

    # on-chain, as the paper uses them: valuation levels plus their momentum
    x["mvrv"] = x["CapMVRVCur"]
    x["nupl"] = 1.0 - 1.0 / x["CapMVRVCur"].replace(0, np.nan)          # identity
    x["realised_cap"] = x["CapMrktCurUSD"] / x["CapMVRVCur"].replace(0, np.nan)
    x["mvrv_chg_7d"] = x["mvrv"].diff(7)
    x["hashrate_logchg_7d"] = np.log1p(x["HashRate"].clip(lower=0)).diff(7)
    x["active_addr_logchg_7d"] = np.log1p(x["AdrActCnt"].clip(lower=0)).diff(7)
    x["txcnt_logchg_7d"] = np.log1p(x["TxCnt"].clip(lower=0)).diff(7)
    x["mcap_real_log"] = np.log(x["CapMrktCurUSD"].clip(lower=eps)
                                / x["realised_cap"].clip(lower=eps))

    # sentiment, levels and dynamics, with a causal z-score only
    x["fg"] = x["fear_greed"]
    x["fg_chg_1d"] = x["fg"].diff()
    x["fg_chg_7d"] = x["fg"].diff(7)
    x["fg_sma7"] = x["fg"].rolling(7).mean()
    x["fg_extreme_fear"] = (x["fg"] <= 25).astype(float)
    x["fg_extreme_greed"] = (x["fg"] >= 75).astype(float)
    roll = x["fg"].rolling(90)
    x["fg_z90"] = (x["fg"] - roll.mean()) / roll.std().replace(0, np.nan)

    # target: the paper's definition -- does price rise more than `threshold`
    # at any point... it uses the forward return, so that is what is used here
    fwd = x["close"].shift(-horizon) / x["close"] - 1.0
    x["target"] = (fwd > threshold).astype(float)
    x.loc[fwd.isna(), "target"] = np.nan

    groups = {
        "market": ["log_return", "range_pct", "volume_logchg", "trades_logchg", "rsi_14"]
                  + [f"{p}_{w}d" for w in (7, 14, 30) for p in ("ret", "vol")]
                  + [f"close_sma{w}_log" for w in (7, 14, 30)],
        "onchain": ["mvrv", "nupl", "mvrv_chg_7d", "hashrate_logchg_7d",
                    "active_addr_logchg_7d", "txcnt_logchg_7d", "mcap_real_log"],
        "sentiment": ["fg", "fg_chg_1d", "fg_chg_7d", "fg_sma7", "fg_z90",
                      "fg_extreme_fear", "fg_extreme_greed"],
    }
    return x, groups


# --------------------------------------------------------------------------- #
# the paper's best model family: dilated causal TCN
# --------------------------------------------------------------------------- #
class TCN(nn.Module):
    """Dilated causal convolutions over a window of days, as in the paper."""

    def __init__(self, n_features: int, channels: int = 32, levels: int = 3, dropout: float = 0.2):
        super().__init__()
        layers: list[nn.Module] = []
        in_ch = n_features
        for level in range(levels):
            dilation = 2 ** level
            layers += [
                nn.Conv1d(in_ch, channels, kernel_size=3, padding=2 * dilation, dilation=dilation),
                nn.ReLU(), nn.BatchNorm1d(channels), nn.Dropout(dropout)]
            in_ch = channels
        self.body = nn.Sequential(*layers)
        self.head = nn.Linear(channels, 1)
        self.levels = levels

    def forward(self, x):                      # x: (batch, time, features)
        h = self.body(x.transpose(1, 2))
        h = h[:, :, : x.size(1)]               # causal crop
        return self.head(h[:, :, -1]).squeeze(-1)


def windows(values: np.ndarray, length: int) -> np.ndarray:
    """Sliding windows ending at each row, left-padded by edge repetition."""
    pad = np.repeat(values[:1], length - 1, axis=0)
    padded = np.vstack([pad, values])
    return np.stack([padded[i: i + length] for i in range(len(values))])


def fit_tcn(X_tr, y_tr, X_te, seq_len: int = 14, epochs: int = 60) -> np.ndarray:
    torch.manual_seed(SEED)
    model = TCN(X_tr.shape[1])
    opt = torch.optim.Adam(model.parameters(), lr=1e-3, weight_decay=1e-4)
    pos = max(y_tr.sum(), 1)
    loss_fn = nn.BCEWithLogitsLoss(
        pos_weight=torch.tensor([(len(y_tr) - pos) / pos], dtype=torch.float32))
    Xw = torch.tensor(windows(X_tr, seq_len), dtype=torch.float32)
    yt = torch.tensor(y_tr, dtype=torch.float32)
    model.train()
    for _ in range(epochs):
        opt.zero_grad()
        loss = loss_fn(model(Xw), yt)
        loss.backward()
        opt.step()
    model.eval()
    with torch.no_grad():
        return torch.sigmoid(model(torch.tensor(windows(X_te, seq_len),
                                                dtype=torch.float32))).numpy().ravel()


# --------------------------------------------------------------------------- #
def models():
    return {
        "LightGBM": LGBMClassifier(n_estimators=300, learning_rate=0.03, num_leaves=15,
                                   min_child_samples=30, subsample=0.8, subsample_freq=1,
                                   colsample_bytree=0.7, reg_lambda=1.0,
                                   random_state=SEED, verbose=-1),
        "XGBoost": XGBClassifier(n_estimators=300, learning_rate=0.03, max_depth=4,
                                 subsample=0.8, colsample_bytree=0.7, reg_lambda=1.0,
                                 random_state=SEED, eval_metric="logloss"),
        "RandomForest": RandomForestClassifier(n_estimators=400, min_samples_leaf=10,
                                               random_state=SEED, n_jobs=-1),
        "TCN": "tcn",
    }


def run(frame: pd.DataFrame, feats: list[str], folds: int, test_size: int,
        paper_threshold: float) -> list[dict]:
    ok = frame[feats + ["target"]].notna().all(axis=1) & np.isfinite(
        frame[feats].to_numpy(dtype=float)).all(axis=1)
    data = frame[ok].reset_index(drop=True)
    X = data[feats].to_numpy(dtype=float)
    y = data["target"].to_numpy().astype(int)

    out = []
    for fold, (tr, te) in enumerate(TimeSeriesSplit(n_splits=folds, test_size=test_size)
                                    .split(X), start=1):
        scaler = StandardScaler().fit(X[tr])
        X_tr, X_te = scaler.transform(X[tr]), scaler.transform(X[te])
        if len(set(y[tr])) < 2 or len(set(y[te])) < 2:
            continue
        for name, model in models().items():
            p = (fit_tcn(X_tr, y[tr], X_te) if model == "tcn"
                 else model.fit(X_tr, y[tr]).predict_proba(X_te)[:, 1])
            for tag, thr in (("0.5", 0.5), ("paper", paper_threshold)):
                pred = (p >= thr).astype(int)
                out.append(dict(
                    fold=fold, model=name, thr_tag=tag, n_test=len(te),
                    start=str(data["date"].iloc[te[0]])[:10], end=str(data["date"].iloc[te[-1]])[:10],
                    pos_rate=float(y[te].mean()),
                    auc=roc_auc_score(y[te], p),
                    precision=precision_score(y[te], pred, zero_division=0),
                    recall=recall_score(y[te], pred, zero_division=0),
                    f1=f1_score(y[te], pred, zero_division=0),
                    pred_rate=float(pred.mean())))
    return out


PAPER = {"ImprovedTCN": 0.6316, "TCN": 0.5692, "ImprovedTCN_GRU": 0.5314,
         "XGBoost": 0.5011, "LSTM": 0.4586, "RandomForest": 0.3627}


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--threshold", type=float, default=0.05, help="forward return threshold")
    ap.add_argument("--horizon", type=int, default=7, help="forward days")
    ap.add_argument("--folds", type=int, default=5)
    ap.add_argument("--test-size", type=int, default=100)
    ap.add_argument("--paper-threshold", type=float, default=0.11)
    args = ap.parse_args()

    frame, groups = build(args.threshold, args.horizon)
    sets = {
        "market only": groups["market"],
        "market + onchain": groups["market"] + groups["onchain"],
        "market + onchain + SENTIMENT": groups["market"] + groups["onchain"] + groups["sentiment"],
    }
    valid = frame["target"].notna()
    print(f"daily rows {int(valid.sum()):,}  "
          f"{str(frame.loc[valid, 'date'].min())[:10]} to {str(frame.loc[valid, 'date'].max())[:10]}")
    print(f"target: forward {args.horizon}d return > {args.threshold:.0%}   "
          f"positive rate {frame.loc[valid, 'target'].mean():.1%}  (paper: 26%)")
    print(f"fear & greed coverage: {frame['fg'].notna().mean():.1%} of rows")
    print(f"protocol: TimeSeriesSplit(n_splits={args.folds}, test_size={args.test_size}) "
          f"= {args.folds * args.test_size} test days\n")

    rows = []
    for label, feats in sets.items():
        print(f"running [{label}]  {len(feats)} features ...")
        for r in run(frame, feats, args.folds, args.test_size, args.paper_threshold):
            r["feature_set"] = label
            rows.append(r)

    df = pd.DataFrame(rows)
    print("\n" + "=" * 100)
    print(f"ROC-AUC by feature set and model (mean of {args.folds} folds)")
    print("  paper's reported AUCs: " + ", ".join(f"{k} {v:.4f}" for k, v in PAPER.items()))
    print("=" * 100)
    auc = (df[df.thr_tag == "0.5"].groupby(["feature_set", "model"])["auc"]
           .agg(["mean", "std", "count"]).reset_index())
    print(f"{'feature set':<32}{'model':<14}{'AUC':>8}{'sd':>8}{'vs paper TCN':>14}")
    for _, r in auc.iterrows():
        ref = PAPER.get(r["model"], PAPER["ImprovedTCN"])
        print(f"{r['feature_set']:<32}{r['model']:<14}{r['mean']:>8.4f}{r['std']:>8.4f}"
              f"{r['mean'] - ref:>+14.4f}")

    print("\n" + "-" * 100)
    print("SENTIMENT ABLATION -- does adding Fear & Greed change the AUC?")
    print("-" * 100)
    base = auc[auc.feature_set == "market + onchain"].set_index("model")["mean"]
    full = auc[auc.feature_set == "market + onchain + SENTIMENT"].set_index("model")["mean"]
    for model in full.index:
        if model in base.index:
            d = full[model] - base[model]
            print(f"  {model:<14} {base[model]:.4f} -> {full[model]:.4f}   {d:+.4f}"
                  f"   {'helps' if d > 0.01 else 'hurts' if d < -0.01 else 'no effect'}")

    print("\n" + "-" * 100)
    print(f"At the paper's operating threshold ({args.paper_threshold}) -- "
          "paper reports precision 19.3%, recall 100.0%, F1 0.3234")
    print("-" * 100)
    op = (df[(df.thr_tag == "paper") & (df.feature_set == "market + onchain + SENTIMENT")]
          .groupby("model")[["precision", "recall", "f1", "pred_rate", "pos_rate"]].mean())
    print(f"{'model':<14}{'precision':>11}{'recall':>9}{'F1':>8}{'says yes':>10}{'base rate':>11}")
    for model, r in op.iterrows():
        print(f"{model:<14}{r['precision'] * 100:>10.1f}%{r['recall'] * 100:>8.1f}%"
              f"{r['f1']:>8.4f}{r['pred_rate'] * 100:>9.1f}%{r['pos_rate'] * 100:>10.1f}%")

    out = ROOT / "outputs" / "sentiment_replication.json"
    out.write_text(json.dumps(rows, indent=2, default=str), encoding="utf-8")
    print(f"\nwritten: {out}")


if __name__ == "__main__":
    main()