portfolio_analyzer/

direction.rs

//! # Direction: what is the best SD can do at calling *up or down*?
//!
//! Chapter 2 of the book found that squeeze and momentum measurements, on
//! their own, call direction no better than a coin flip once the market's
//! upward drift is removed.  This module widens the search to every family
//! of measurement that the finance literature says carries *some* direction
//! information, and lets a walk-forward test (`sd_direction`) decide which
//! of them earn their place.
//!
//! The families, and the one-line reason each is here:
//!
//! | Group | Measurements | Why it might know direction |
//! | --- | --- | --- |
//! | `sd` | blend, agreement, fired×decay×blend, tier, energy, stretch, semi | SD's own reading |
//! | `mom` | 5/21/63/126-bar returns, 12-minus-1-month return | trends persist; very short moves reverse |
//! | `trend` | distance from 50/200-bar averages, drawdown from 1-year high | trend-following, anchoring |
//! | `vol` | 20/126 vol ratio, vol level | calm markets drift up; panics overshoot |
//! | `rel` | 63-bar return minus SPY's | relative strength |
//! | `mkt` | SPY trend, momentum, vol ratio, 5/21-bar return | the tide lifts (and sinks) most boats |
//! | `vix` | VIX vs its 1-year average, VIX level, VIX minus realized vol | fear is paid for: high implied vol precedes better returns |
//! | `season` | November–April flag, month sine/cosine | the "Halloween" seasonal |
//! | `sdw`, `sdm` | SD's reading on weekly and monthly bars | Extended SD's idea: slower squeezes, slower trends |
//!
//! Every value at bar `t` uses only bars `<= t`.

use crate::history::DatedBar;
use crate::indicators::{rolling_std, sma, OHLCV};
use crate::sd2::features::FeatureFrame;
use crate::sd2::primitives::{highest, squash};
use chrono::{Datelike, NaiveDate};

/// Feature groups, in column order.
pub const GROUPS: &[(&str, &[&str])] = &[
    ("sd", &["sd_blend", "sd_agree", "sd_fire_blend", "sd_comp", "sd_energy", "sd_stretch", "sd_semi"]),
    ("mom", &["r5", "r21", "r63", "r126", "mom12_1"]),
    ("trend", &["d50", "d200", "dd252"]),
    ("vol", &["vol_ratio", "vol_level"]),
    ("rel", &["rel63"]),
    ("mkt", &["spy_trend200", "spy_tsmom", "spy_volratio", "spy_r5", "spy_r21"]),
    ("vix", &["vix_rel", "vix_level", "vix_vrp"]),
    ("season", &["halloween", "month_sin", "month_cos"]),
    ("sdw", &["w_blend", "w_agree", "w_fire_blend", "w_comp", "w_energy", "w_stretch"]),
    ("sdm", &["m_blend", "m_agree", "m_fire_blend", "m_comp", "m_energy", "m_stretch"]),
];

pub fn n_features() -> usize {
    GROUPS.iter().map(|g| g.1.len()).sum()
}

/// Column range of a group.
pub fn group_range(name: &str) -> std::ops::Range<usize> {
    let mut start = 0;
    for (g, cols) in GROUPS {
        if *g == name {
            return start..start + cols.len();
        }
        start += cols.len();
    }
    panic!("unknown feature group {name}");
}

fn log_returns(close: &[f64]) -> Vec<f64> {
    (0..close.len()).map(|i| if i == 0 { 0.0 } else { (close[i] / close[i - 1]).ln() }).collect()
}

/// A move of `h` bars in units of its normal size (`σ₁₂₆·√h`), squashed.
fn norm_ret(close: &[f64], v126: &[f64], i: usize, from: usize, h: usize) -> f64 {
    if i < from + h || !v126[i].is_finite() || v126[i] <= 0.0 {
        return f64::NAN;
    }
    squash((close[i - from] / close[i - from - h]).ln() / (v126[i] * (h as f64).sqrt()), 2.0)
}

/// Market weather on each SPY date: SPY's own trend and VIX's fear gauge.
pub struct MarketWeather {
    pub dates: Vec<NaiveDate>,
    pub rows: Vec<[f64; 8]>,
    pub close: Vec<f64>,
}

