portfolio_analyzer/

research.rs

//! Research universe and shared helpers for the backtesters.

/// A deliberately mixed research universe: index and sector ETFs, mega caps,
/// cyclicals, defensives, laggards and a few high-volatility names.  It is
/// *not* the investable portfolio; it exists so the indicator is judged on
/// many kinds of price behaviour, including stocks that went nowhere.
pub const UNIVERSE: &[&str] = &[
    // Broad market, size, style, rates, commodities
    "SPY", "QQQ", "IWM", "DIA", "MDY", "EFA", "EEM", "TLT", "IEF", "GLD", "SLV", "USO", "UNG",
    // Sector ETFs
    "XLK", "XLF", "XLE", "XLV", "XLI", "XLY", "XLP", "XLU", "XLB", "XLRE", "SMH", "XBI", "KRE",
    // Technology / communication
    "AAPL", "MSFT", "NVDA", "GOOGL", "META", "AMZN", "ORCL", "CSCO", "INTC", "IBM", "ADBE", "CRM",
    "QCOM", "TXN", "AMD", "MU", "NFLX", "AVGO", "AMAT", "T", "VZ", "DIS",
    // Financials
    "JPM", "BAC", "WFC", "C", "GS", "MS", "AXP", "BLK", "SCHW", "USB", "PNC", "AIG", "MET",
    // Health care
    "JNJ", "PFE", "MRK", "UNH", "ABT", "LLY", "BMY", "AMGN", "GILD", "MDT", "CVS", "TMO",
    // Consumer
    "WMT", "PG", "KO", "PEP", "MCD", "NKE", "SBUX", "HD", "LOW", "TGT", "COST", "CL", "KHC", "F", "GM",
    // Industrials / materials / energy / utilities / REIT
    "CAT", "DE", "BA", "GE", "HON", "MMM", "UPS", "FDX", "LMT", "RTX", "UNP", "XOM", "CVX", "COP",
    "SLB", "OXY", "DOW", "NEM", "FCX", "NEE", "DUK", "SO", "AMT", "SPG", "O",
    // Higher-volatility / mid caps
    "TSLA", "SHOP", "UBER", "PYPL", "XYZ", "ROKU", "ETSY", "AAL", "CCL", "MRNA", "PLTR", "SNAP",
];

// ════════════════════════════════════════════════════════════════════════════
// The steady backtester
// ════════════════════════════════════════════════════════════════════════════
//
// "Steady" means three things here:
//
// 1. **Deterministic.**  Every random draw comes from a seeded generator, so
//    the same command prints the same numbers tomorrow.
// 2. **Honest about time.**  Models learn only from bars before
//    `TRAIN_END`, are tuned on bars before `VALID_END`, and are graded on
//    the bars after.  A gap of `h` bars is dropped at each boundary so that a
//    label never peeks across it.
// 3. **Shaken, not stirred.**  Gaussian noise is added to the *history the
//    indicator sees*, while outcomes are always measured on the *real*
//    future.  If SD only works on the exact historical prints, the noise
//    will expose it.

use crate::history::{self, DatedBar};
use crate::indicators::OHLCV;
use crate::sd2::classic::classic_series;
use crate::sd2::features::{FeatureConfig, FeatureFrame};
use chrono::NaiveDate;
use std::path::Path;

/// Forecast horizons in trading days: one week, two weeks, one month.
pub const HORIZONS: [usize; 4] = [5, 10, 21, 63];

/// Triple-barrier settings: which of ±`BARRIER_ATR`·ATR is touched first
/// within `BARRIER_BARS` bars.  This is the cleanest form of the question
/// "which way does the squeeze break?".
pub const BARRIER_ATR: f64 = 1.5;
pub const BARRIER_BARS: usize = 21;

pub fn train_end() -> NaiveDate {
    NaiveDate::from_ymd_opt(2014, 1, 1).unwrap()
}
pub fn valid_end() -> NaiveDate {
    NaiveDate::from_ymd_opt(2019, 1, 1).unwrap()
}

