portfolio_analyzer/sd2/

engine.rs

//! # Engine: from measurements to the SD number
//!
//! This is the file to open if you want to know what the number on the chart
//! *means*.  Everything upstream measures; this file decides.
//!
//! ## What the backtests taught us (and why the engine looks like this)
//!
//! `sd_backtest study` and `sd_backtest grid` asked, over 125 symbols and
//! 25 years, trained before 2014 and graded after 2019:
//!
//! 1. **Can squeeze + momentum features call direction?**  Barely.  Every
//!    feature, alone or learned in combination, lands within a point or two
//!    of a coin flip once the market's upward drift is removed.  Classic SD's
//!    `+`/`-` calls hit 50.5%, while "always up" hit 54-61%.
//! 2. **Can they call the size of the move?**  Yes.  Predicted versus
//!    realized move size has a rank correlation of about 0.22 out of
//!    sample, and the effect is monotone in squeeze tier, stored energy and
//!    squeeze duration in *both* eras.
//!
//! So SD v2 is built around the thing it can actually do.  The SD number
//! keeps its familiar shape (0 = bearish, 0.5 = nothing, 1 = bullish), but:
//!
//! * **how far** SD sits from 0.5 is set by the *expected movement*
//!   (bigger expected move → further from 0.5), and
//! * **which side** of 0.5 it sits on is the best direction call the
//!   walk-forward search (`sd_direction`, chapter 8) could find: P(up) from
//!   drift plus the market's weather.  That is right about 58% of the time
//!   at 21 bars, the same as "always up", because nothing beat drift.  The
//!   momentum lean is still reported, but it no longer picks the side.
//!
//! ```text
//!   SD = 0.5 + 0.5 · side · strength · intensity
//!        side      ∈ {-1, +1}  sign of P(up) − 0.5
//!        strength  ∈ [ 0, 1]   time-decay strength of the squeeze
//!        intensity ∈ [ 0, 1]   where this bar's expected move ranks
//! ```

use super::features::{decay_kernel, FeatureConfig, FeatureFrame};
use super::model::{sigmoid, Logistic, Ridge};
use super::primitives::OHLCV;
use super::squeeze::group_c;
use crate::research::{magnitude_inputs, MAGNITUDE_NAMES};
use serde::{Deserialize, Serialize};

/// The frozen, fitted pieces SD v2 needs at run time.  Produced by
/// `sd_backtest fit` and embedded in the binary from `sd2_model.json`.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Sd2Model {
    /// Horizon, in bars, that `magnitude` was trained for.
    pub horizon: usize,
    /// Predicts |move| / (σ·√h) from [`MAGNITUDE_NAMES`] inputs.
    pub magnitude: Ridge,
    /// Expected-move values (σ units) at the 10th and 90th percentile of the
    /// training data; intensity maps this range onto [0, 1].
    pub intensity_lo: f64,
    pub intensity_hi: f64,
    /// Observed share of up moves by SD bucket (0.0-0.1, …, 0.9-1.0) on the
    /// training data.  This is the *honest* probability to show a person.
    pub calibration: Vec<f64>,
    /// Implied-volatility proxy for symbols without a listed vol index:
    /// `ln IV = a + b·ln RV20 + c·ln RV126` (annualized), fitted on SPY,
    /// QQQ, DIA, GLD and USO against their Cboe indices.
    pub iv_proxy: [f64; 3],
    pub feature_names: Vec<String>,
    pub fitted_on: String,
    /// Which side of 0.5 SD sits on: P(close higher in `horizon` bars) from
    /// the market's weather (SPY trend, momentum, vol ratio, 5- and 21-bar
    /// return; see `direction.rs`).  Written by `sd_direction --ship`.
    #[serde(default)]
    pub side: Option<Logistic>,
}

impl Sd2Model {
    /// The model shipped with this crate.
    pub fn embedded() -> Self {
        serde_json::from_str(include_str!("sd2_model.json")).expect("valid embedded sd2_model.json")
    }

    pub fn intensity(&self, expected_sigma: f64) -> f64 {
        ((expected_sigma - self.intensity_lo) / (self.intensity_hi - self.intensity_lo)).clamp(0.0, 1.0)
    }