impl MarketWeather {
    /// `vix` may be empty (then the VIX columns are neutral zeros).
    pub fn new(spy: &[DatedBar], vix: &[DatedBar]) -> Self {
        let close: Vec<f64> = spy.iter().map(|b| b.bar.close).collect();
        let n = close.len();
        let lr = log_returns(&close);
        let v20 = rolling_std(&lr, 20);
        let v126 = rolling_std(&lr, 126);
        let s200 = sma(&close, 200);
        let vix_dates: Vec<NaiveDate> = vix.iter().map(|b| b.date).collect();
        let vix_close: Vec<f64> = vix.iter().map(|b| b.bar.close).collect();
        let vix_mean = sma(&vix_close, 252);
        let rows = (0..n)
            .map(|i| {
                if i < 252 {
                    return [f64::NAN; 8];
                }
                let j = match vix_dates.binary_search(&spy[i].date) {
                    Ok(j) => Some(j),
                    Err(0) => None,
                    Err(j) => Some(j - 1),
                };
                let (vrel, vlvl, vrp) = match j {
                    Some(j) if j >= 252 && vix_close[j] > 0.0 => {
                        let v = vix_close[j];
                        (
                            squash((v / vix_mean[j]).ln(), 0.5),
                            squash((v / 20.0).ln(), 0.5),
                            squash((v / 100.0 / (v20[i] * 252f64.sqrt())).ln(), 0.5),
                        )
                    }
                    _ => (0.0, 0.0, 0.0),
                };
                [
                    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),
                    norm_ret(&close, &v126, i, 0, 5),
                    norm_ret(&close, &v126, i, 0, 21),
                    vrel,
                    vlvl,
                    vrp,
                ]
            })
            .collect();
        Self { dates: spy.iter().map(|b| b.date).collect(), rows, close }
    }

    /// Latest row at or before `d`.
    pub fn at(&self, d: NaiveDate) -> Option<(usize, [f64; 8])> {
        let i = match self.dates.binary_search(&d) {
            Ok(i) => i,
            Err(0) => return None,
            Err(i) => i - 1,
        };
        Some((i, self.rows[i]))
    }
}

/// SD's reading on slower bars (weekly or monthly), mapped back to days.
///
/// Day `i` sees only *completed* slow bars: the last week (or month) that
/// ended before the one day `i` belongs to.  Unusable values become 0
/// (neutral), so a young listing is not dropped, merely uninformed.
pub fn slow_sd(bars: &[OHLCV], dates: &[NaiveDate], monthly: bool) -> Vec<[f64; 6]> {
    let key = |d: &NaiveDate| if monthly { d.year() * 12 + d.month() as i32 } else { d.iso_week().year() * 53 + d.iso_week().week() as i32 };
    let mut agg: Vec<OHLCV> = Vec::new();
    let mut period_of_day = Vec::with_capacity(bars.len());
    let mut last_key = None;
    for (b, d) in bars.iter().zip(dates) {
        let k = key(d);
        if last_key != Some(k) {
            agg.push(b.clone());
            last_key = Some(k);
        } else {
            let a = agg.last_mut().unwrap();
            a.high = a.high.max(b.high);
            a.low = a.low.min(b.low);
            a.close = b.close;
            a.volume += b.volume;
        }
        period_of_day.push(agg.len() - 1);
    }
    let f = FeatureFrame::compute(&agg, crate::sd2::features::FeatureConfig::default());
    let clean = |v: f64| if v.is_finite() { v } else { 0.0 };
    period_of_day
        .into_iter()
        .map(|p| {
            if p == 0 || p - 1 < 30 {
                return [0.0; 6];
            }
            let j = p - 1;
            let alive = if f.fire_alive(j) { 1.0 } else { 0.0 };
            [
                clean(f.blend[j]),
                clean(f.agree[j]),
                clean(alive * f.blend[j] * f.decay[j]),
                f.comp[j] as f64 / 3.0,
                clean(squash(f.energy[j], 8.0)),
                clean(f.stretch[j]),
            ]
        })
        .collect()
}