/// Which slice of history a sample belongs to.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum Split {
    Train,
    Valid,
    Test,
}

/// One (symbol, bar) observation.
#[derive(Debug, Clone)]
pub struct Sample {
    pub symbol: u16,
    pub date: NaiveDate,
    pub split: Split,
    pub x: Vec<f64>,
    /// Direction-free inputs for the movement model.
    pub mag_x: Vec<f64>,
    /// Real forward log return for each of [`HORIZONS`]; NaN if unavailable.
    pub fwd: [f64; 4],
    /// +1 if the upper barrier was touched first, -1 lower, 0 neither / both.
    pub barrier: f64,
    /// Forward log return of the symbol minus SPY's, per horizon.
    pub excess: [f64; 4],
    /// Daily volatility at the cutoff, so moves can be expressed in σ units.
    pub vol: f64,
    pub close: f64,
    pub classic_score: f64,
    pub alive: bool,
    pub comp: u8,
    pub fired_start: bool,
    /// τ/θ: how far through its expected life the current fire is (NaN if not fired).
    pub phase: f64,
    pub energy: f64,
    pub tier21: u8,
    pub tier10: u8,
    pub bars_in_squeeze: usize,
}

/// A loaded symbol: name plus dated bars.
pub struct Series {
    pub symbol: String,
    pub bars: Vec<DatedBar>,
}

pub fn load_universe(dir: &Path) -> Vec<Series> {
    let mut out = Vec::new();
    for s in UNIVERSE {
        if let Ok(bars) = history::load(&dir.join(format!("{s}.csv"))) {
            if bars.len() > 600 {
                out.push(Series { symbol: s.to_string(), bars });
            }
        }
    }
    out
}

