portfolio_analyzer/sd2/

model.rs

//! # Model: how much is each measurement worth?
//!
//! Classic SD hard-coded its opinions: Group A `+` is worth 1.0, `(+)` is
//! worth 0.75, every `*` adds 0.05.  Those numbers were sensible guesses, but
//! nobody had checked them against what prices actually did next.
//!
//! SD v2 keeps the *structure* (the features in `features.rs` are all squeeze
//! and momentum ideas from the paper) but lets history set the *weights*.
//! Two tiny, transparent models do this:
//!
//! * [`Logistic`]: "given these measurements, what is the chance the next
//!   `h` bars close higher?"  Its output is the SD number, in `[0, 1]`.
//! * [`Ridge`]: "how big will the move be, in units of normal daily noise?"
//!   This is the *expected movement* that SD's conviction is scaled by.
//!
//! Both are linear, so every weight can be printed, read, argued with and
//! copied into ThinkScript.  Both use an L2 ("ridge") penalty, which pulls
//! weights toward zero unless the data insists, the cheapest defence there
//! is against fitting noise.

use serde::{Deserialize, Serialize};

/// Solve `A x = b` for a small symmetric positive-definite `A` (Cholesky).
fn solve_spd(a: &[Vec<f64>], b: &[f64]) -> Vec<f64> {
    let k = b.len();
    let mut l = vec![vec![0.0; k]; k];
    for i in 0..k {
        for j in 0..=i {
            let s: f64 = (0..j).map(|p| l[i][p] * l[j][p]).sum();
            if i == j {
                l[i][j] = (a[i][i] - s).max(1e-12).sqrt();
            } else {
                l[i][j] = (a[i][j] - s) / l[j][j];
            }
        }
    }
    let mut y = vec![0.0; k];
    for i in 0..k {
        y[i] = (b[i] - (0..i).map(|p| l[i][p] * y[p]).sum::<f64>()) / l[i][i];
    }
    let mut x = vec![0.0; k];
    for i in (0..k).rev() {
        x[i] = (y[i] - (i + 1..k).map(|p| l[p][i] * x[p]).sum::<f64>()) / l[i][i];
    }
    x
}

/// Per-column mean and standard deviation, so weights are comparable.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Scaler {
    pub mean: Vec<f64>,
    pub std: Vec<f64>,
}

impl Scaler {
    pub fn fit(x: &[Vec<f64>]) -> Self {
        let k = x[0].len();
        let n = x.len() as f64;
        let mean: Vec<f64> = (0..k).map(|j| x.iter().map(|r| r[j]).sum::<f64>() / n).collect();
        let std = (0..k)
            .map(|j| {
                let v = x.iter().map(|r| (r[j] - mean[j]).powi(2)).sum::<f64>() / n;
                v.sqrt().max(1e-9)
            })
            .collect();
        Self { mean, std }
    }
    pub fn apply(&self, r: &[f64]) -> Vec<f64> {
        r.iter().enumerate().map(|(j, v)| (v - self.mean[j]) / self.std[j]).collect()
    }
}

/// Logistic regression `p = σ(b₀ + w·z)` on standardized inputs `z`.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Logistic {
    pub scaler: Scaler,
    pub intercept: f64,
    pub weights: Vec<f64>,
}

pub fn sigmoid(x: f64) -> f64 {
    1.0 / (1.0 + (-x).exp())
}