/// Every direction feature for every bar of one symbol (`NaN` = unusable).
pub fn direction_frame(bars: &[OHLCV], dates: &[NaiveDate], f: &FeatureFrame, mkt: &MarketWeather) -> Vec<Vec<f64>> {
    let close: Vec<f64> = bars.iter().map(|b| b.close).collect();
    let n = close.len();
    let lr = log_returns(&close);
    let v20 = rolling_std(&lr, 20);
    let v126 = rolling_std(&lr, 126);
    let s50 = sma(&close, 50);
    let s200 = sma(&close, 200);
    let hi252 = highest(&close, 252);
    let weekly = slow_sd(bars, dates, false);
    let monthly = slow_sd(bars, dates, true);
    (0..n)
        .map(|i| {
            if i < 260 || !v126[i].is_finite() || v126[i] <= 0.0 {
                return vec![f64::NAN; n_features()];
            }
            let alive = if f.fire_alive(i) { 1.0 } else { 0.0 };
            let mut r = vec![
                f.blend[i],
                f.agree[i],
                alive * f.blend[i] * f.decay[i],
                f.comp[i] as f64 / 3.0,
                squash(f.energy[i], 8.0),
                f.stretch[i],
                f.semi[i],
            ];
            r.extend([
                norm_ret(&close, &v126, i, 0, 5),
                norm_ret(&close, &v126, i, 0, 21),
                norm_ret(&close, &v126, i, 0, 63),
                norm_ret(&close, &v126, i, 0, 126),
                norm_ret(&close, &v126, i, 21, 231),
            ]);
            r.extend([
                squash((close[i] / s50[i]).ln() / (v126[i] * 25f64.sqrt()), 2.0),
                squash((close[i] / s200[i]).ln() / (v126[i] * 100f64.sqrt()), 2.0),
                squash((close[i] / hi252[i]).ln() / (v126[i] * 252f64.sqrt()), 1.0),
            ]);
            r.extend([squash((v20[i] / v126[i]).ln(), 0.5), squash((v126[i] * 252f64.sqrt() / 0.25).ln(), 1.0)]);
            let (mi, m) = match mkt.at(dates[i]) {
                Some(x) => x,
                None => return vec![f64::NAN; n_features()],
            };
            let spy63 = if mi >= 63 { (mkt.close[mi] / mkt.close[mi - 63]).ln() } else { f64::NAN };
            r.push(squash(((close[i] / close[i - 63]).ln() - spy63) / (v126[i] * 63f64.sqrt()), 2.0));
            r.extend(&m[0..5]);
            r.extend(&m[5..8]);
            let month = dates[i].month() as f64;
            let angle = std::f64::consts::TAU * (month - 0.5) / 12.0;
            r.extend([if (5..=10).contains(&dates[i].month()) { 0.0 } else { 1.0 }, angle.sin(), angle.cos()]);
            r.extend(weekly[i]);
            r.extend(monthly[i]);
            r
        })
        .collect()
}

/// A small gradient-boosted tree ensemble for "will it close higher?".
///
/// The linear model can only say "more of X, more up".  Direction might
/// instead live in *regimes*: "falling market **and** fear spiking **and**
/// this stock far below its average".  Shallow trees find such corners.
/// Kept tiny (depth 3, coarse 32-bin thresholds, big leaves) because
/// direction data is mostly noise and big trees would memorise it.
#[derive(Debug, Clone)]
pub struct Gbm {
    pub base: f64,
    pub edges: Vec<Vec<f64>>,
    pub trees: Vec<Vec<GbmNode>>,
}

#[derive(Debug, Clone)]
pub enum GbmNode {
    /// Rows with `bin(feature) <= bin` go to `left`, others to `right`.
    Split { feature: usize, bin: u8, left: usize, right: usize },
    Leaf(f64),
}

pub struct GbmParams {
    pub rounds: usize,
    pub learning_rate: f64,
    pub depth: usize,
    /// Smallest leaf, as a share of training rows.
    pub min_leaf: f64,
}

impl Default for GbmParams {
    fn default() -> Self {
        Self { rounds: 150, learning_rate: 0.05, depth: 3, min_leaf: 0.01 }
    }
}

const BINS: usize = 32;

fn bin_of(edges: &[f64], v: f64) -> u8 {
    edges.partition_point(|&e| e < v) as u8
}