/// Small, fast, seedable PRNG (SplitMix64) with Gaussian draws.
pub struct Rng(u64);
impl Rng {
    pub fn new(seed: u64) -> Self {
        Self(seed ^ 0x9E37_79B9_7F4A_7C15)
    }
    pub fn next_u64(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }
    pub fn uniform(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }
    /// Standard normal via Box–Muller.
    pub fn gauss(&mut self) -> f64 {
        let u1 = self.uniform().max(1e-300);
        let u2 = self.uniform();
        (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos()
    }
}

/// Add Gaussian *measurement* noise to a price history.
///
/// Each bar is multiplied by `exp(ε)`, `ε ~ N(0, (level · σ)²)`, where `σ` is
/// that symbol's typical daily log-return volatility.  `level = 0.5` means
/// "every print is wrong by about half a normal day's move", which is a lot
/// more than any real data error.  On top of that, high and low each get a
/// smaller independent wobble so the bar's *shape* is perturbed too (that is
/// what ATR and the Keltner channels feel).
pub fn add_noise(bars: &[OHLCV], level: f64, rng: &mut Rng) -> Vec<OHLCV> {
    if level <= 0.0 {
        return bars.to_vec();
    }
    let mut rets: Vec<f64> = bars.windows(2).map(|w| (w[1].close / w[0].close).ln().abs()).collect();
    rets.sort_by(|a, b| a.partial_cmp(b).unwrap());
    // Median absolute return × 1.4826 is a robust σ that ignores crash days.
    let sigma = rets.get(rets.len() / 2).copied().unwrap_or(0.01) * 1.4826;
    bars.iter()
        .map(|b| {
            let common = (level * sigma * rng.gauss()).exp();
            let wob = |rng: &mut Rng| (0.5 * level * sigma * rng.gauss()).abs();
            let open = b.open * common;
            let close = b.close * common;
            let high = (b.high * common * (1.0 + wob(rng))).max(open).max(close);
            let low = (b.low * common * (1.0 - wob(rng))).min(open).min(close);
            OHLCV { open, high, low, close, volume: b.volume }
        })
        .collect()
}

/// Names of the expected-movement inputs, in [`magnitude_inputs`] order.
pub const MAGNITUDE_NAMES: &[&str] = &[
    // Baseline: volatility mean reversion alone.
    "ln_vol_ratio_20_126",
    "ln_atr_vs_close_vol",
    // Squeeze structure.
    "energy",
    "depth21",
    "depth10",
    "tier",
    "ln_bars_in_squeeze",
    "fired_x_decay",
    "width_pct",
    // Momentum strength (unsigned).
    "abs_blend",
    "abs_agree",
];

/// How many leading [`MAGNITUDE_NAMES`] form the volatility-only baseline.
pub const MAGNITUDE_BASELINE: usize = 2;

/// Direction-free inputs for the expected-movement model.
pub fn magnitude_inputs(f: &FeatureFrame, i: usize) -> Vec<f64> {
    use crate::sd2::primitives::squash;
    let fired = if f.clock.is_fired[i] { 1.0 } else { 0.0 };
    vec![
        (f.vol20[i] / f.vol126[i]).ln(),
        ((f.atr[i] / f.close[i]) / f.vol20[i].max(1e-6)).ln(),
        squash(f.energy[i], 8.0),
        (2.0 - f.q21[i]).clamp(0.0, 2.0),
        (2.0 - f.q10[i]).clamp(0.0, 2.0),
        f.comp[i] as f64,
        (1.0 + f.clock.bars_since_start[i] as f64).ln(),
        fired * f.decay[i],
        f.width_pct[i],
        f.blend[i].abs(),
        f.agree[i].abs(),
    ]
}

/// Build samples for one symbol.  `signal_bars` is what the indicator sees
/// (maybe noisy); `truth` is the real history used for outcomes.
pub fn build_samples(
    symbol: u16,
    truth: &[DatedBar],
    signal_bars: &[OHLCV],
    cfg: FeatureConfig,
    market: &Market,
    out: &mut Vec<Sample>,
) {
    let f = FeatureFrame::compute(signal_bars, cfg);
    let classic = classic_series(signal_bars);
    let n = truth.len();
    let max_h = *HORIZONS.iter().max().unwrap();
    for i in FeatureFrame::WARMUP..n {
        if !f.row_is_valid(i) || !f.vol20[i].is_finite() || f.vol20[i] <= 0.0 {
            continue;
        }
        let mut x = f.row(i);
        x.extend(market.context(truth[i].date));
        if x.iter().any(|v| !v.is_finite()) {
            continue;
        }
        let mag_x = magnitude_inputs(&f, i);
        if mag_x.iter().any(|v| !v.is_finite()) {
            continue;
        }
        let date = truth[i].date;
        let mut fwd = [f64::NAN; 4];
        let mut excess = [f64::NAN; 4];
        for (k, h) in HORIZONS.iter().enumerate() {
            if i + h < n {
                fwd[k] = (truth[i + h].bar.close / truth[i].bar.close).ln();
                if let (Some(a), Some(b)) = (market.close_on(date), market.close_on(truth[i + h].date)) {
                    excess[k] = fwd[k] - (b / a).ln();
                }
            }
        }
        let barrier = barrier_label(truth, i, f.atr[i]);
        // Purge: drop samples whose longest label would cross a split edge.
        let end_date = truth.get(i + max_h).map(|b| b.date);
        let split = if date < train_end() {
            if end_date.is_none_or(|d| d >= train_end()) { continue } else { Split::Train }
        } else if date < valid_end() {
            if end_date.is_none_or(|d| d >= valid_end()) && i + max_h < n { continue } else { Split::Valid }
        } else {
            Split::Test
        };
        out.push(Sample {
            symbol,
            date,
            split,
            x,
            mag_x,
            fwd,
            barrier,
            excess,
            vol: f.vol20[i],
            close: truth[i].bar.close,
            classic_score: classic.score[i],
            alive: f.fire_alive(i),
            comp: f.comp[i],
            fired_start: f.clock.fired_start[i],
            phase: if f.clock.is_fired[i] { f.clock.bars_since_fired[i] as f64 / f.clock.theta[i].max(1.0) } else { f64::NAN },
            energy: f.energy[i],
            tier21: crate::sd2::features::tier(f.q21[i]),
            tier10: crate::sd2::features::tier(f.q10[i]),
            bars_in_squeeze: f.clock.bars_since_start[i],
        });
    }
}

/// Spearman rank correlation.
pub fn spearman(a: &[f64], b: &[f64]) -> f64 {
    fn ranks(v: &[f64]) -> Vec<f64> {
        let mut idx: Vec<usize> = (0..v.len()).collect();
        idx.sort_by(|&i, &j| v[i].partial_cmp(&v[j]).unwrap_or(std::cmp::Ordering::Equal));
        let mut r = vec![0.0; v.len()];
        let mut k = 0;
        while k < idx.len() {
            let mut e = k;
            while e + 1 < idx.len() && v[idx[e + 1]] == v[idx[k]] {
                e += 1;
            }
            let avg = (k + e) as f64 / 2.0;
            for &i in &idx[k..=e] {
                r[i] = avg;
            }
            k = e + 1;
        }
        r
    }
    pearson(&ranks(a), &ranks(b))
}

pub fn pearson(a: &[f64], b: &[f64]) -> f64 {
    let n = a.len() as f64;
    if n < 3.0 {
        return f64::NAN;
    }
    let (ma, mb) = (a.iter().sum::<f64>() / n, b.iter().sum::<f64>() / n);
    let (mut sab, mut saa, mut sbb) = (0.0, 0.0, 0.0);
    for (x, y) in a.iter().zip(b) {
        sab += (x - ma) * (y - mb);
        saa += (x - ma).powi(2);
        sbb += (y - mb).powi(2);
    }
    sab / (saa * sbb).sqrt()
}

/// Directional scorecard for a set of (signed call strength, realized return).
#[derive(Debug, Clone, Default, serde::Serialize)]
pub struct Scorecard {
    pub calls: usize,
    pub coverage: f64,
    pub hit_rate: f64,
    /// Hit rate of "always call the majority direction" on the same rows.
    pub base_rate: f64,
    pub long_share: f64,
    /// Average of the up-call hit rate and the down-call hit rate.
    pub balanced: f64,
    pub spearman: f64,
    pub mean_signed_return: f64,
}

/// Grade calls.  `strength[i] > 0` is an up call, `< 0` a down call, `0` no
/// call.  Only rows with `|strength| >= threshold` are graded.
pub fn scorecard(strength: &[f64], ret: &[f64], threshold: f64) -> Scorecard {
    let total = strength.len();
    let (mut calls, mut hits, mut ups, mut up_hits, mut downs, mut down_hits, mut real_up) =
        (0, 0, 0, 0, 0, 0, 0);
    let mut signed = 0.0;
    let (mut s_used, mut r_used) = (Vec::new(), Vec::new());
    for (&s, &r) in strength.iter().zip(ret) {
        if !r.is_finite() || s.abs() < threshold || s == 0.0 || r == 0.0 {
            continue;
        }
        calls += 1;
        s_used.push(s);
        r_used.push(r);
        if r > 0.0 {
            real_up += 1;
        }
        let hit = (s > 0.0) == (r > 0.0);
        if hit {
            hits += 1;
        }
        signed += s.signum() * r;
        if s > 0.0 {
            ups += 1;
            up_hits += hit as usize;
        } else {
            downs += 1;
            down_hits += hit as usize;
        }
    }
    let frac = |a: usize, b: usize| if b == 0 { f64::NAN } else { a as f64 / b as f64 };
    let up_rate = frac(up_hits, ups);
    let down_rate = frac(down_hits, downs);
    Scorecard {
        calls,
        coverage: frac(calls, total),
        hit_rate: frac(hits, calls),
        base_rate: frac(real_up.max(calls - real_up), calls),
        long_share: frac(ups, calls),
        balanced: match (up_rate.is_nan(), down_rate.is_nan()) {
            (false, false) => (up_rate + down_rate) / 2.0,
            (false, true) => up_rate,
            (true, false) => down_rate,
            _ => f64::NAN,
        },
        spearman: spearman(&s_used, &r_used),
        mean_signed_return: if calls > 0 { signed / calls as f64 } else { f64::NAN },
    }
}

/// First-touch label.  Uses the *real* future highs and lows.
pub fn barrier_label(truth: &[DatedBar], i: usize, atr: f64) -> f64 {
    let c = truth[i].bar.close;
    let (up, dn) = (c + BARRIER_ATR * atr, c - BARRIER_ATR * atr);
    for b in truth.iter().skip(i + 1).take(BARRIER_BARS) {
        let (hit_up, hit_dn) = (b.bar.high >= up, b.bar.low <= dn);
        match (hit_up, hit_dn) {
            (true, false) => return 1.0,
            (false, true) => return -1.0,
            (true, true) => return 0.0,
            _ => {}
        }
    }
    0.0
}

/// Market context from SPY, joined by date.
///
/// A single stock's squeeze does not fire in a vacuum: in a calm, rising
/// market, breakouts tend to resolve upward; in a falling, volatile one the
/// downside breaks are fast and large.  These three numbers let the model
/// know which weather it is in.
pub struct Market {
    pub dates: Vec<NaiveDate>,
    pub close: Vec<f64>,
    /// [trend vs 200-day average in σ, 126-day time-series momentum, vol ratio 20/126]
    pub ctx: Vec<[f64; 3]>,
}

pub const MARKET_FEATURES: &[&str] = &["mkt_trend200", "mkt_tsmom", "mkt_volratio"];

impl Market {
    pub fn new(spy: &[DatedBar]) -> Self {
        use crate::sd2::primitives::{rolling_std, sma, squash};
        let close: Vec<f64> = spy.iter().map(|b| b.bar.close).collect();
        let n = close.len();
        let lr: Vec<f64> = (0..n).map(|i| if i == 0 { 0.0 } else { (close[i] / close[i - 1]).ln() }).collect();
        let v20 = rolling_std(&lr, 20);
        let v126 = rolling_std(&lr, 126);
        let s200 = sma(&close, 200);
        let ctx = (0..n)
            .map(|i| {
                if i < 200 {
                    return [f64::NAN; 3];
                }
                [
                    squash((close[i] / s200[i]).ln() / (v126[i] * 200f64.sqrt() * 0.5), 1.0),
                    squash((close[i] / close[i - 126]).ln() / (v126[i] * 126f64.sqrt()), 1.0),
                    squash(v20[i] / v126[i] - 1.0, 0.5),
                ]
            })
            .collect();
        Self { dates: spy.iter().map(|b| b.date).collect(), close, ctx }
    }
    fn idx(&self, d: NaiveDate) -> Option<usize> {
        self.dates.binary_search(&d).ok()
    }
    pub fn close_on(&self, d: NaiveDate) -> Option<f64> {
        self.idx(d).map(|i| self.close[i])
    }
    pub fn context(&self, d: NaiveDate) -> [f64; 3] {
        // Use the latest SPY bar at or before `d`.
        match self.dates.binary_search(&d) {
            Ok(i) => self.ctx[i],
            Err(0) => [f64::NAN; 3],
            Err(i) => self.ctx[i - 1],
        }
    }
}

/// Area under the ROC curve: the chance a random "up" row scores higher
/// than a random "down" row.  0.5 = no skill, regardless of base rate.
pub fn auc(score: &[f64], label: &[bool]) -> f64 {
    let mut idx: Vec<usize> = (0..score.len()).collect();
    idx.sort_by(|&a, &b| score[a].partial_cmp(&score[b]).unwrap_or(std::cmp::Ordering::Equal));
    let (mut rank_sum, mut pos) = (0.0, 0.0);
    for (r, &i) in idx.iter().enumerate() {
        if label[i] {
            rank_sum += (r + 1) as f64;
            pos += 1.0;
        }
    }
    let neg = score.len() as f64 - pos;
    if pos == 0.0 || neg == 0.0 {
        return f64::NAN;
    }
    (rank_sum - pos * (pos + 1.0) / 2.0) / (pos * neg)
}