    /// P(up) from market weather; without weather, the average day's answer.
    pub fn p_up(&self, weather: Option<&[f64]>) -> f64 {
        match (&self.side, weather) {
            (Some(m), Some(w)) if w.iter().all(|v| v.is_finite()) => m.predict(w),
            (Some(m), _) => sigmoid(m.intercept),
            (None, _) => 0.5,
        }
    }

    pub fn calibrated_up(&self, sd: f64) -> f64 {
        let b = ((sd * 10.0).floor() as usize).min(self.calibration.len().saturating_sub(1));
        self.calibration.get(b).copied().unwrap_or(0.5)
    }

    /// Implied-volatility proxy (annualized) from daily realized vols.
    pub fn iv_proxy(&self, vol20_daily: f64, vol126_daily: f64) -> f64 {
        let a = 252f64.sqrt();
        let [c0, c1, c2] = self.iv_proxy;
        (c0 + c1 * (vol20_daily * a).ln() + c2 * (vol126_daily * a).ln()).exp()
    }
}

/// One SD v2 reading.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Sd2Reading {
    /// The SD number, `[0, 1]`.
    pub sd: f64,
    /// Momentum lean after time decay, `[-1, 1]`.  Shown for context; it
    /// no longer sets SD's side (it calls direction at coin-flip rates).
    pub lean: f64,
    /// Probability of a higher close in `horizon` bars (market weather + drift).
    pub p_up: f64,
    /// SD's side: +1 above 0.5, −1 below (from `p_up`).
    pub side: f64,
    /// Expected-move rank, `[0, 1]`.
    pub intensity: f64,
    /// Expected |move| over `horizon` bars, in σ·√h units (≈0.8 is ordinary).
    pub expected_move_sigma: f64,
    /// The same, as a fraction of price (0.05 = 5%).
    pub expected_move_pct: f64,
    pub horizon: usize,
    /// Historical share of up moves for readings like this one.
    pub calibrated_up_probability: f64,
    /// Time-decay strength of the current fire, `[0, 1]`.
    pub decay_strength: f64,
    /// τ/θ, or `None` if no squeeze has fired.
    pub phase: Option<f64>,
    pub compression: u8,
    pub energy: f64,
    pub group_c: u8,
    pub is_active_fire: bool,
    pub in_squeeze: bool,
    /// Daily volatility (20 and 126 bars) at this bar.
    pub vol20: f64,
    pub vol126: f64,
    pub close: f64,
    /// How far price is from its 21-bar EMA, squashed to (−1, 1).
    #[serde(default)]
    pub stretch: f64,
}

/// How much to trust the momentum lean at this point in the squeeze's life.
///
/// * A live fire gets the full decay curve: 1.0 while young, fading after.
/// * Inside a squeeze (not yet fired), direction is genuinely unknown; the
///   lean is shown at half strength so the chart still hints which way
///   momentum is leaning while the coil builds.
/// * Otherwise the lean is shown at a quarter strength: momentum without a
///   squeeze behind it is the weakest kind of evidence.
pub fn time_strength(f: &FeatureFrame, i: usize) -> f64 {
    if f.clock.is_active[i] {
        decay_kernel(f.clock.bars_since_fired[i] as f64 / f.clock.theta[i].max(1.0))
    } else if f.comp[i] > 0 {
        0.5
    } else if f.clock.is_fired[i] {
        0.25 * decay_kernel(f.clock.bars_since_fired[i] as f64 / f.clock.theta[i].max(1.0))
    } else {
        0.25
    }
}

/// Which daily volatility turns σ-units into a % move.
///
/// Measurement noise can only *inflate* the 20-bar estimate, so taking the
/// smaller of it and the noise-robust estimate guards against bad prints.
/// `SD_VOL=vol20` / `SD_VOL=robust` force one estimator for A/B tests.
fn move_vol(f: &FeatureFrame, i: usize) -> f64 {
    let robust = f.vol_robust[i];
    match std::env::var("SD_VOL").as_deref() {
        _ if !robust.is_finite() => f.vol20[i],
        Ok("vol20") => f.vol20[i],
        Ok("robust") => robust,
        _ => f.vol20[i].min(robust),
    }
}