impl Gbm {
    pub fn fit(x: &[Vec<f64>], y: &[f64], p: &GbmParams) -> Self {
        let n = x.len();
        let d = x[0].len();
        // Quantile thresholds from a strided sample.
        let step = (n / 20_000).max(1);
        let edges: Vec<Vec<f64>> = (0..d)
            .map(|j| {
                let mut v: Vec<f64> = x.iter().step_by(step).map(|r| r[j]).collect();
                v.sort_by(|a, b| a.partial_cmp(b).unwrap());
                let mut e: Vec<f64> = (1..BINS).map(|k| v[k * (v.len() - 1) / BINS]).collect();
                e.dedup();
                e
            })
            .collect();
        let binned: Vec<Vec<u8>> = x.iter().map(|r| r.iter().enumerate().map(|(j, &v)| bin_of(&edges[j], v)).collect()).collect();
        let mean = (y.iter().sum::<f64>() / n as f64).clamp(1e-6, 1.0 - 1e-6);
        let base = (mean / (1.0 - mean)).ln();
        let mut f = vec![base; n];
        let min_leaf = ((n as f64 * p.min_leaf) as usize).max(50);
        let mut trees = Vec::with_capacity(p.rounds);
        for _ in 0..p.rounds {
            let prob: Vec<f64> = f.iter().map(|&v| 1.0 / (1.0 + (-v).exp())).collect();
            let g: Vec<f64> = (0..n).map(|i| y[i] - prob[i]).collect();
            let h: Vec<f64> = prob.iter().map(|&q| (q * (1.0 - q)).max(1e-6)).collect();
            let mut nodes = Vec::new();
            let idx: Vec<u32> = (0..n as u32).collect();
            grow(&binned, &g, &h, idx, 0, p, min_leaf, d, &mut nodes);
            for i in 0..n {
                f[i] += leaf_value(&nodes, &binned[i]);
            }
            trees.push(nodes);
        }
        Self { base, edges, trees }
    }

    pub fn predict(&self, r: &[f64]) -> f64 {
        let b: Vec<u8> = r.iter().enumerate().map(|(j, &v)| bin_of(&self.edges[j], v)).collect();
        let z = self.base + self.trees.iter().map(|t| leaf_value(t, &b)).sum::<f64>();
        1.0 / (1.0 + (-z).exp())
    }
}

fn leaf_value(nodes: &[GbmNode], b: &[u8]) -> f64 {
    let mut k = 0;
    loop {
        match nodes[k] {
            GbmNode::Leaf(v) => return v,
            GbmNode::Split { feature, bin, left, right } => k = if b[feature] <= bin { left } else { right },
        }
    }
}

#[allow(clippy::too_many_arguments)]
fn grow(
    x: &[Vec<u8>],
    g: &[f64],
    h: &[f64],
    idx: Vec<u32>,
    depth: usize,
    p: &GbmParams,
    min_leaf: usize,
    d: usize,
    nodes: &mut Vec<GbmNode>,
) -> usize {
    const L2: f64 = 1.0;
    let me = nodes.len();
    let (gs, hs) = idx.iter().fold((0.0, 0.0), |(a, b), &i| (a + g[i as usize], b + h[i as usize]));
    nodes.push(GbmNode::Leaf(p.learning_rate * gs / (hs + L2)));
    if depth >= p.depth || idx.len() < 2 * min_leaf {
        return me;
    }
    let parent = gs * gs / (hs + L2);
    let mut best: Option<(f64, usize, u8)> = None;
    for j in 0..d {
        let mut hist = [(0.0f64, 0.0f64, 0usize); BINS];
        for &i in &idx {
            let e = &mut hist[x[i as usize][j] as usize];
            e.0 += g[i as usize];
            e.1 += h[i as usize];
            e.2 += 1;
        }
        let (mut gl, mut hl, mut nl) = (0.0, 0.0, 0usize);
        for (b, e) in hist.iter().enumerate().take(BINS - 1) {
            gl += e.0;
            hl += e.1;
            nl += e.2;
            let nr = idx.len() - nl;
            if nl < min_leaf || nr < min_leaf {
                continue;
            }
            let (gr, hr) = (gs - gl, hs - hl);
            let gain = gl * gl / (hl + L2) + gr * gr / (hr + L2) - parent;
            if best.is_none_or(|(bg, _, _)| gain > bg) {
                best = Some((gain, j, b as u8));
            }
        }
    }
    let Some((gain, feature, bin)) = best else { return me };
    if gain <= 0.0 {
        return me;
    }
    let (li, ri): (Vec<u32>, Vec<u32>) = idx.into_iter().partition(|&i| x[i as usize][feature] <= bin);
    let left = grow(x, g, h, li, depth + 1, p, min_leaf, d, nodes);
    let right = grow(x, g, h, ri, depth + 1, p, min_leaf, d, nodes);
    nodes[me] = GbmNode::Split { feature, bin, left, right };
    me
}