impl Logistic {
    /// Fit by Newton's method (IRLS) with ridge penalty `lambda` per row.
    /// Samples can be weighted (e.g. to stop one symbol dominating).
    pub fn fit(x: &[Vec<f64>], y: &[f64], w: &[f64], lambda: f64) -> Self {
        let scaler = Scaler::fit(x);
        let z: Vec<Vec<f64>> = x.iter().map(|r| scaler.apply(r)).collect();
        let k = z[0].len() + 1;
        let total_w: f64 = w.iter().sum();
        let mut beta = vec![0.0; k];
        for _ in 0..25 {
            let mut h = vec![vec![0.0; k]; k];
            let mut g = vec![0.0; k];
            for ((row, &yi), &wi) in z.iter().zip(y).zip(w) {
                let eta = beta[0] + row.iter().zip(&beta[1..]).map(|(a, b)| a * b).sum::<f64>();
                let p = sigmoid(eta);
                let s = wi * p * (1.0 - p);
                let r = wi * (yi - p);
                let f = |j: usize| if j == 0 { 1.0 } else { row[j - 1] };
                for a in 0..k {
                    g[a] += r * f(a);
                    for b in 0..=a {
                        h[a][b] += s * f(a) * f(b);
                    }
                }
            }
            for a in 0..k {
                for b in 0..a {
                    h[b][a] = h[a][b];
                }
                if a > 0 {
                    h[a][a] += lambda * total_w;
                    g[a] -= lambda * total_w * beta[a];
                }
            }
            let step = solve_spd(&h, &g);
            let size: f64 = step.iter().map(|s| s * s).sum::<f64>().sqrt();
            for (b, s) in beta.iter_mut().zip(&step) {
                *b += s;
            }
            if size < 1e-8 {
                break;
            }
        }
        Self { scaler, intercept: beta[0], weights: beta[1..].to_vec() }
    }

    pub fn logit(&self, r: &[f64]) -> f64 {
        let z = self.scaler.apply(r);
        self.intercept + z.iter().zip(&self.weights).map(|(a, b)| a * b).sum::<f64>()
    }

    pub fn predict(&self, r: &[f64]) -> f64 {
        sigmoid(self.logit(r))
    }

    /// The part of the logit that comes from the features, not the base rate.
    /// Positive means "the measurements lean up".
    pub fn lean(&self, r: &[f64]) -> f64 {
        self.logit(r) - self.intercept
    }
}

/// Ridge regression `y = b₀ + w·z`.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Ridge {
    pub scaler: Scaler,
    pub intercept: f64,
    pub weights: Vec<f64>,
}

impl Ridge {
    pub fn fit(x: &[Vec<f64>], y: &[f64], lambda: f64) -> Self {
        let scaler = Scaler::fit(x);
        let k = x[0].len() + 1;
        let n = x.len() as f64;
        let mut a = vec![vec![0.0; k]; k];
        let mut b = vec![0.0; k];
        for (r, &yi) in x.iter().zip(y) {
            let z = scaler.apply(r);
            let f = |j: usize| if j == 0 { 1.0 } else { z[j - 1] };
            for i in 0..k {
                b[i] += f(i) * yi;
                for j in 0..=i {
                    a[i][j] += f(i) * f(j);
                }
            }
        }
        for i in 0..k {
            for j in 0..i {
                a[j][i] = a[i][j];
            }
            if i > 0 {
                a[i][i] += lambda * n;
            }
        }
        let beta = solve_spd(&a, &b);
        Self { scaler, intercept: beta[0], weights: beta[1..].to_vec() }
    }

    pub fn predict(&self, r: &[f64]) -> f64 {
        let z = self.scaler.apply(r);
        self.intercept + z.iter().zip(&self.weights).map(|(a, b)| a * b).sum::<f64>()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn logistic_recovers_a_clear_signal() {
        let x: Vec<Vec<f64>> = (0..2000).map(|i| vec![((i * 37) % 200) as f64 / 100.0 - 1.0]).collect();
        let y: Vec<f64> = x.iter().map(|r| if r[0] > 0.0 { 1.0 } else { 0.0 }).collect();
        let m = Logistic::fit(&x, &y, &vec![1.0; 2000], 1e-4);
        assert!(m.predict(&[0.8]) > 0.9 && m.predict(&[-0.8]) < 0.1);
    }

    #[test]
    fn ridge_recovers_a_line() {
        let x: Vec<Vec<f64>> = (0..100).map(|i| vec![i as f64]).collect();
        let y: Vec<f64> = x.iter().map(|r| 3.0 * r[0] + 2.0).collect();
        let m = Ridge::fit(&x, &y, 0.0);
        assert!((m.predict(&[10.0]) - 32.0).abs() < 1e-6);
    }
}