/// The SD number: side × how much to trust it now × how big the move will be.
pub fn side_sd(side: f64, strength: f64, intensity: f64) -> f64 {
    (0.5 + 0.5 * side * strength * intensity).clamp(0.0, 1.0)
}

/// Re-decide a reading's side with today's market weather (SPY row from
/// [`crate::direction::MarketWeather`]; the first five columns).
pub fn apply_weather(r: &mut Sd2Reading, model: &Sd2Model, weather: Option<&[f64]>) {
    r.p_up = model.p_up(weather);
    r.side = if r.p_up >= 0.5 { 1.0 } else { -1.0 };
    r.sd = side_sd(r.side, r.decay_strength, r.intensity);
    r.calibrated_up_probability = r.p_up;
}

/// SD v2 for every bar, with the market's weather deciding the side.
pub fn sd2_series_market(
    bars: &[OHLCV],
    dates: &[chrono::NaiveDate],
    model: &Sd2Model,
    mkt: &crate::direction::MarketWeather,
) -> Vec<Option<Sd2Reading>> {
    let mut out = sd2_series(bars, model);
    for (r, d) in out.iter_mut().zip(dates) {
        if let Some(r) = r {
            let w = mkt.at(*d).map(|(_, row)| row[..5].to_vec());
            apply_weather(r, model, w.as_deref());
        }
    }
    out
}

/// Compute SD v2 for every bar of `bars`.
pub fn sd2_series(bars: &[OHLCV], model: &Sd2Model) -> Vec<Option<Sd2Reading>> {
    let f = FeatureFrame::compute(bars, FeatureConfig::default());
    (0..f.n).map(|i| reading_at(&f, i, model)).collect()
}

/// Compute SD v2 for the last bar only.
pub fn sd2_latest(bars: &[OHLCV], model: &Sd2Model) -> Option<Sd2Reading> {
    let f = FeatureFrame::compute(bars, FeatureConfig::default());
    reading_at(&f, f.n.checked_sub(1)?, model)
}

pub fn reading_at(f: &FeatureFrame, i: usize, model: &Sd2Model) -> Option<Sd2Reading> {
    if i < 140 || !f.atr[i].is_finite() || !f.vol20[i].is_finite() || !f.vol126[i].is_finite() {
        return None;
    }
    let mut mag_x = magnitude_inputs(f, i);
    // Under a year of history the bandwidth percentile is unknown; use the
    // neutral middle rather than refusing to answer.
    if !mag_x[8].is_finite() {
        mag_x[8] = 0.5;
    }
    if mag_x.iter().any(|v| !v.is_finite()) || mag_x.len() != MAGNITUDE_NAMES.len() {
        return None;
    }
    let expected = model.magnitude.predict(&mag_x).max(0.0);
    let intensity = model.intensity(expected);
    let strength = time_strength(f, i);
    let lean = (2.0 * f.blend[i]).tanh() * strength;
    let p_up = model.p_up(None);
    let side = if p_up >= 0.5 { 1.0 } else { -1.0 };
    let sd = side_sd(side, strength, intensity);
    let active = f.clock.is_active[i];
    let eff = if active { f.clock.highest_comp[i] } else { f.comp[i] };
    Some(Sd2Reading {
        sd,
        lean,
        p_up,
        side,
        intensity,
        expected_move_sigma: expected,
        expected_move_pct: expected * move_vol(f, i) * (model.horizon as f64).sqrt(),
        horizon: model.horizon,
        calibrated_up_probability: model.calibrated_up(sd),
        decay_strength: strength,
        phase: f.clock.is_fired[i].then(|| f.clock.bars_since_fired[i] as f64 / f.clock.theta[i].max(1.0)),
        compression: f.comp[i],
        energy: f.energy[i],
        group_c: group_c(eff, f.comp[i]),
        is_active_fire: active,
        in_squeeze: f.comp[i] > 0,
        vol20: f.vol20[i],
        vol126: f.vol126[i],
        close: f.close[i],
        stretch: f.stretch[i],
    })
}