chore: prepare v1.1.0 release

Update version numbers across Rust, Python, and documentation files to 1.1.0. Enhance the .gitignore to include macOS dSYM files and plans directory. Introduce new dependencies in the Rust core library and update the README to reflect recent performance benchmarks and backtesting engine capabilities. Add new artifacts to the benchmarks manifest and improve documentation for the backtesting engine API.

crates/ferro_ta_core/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ferro_ta_core"
-version = "1.0.6"
+version = "1.2.0"
 edition = "2021"
 description = "Pure Rust core indicator library — no PyO3, no numpy dependency"
 license = "MIT"
@@ -17,6 +17,8 @@ crate-type = ["lib"]
 
 [dependencies]
 wide = { version = "1.1.1", optional = true }
+serde = { version = "1.0", features = ["derive"], optional = true }
+serde_json = { version = "1.0", optional = true }
 
 [dev-dependencies]
 criterion = { version = "0.8", features = ["html_reports"] }
@@ -28,3 +30,4 @@ harness = false
 [features]
 wide = ["dep:wide"]
 simd = ["wide"]
+serde = ["dep:serde", "dep:serde_json"]

crates/ferro_ta_core/README.md

@@ -13,7 +13,7 @@ PyO3, NumPy, or Python runtime dependency, which makes it a good fit for:
 
 ```toml
 [dependencies]
-ferro_ta_core = "1.0.6"
+ferro_ta_core = "1.2.0"
 ```
 
 ## Design

crates/ferro_ta_core/src/aggregation.rs

@@ -0,0 +1,346 @@
+//! Tick / Trade Aggregation Pipeline — pure Rust, no PyO3.
+//!
+//! Aggregates raw tick/trade data into OHLCV bars:
+//! - **tick bars**   — fixed number of ticks per bar
+//! - **volume bars** — fixed volume threshold per bar
+//! - **time bars**   — label-based grouping (labels from Python timestamps)
+
+/// OHLCV 5-tuple return type alias.
+type Ohlcv5 = (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>);
+
+/// OHLCV 5-tuple plus labels return type alias.
+type Ohlcv5AndLabels = (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<i64>);
+
+// ---------------------------------------------------------------------------
+// aggregate_tick_bars
+// ---------------------------------------------------------------------------
+
+/// Aggregate tick/trade data into tick bars (every N ticks become one bar).
+///
+/// Returns `(open, high, low, close, volume)` where volume = sum of sizes.
+///
+/// # Panics
+/// Panics if `ticks_per_bar == 0`, arrays are empty, or lengths differ.
+pub fn aggregate_tick_bars(
+    price: &[f64],
+    size: &[f64],
+    ticks_per_bar: usize,
+) -> Ohlcv5 {
+    assert!(ticks_per_bar >= 1, "ticks_per_bar must be >= 1");
+    let n = price.len();
+    assert!(n > 0 && size.len() == n, "price and size must be non-empty and equal length");
+
+    let n_bars = n.div_ceil(ticks_per_bar);
+    let mut out_open = Vec::with_capacity(n_bars);
+    let mut out_high = Vec::with_capacity(n_bars);
+    let mut out_low = Vec::with_capacity(n_bars);
+    let mut out_close = Vec::with_capacity(n_bars);
+    let mut out_vol = Vec::with_capacity(n_bars);
+
+    let mut i = 0;
+    while i < n {
+        let end = (i + ticks_per_bar).min(n);
+        let bar_p = &price[i..end];
+        let bar_s = &size[i..end];
+        let bar_open = bar_p[0];
+        let bar_high = bar_p.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+        let bar_low = bar_p.iter().cloned().fold(f64::INFINITY, f64::min);
+        let bar_close = *bar_p.last().expect("slice cannot be empty");
+        let bar_vol: f64 = bar_s.iter().sum();
+        out_open.push(bar_open);
+        out_high.push(bar_high);
+        out_low.push(bar_low);
+        out_close.push(bar_close);
+        out_vol.push(bar_vol);
+        i = end;
+    }
+
+    (out_open, out_high, out_low, out_close, out_vol)
+}
+
+// ---------------------------------------------------------------------------
+// aggregate_volume_bars_ticks
+// ---------------------------------------------------------------------------
+
+/// Aggregate tick data into volume bars (fixed volume threshold).
+///
+/// Accumulates ticks until cumulative size >= `volume_threshold`, then emits
+/// a bar. Any remaining partial bar is also emitted.
+///
+/// Returns `(open, high, low, close, volume)`.
+///
+/// # Panics
+/// Panics if `volume_threshold <= 0`, arrays are empty, or lengths differ.
+pub fn aggregate_volume_bars_ticks(
+    price: &[f64],
+    size: &[f64],
+    volume_threshold: f64,
+) -> Ohlcv5 {
+    assert!(volume_threshold > 0.0, "volume_threshold must be > 0");
+    let n = price.len();
+    assert!(n > 0 && size.len() == n, "price and size must be non-empty and equal length");
+
+    let mut out_open: Vec<f64> = Vec::new();
+    let mut out_high: Vec<f64> = Vec::new();
+    let mut out_low: Vec<f64> = Vec::new();
+    let mut out_close: Vec<f64> = Vec::new();
+    let mut out_vol: Vec<f64> = Vec::new();
+
+    let mut bar_open = price[0];
+    let mut bar_high = price[0];
+    let mut bar_low = price[0];
+    let mut bar_close = price[0];
+    let mut bar_vol = size[0];
+
+    for i in 1..n {
+        bar_high = bar_high.max(price[i]);
+        bar_low = bar_low.min(price[i]);
+        bar_close = price[i];
+        bar_vol += size[i];
+
+        if bar_vol >= volume_threshold {
+            out_open.push(bar_open);
+            out_high.push(bar_high);
+            out_low.push(bar_low);
+            out_close.push(bar_close);
+            out_vol.push(bar_vol);
+            if i + 1 < n {
+                bar_open = price[i + 1];
+                bar_high = price[i + 1];
+                bar_low = price[i + 1];
+                bar_close = price[i + 1];
+                bar_vol = size[i + 1];
+            } else {
+                bar_vol = 0.0;
+            }
+        }
+    }
+    // Push remaining partial bar
+    if bar_vol > 0.0 {
+        out_open.push(bar_open);
+        out_high.push(bar_high);
+        out_low.push(bar_low);
+        out_close.push(bar_close);
+        out_vol.push(bar_vol);
+    }
+
+    (out_open, out_high, out_low, out_close, out_vol)
+}
+
+// ---------------------------------------------------------------------------
+// aggregate_time_bars
+// ---------------------------------------------------------------------------
+
+/// Aggregate tick data into time bars using pre-computed integer bucket labels.
+///
+/// Each tick is assigned a `label` (e.g. unix_ts // period_secs). Ticks with
+/// the same label are accumulated into one bar. Labels must be non-decreasing.
+///
+/// Returns `(open, high, low, close, volume, unique_labels)`.
+///
+/// # Panics
+/// Panics if arrays are empty or have unequal lengths.
+pub fn aggregate_time_bars(
+    price: &[f64],
+    size: &[f64],
+    labels: &[i64],
+) -> Ohlcv5AndLabels {
+    let n = price.len();
+    assert!(
+        n > 0 && size.len() == n && labels.len() == n,
+        "price, size, and labels must be non-empty and equal length"
+    );
+
+    let mut out_open: Vec<f64> = Vec::new();
+    let mut out_high: Vec<f64> = Vec::new();
+    let mut out_low: Vec<f64> = Vec::new();
+    let mut out_close: Vec<f64> = Vec::new();
+    let mut out_vol: Vec<f64> = Vec::new();
+    let mut out_labels: Vec<i64> = Vec::new();
+
+    let mut cur_label = labels[0];
+    let mut bar_open = price[0];
+    let mut bar_high = price[0];
+    let mut bar_low = price[0];
+    let mut bar_close = price[0];
+    let mut bar_vol = size[0];
+
+    for i in 1..n {
+        if labels[i] != cur_label {
+            out_open.push(bar_open);
+            out_high.push(bar_high);
+            out_low.push(bar_low);
+            out_close.push(bar_close);
+            out_vol.push(bar_vol);
+            out_labels.push(cur_label);
+            cur_label = labels[i];
+            bar_open = price[i];
+            bar_high = price[i];
+            bar_low = price[i];
+            bar_close = price[i];
+            bar_vol = size[i];
+        } else {
+            bar_high = bar_high.max(price[i]);
+            bar_low = bar_low.min(price[i]);
+            bar_close = price[i];
+            bar_vol += size[i];
+        }
+    }
+    out_open.push(bar_open);
+    out_high.push(bar_high);
+    out_low.push(bar_low);
+    out_close.push(bar_close);
+    out_vol.push(bar_vol);
+    out_labels.push(cur_label);
+
+    (out_open, out_high, out_low, out_close, out_vol, out_labels)
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // -- aggregate_tick_bars -------------------------------------------------
+
+    #[test]
+    fn test_tick_bars_exact_division() {
+        let price = [10.0, 11.0, 12.0, 13.0, 14.0, 15.0];
+        let size = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
+        let (o, h, l, c, v) = aggregate_tick_bars(&price, &size, 3);
+        assert_eq!(o.len(), 2);
+        // Bar 0: ticks 0..3
+        assert!((o[0] - 10.0).abs() < 1e-10);
+        assert!((h[0] - 12.0).abs() < 1e-10);
+        assert!((l[0] - 10.0).abs() < 1e-10);
+        assert!((c[0] - 12.0).abs() < 1e-10);
+        assert!((v[0] - 6.0).abs() < 1e-10);
+        // Bar 1: ticks 3..6
+        assert!((o[1] - 13.0).abs() < 1e-10);
+        assert!((h[1] - 15.0).abs() < 1e-10);
+        assert!((l[1] - 13.0).abs() < 1e-10);
+        assert!((c[1] - 15.0).abs() < 1e-10);
+        assert!((v[1] - 15.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_tick_bars_partial_last_bar() {
+        let price = [10.0, 11.0, 12.0, 13.0, 14.0];
+        let size = [1.0, 2.0, 3.0, 4.0, 5.0];
+        let (o, _h, _l, c, v) = aggregate_tick_bars(&price, &size, 3);
+        assert_eq!(o.len(), 2);
+        // Partial bar: ticks 3..5
+        assert!((o[1] - 13.0).abs() < 1e-10);
+        assert!((c[1] - 14.0).abs() < 1e-10);
+        assert!((v[1] - 9.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_tick_bars_single_tick() {
+        let (o, h, l, c, v) = aggregate_tick_bars(&[42.0], &[100.0], 5);
+        assert_eq!(o.len(), 1);
+        assert!((o[0] - 42.0).abs() < 1e-10);
+        assert!((h[0] - 42.0).abs() < 1e-10);
+        assert!((l[0] - 42.0).abs() < 1e-10);
+        assert!((c[0] - 42.0).abs() < 1e-10);
+        assert!((v[0] - 100.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "ticks_per_bar must be >= 1")]
+    fn test_tick_bars_zero_ticks() {
+        aggregate_tick_bars(&[1.0], &[1.0], 0);
+    }
+
+    // -- aggregate_volume_bars_ticks -----------------------------------------
+
+    #[test]
+    fn test_volume_bars_ticks_basic() {
+        let price = [10.0, 11.0, 12.0, 13.0, 14.0];
+        let size = [30.0, 40.0, 50.0, 20.0, 60.0];
+        // threshold=70: bar0 = ticks 0+1 (vol=70), bar1 = tick2 (vol=50) + tick3 (vol=70),
+        // then tick4 as partial
+        let (o, h, l, c, v) = aggregate_volume_bars_ticks(&price, &size, 70.0);
+        // First bar: 30+40=70 >= 70
+        assert!((o[0] - 10.0).abs() < 1e-10);
+        assert!((c[0] - 11.0).abs() < 1e-10);
+        assert!((v[0] - 70.0).abs() < 1e-10);
+        assert!((h[0] - 11.0).abs() < 1e-10);
+        assert!((l[0] - 10.0).abs() < 1e-10);
+        assert!(v.len() >= 2);
+    }
+
+    #[test]
+    fn test_volume_bars_ticks_single() {
+        let (o, _h, _l, _c, v) = aggregate_volume_bars_ticks(&[5.0], &[10.0], 100.0);
+        assert_eq!(o.len(), 1);
+        assert!((v[0] - 10.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "volume_threshold must be > 0")]
+    fn test_volume_bars_ticks_zero_threshold() {
+        aggregate_volume_bars_ticks(&[1.0], &[1.0], 0.0);
+    }
+
+    // -- aggregate_time_bars -------------------------------------------------
+
+    #[test]
+    fn test_time_bars_basic() {
+        let price = [10.0, 11.0, 12.0, 13.0, 14.0];
+        let size = [1.0, 2.0, 3.0, 4.0, 5.0];
+        let labels: [i64; 5] = [0, 0, 1, 1, 1];
+        let (o, h, l, c, v, out_lbl) = aggregate_time_bars(&price, &size, &labels);
+        assert_eq!(o.len(), 2);
+        assert_eq!(out_lbl, vec![0, 1]);
+        // Group 0: ticks 0,1
+        assert!((o[0] - 10.0).abs() < 1e-10);
+        assert!((h[0] - 11.0).abs() < 1e-10);
+        assert!((l[0] - 10.0).abs() < 1e-10);
+        assert!((c[0] - 11.0).abs() < 1e-10);
+        assert!((v[0] - 3.0).abs() < 1e-10);
+        // Group 1: ticks 2,3,4
+        assert!((o[1] - 12.0).abs() < 1e-10);
+        assert!((h[1] - 14.0).abs() < 1e-10);
+        assert!((l[1] - 12.0).abs() < 1e-10);
+        assert!((c[1] - 14.0).abs() < 1e-10);
+        assert!((v[1] - 12.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_time_bars_all_same_label() {
+        let price = [5.0, 6.0, 4.0];
+        let size = [10.0, 20.0, 30.0];
+        let labels: [i64; 3] = [42, 42, 42];
+        let (o, h, l, c, v, out_lbl) = aggregate_time_bars(&price, &size, &labels);
+        assert_eq!(o.len(), 1);
+        assert_eq!(out_lbl, vec![42]);
+        assert!((o[0] - 5.0).abs() < 1e-10);
+        assert!((h[0] - 6.0).abs() < 1e-10);
+        assert!((l[0] - 4.0).abs() < 1e-10);
+        assert!((c[0] - 4.0).abs() < 1e-10);
+        assert!((v[0] - 60.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_time_bars_each_tick_own_label() {
+        let price = [10.0, 20.0, 30.0];
+        let size = [1.0, 2.0, 3.0];
+        let labels: [i64; 3] = [0, 1, 2];
+        let (o, _h, _l, _c, v, out_lbl) = aggregate_time_bars(&price, &size, &labels);
+        assert_eq!(o.len(), 3);
+        assert_eq!(out_lbl, vec![0, 1, 2]);
+        assert!((v[0] - 1.0).abs() < 1e-10);
+        assert!((v[1] - 2.0).abs() < 1e-10);
+        assert!((v[2] - 3.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "price, size, and labels must be non-empty and equal length")]
+    fn test_time_bars_empty() {
+        aggregate_time_bars(&[], &[], &[]);
+    }
+}

crates/ferro_ta_core/src/alerts.rs

@@ -0,0 +1,130 @@
+//! Alerts — condition evaluation helpers.
+//!
+//! - `check_threshold` — fires when a series crosses above/below a level
+//! - `check_cross`     — fires when *fast* crosses above or below *slow*
+//! - `collect_alert_bars` — returns indices of bars where a mask is non-zero
+
+/// Fire an alert when `series` crosses a threshold level.
+///
+/// `direction`: `1` = cross above, `-1` = cross below.
+///
+/// Returns a `Vec<i8>` with `1` at crossing bars, `0` elsewhere.
+/// Element 0 is always 0.
+pub fn check_threshold(series: &[f64], level: f64, direction: i32) -> Vec<i8> {
+    let n = series.len();
+    let mut out = vec![0i8; n];
+    if n < 2 {
+        return out;
+    }
+    for i in 1..n {
+        let prev = series[i - 1];
+        let curr = series[i];
+        if prev.is_nan() || curr.is_nan() {
+            continue;
+        }
+        if direction == 1 {
+            if prev <= level && curr > level {
+                out[i] = 1;
+            }
+        } else if direction == -1 {
+            if prev >= level && curr < level {
+                out[i] = 1;
+            }
+        }
+    }
+    out
+}
+
+/// Detect cross-over / cross-under events between two series.
+///
+/// Returns `Vec<i8>`: `1` = bullish cross (fast above slow), `-1` = bearish, `0` = none.
+/// Element 0 is always 0.
+pub fn check_cross(fast: &[f64], slow: &[f64]) -> Vec<i8> {
+    let n = fast.len();
+    let mut out = vec![0i8; n];
+    if n < 2 {
+        return out;
+    }
+    for i in 1..n {
+        let fp = fast[i - 1];
+        let fc = fast[i];
+        let sp = slow[i - 1];
+        let sc = slow[i];
+        if fp.is_nan() || fc.is_nan() || sp.is_nan() || sc.is_nan() {
+            continue;
+        }
+        if fp <= sp && fc > sc {
+            out[i] = 1;
+        } else if fp >= sp && fc < sc {
+            out[i] = -1;
+        }
+    }
+    out
+}
+
+/// Collect bar indices where `mask` is non-zero.
+pub fn collect_alert_bars(mask: &[i8]) -> Vec<i64> {
+    mask.iter()
+        .enumerate()
+        .filter(|(_, &v)| v != 0)
+        .map(|(i, _)| i as i64)
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_check_threshold_cross_above() {
+        let series = vec![10.0, 20.0, 30.0, 40.0, 50.0];
+        let result = check_threshold(&series, 25.0, 1);
+        assert_eq!(result, vec![0, 0, 1, 0, 0]);
+    }
+
+    #[test]
+    fn test_check_threshold_cross_below() {
+        let series = vec![50.0, 40.0, 30.0, 20.0, 10.0];
+        let result = check_threshold(&series, 25.0, -1);
+        assert_eq!(result, vec![0, 0, 0, 1, 0]);
+    }
+
+    #[test]
+    fn test_check_cross_bullish() {
+        let fast = vec![1.0, 2.0, 5.0];
+        let slow = vec![3.0, 3.0, 3.0];
+        let result = check_cross(&fast, &slow);
+        assert_eq!(result, vec![0, 0, 1]);
+    }
+
+    #[test]
+    fn test_check_cross_bearish() {
+        let fast = vec![5.0, 4.0, 1.0];
+        let slow = vec![3.0, 3.0, 3.0];
+        let result = check_cross(&fast, &slow);
+        assert_eq!(result, vec![0, 0, -1]);
+    }
+
+    #[test]
+    fn test_collect_alert_bars() {
+        let mask = vec![0i8, 1, 0, -1, 0, 1];
+        let result = collect_alert_bars(&mask);
+        assert_eq!(result, vec![1, 3, 5]);
+    }
+
+    #[test]
+    fn test_empty() {
+        assert_eq!(check_threshold(&[], 0.0, 1), Vec::<i8>::new());
+        assert_eq!(check_cross(&[], &[]), Vec::<i8>::new());
+        assert_eq!(collect_alert_bars(&[]), Vec::<i64>::new());
+    }
+
+    #[test]
+    fn test_nan_handling() {
+        let series = vec![10.0, f64::NAN, 30.0, 40.0];
+        let result = check_threshold(&series, 25.0, 1);
+        // NaN bars are skipped
+        assert_eq!(result[1], 0);
+        assert_eq!(result[2], 0); // prev is NaN
+    }
+}

crates/ferro_ta_core/src/attribution.rs

@@ -0,0 +1,329 @@
+//! Performance attribution and trade analysis — pure Rust, no PyO3.
+//!
+//! Functions
+//! ---------
+//! - `trade_stats`          — win rate, avg win/loss, profit factor, avg hold
+//! - `monthly_contribution` — group bar returns by month index and sum
+//! - `signal_attribution`   — group bar returns by signal label and sum
+//! - `extract_trades`       — extract trade pnl and hold durations from positions
+
+use std::collections::HashMap;
+
+// ---------------------------------------------------------------------------
+// trade_stats
+// ---------------------------------------------------------------------------
+
+/// Compute trade-level statistics from trade PnL and hold durations.
+///
+/// Returns `(win_rate, avg_win, avg_loss, profit_factor, avg_hold_bars)`.
+///
+/// - **win_rate**      : fraction of trades with PnL > 0
+/// - **avg_win**       : mean PnL of winning trades (0 if none)
+/// - **avg_loss**      : mean PnL of losing trades (negative; 0 if none)
+/// - **profit_factor** : gross profit / |gross loss| (inf if no losses)
+/// - **avg_hold_bars** : mean hold duration across all trades
+///
+/// # Panics
+/// Panics if `pnl` is empty or `pnl.len() != hold_bars.len()`.
+pub fn trade_stats(pnl: &[f64], hold_bars: &[f64]) -> (f64, f64, f64, f64, f64) {
+    let n = pnl.len();
+    assert!(n > 0, "pnl must be non-empty");
+    assert_eq!(n, hold_bars.len(), "pnl and hold_bars must have equal length");
+
+    let mut wins: Vec<f64> = Vec::new();
+    let mut losses: Vec<f64> = Vec::new();
+    for &v in pnl.iter() {
+        if v > 0.0 {
+            wins.push(v);
+        } else if v < 0.0 {
+            losses.push(v);
+        }
+    }
+
+    let win_rate = wins.len() as f64 / n as f64;
+    let avg_win = if wins.is_empty() {
+        0.0
+    } else {
+        wins.iter().sum::<f64>() / wins.len() as f64
+    };
+    let avg_loss = if losses.is_empty() {
+        0.0
+    } else {
+        losses.iter().sum::<f64>() / losses.len() as f64
+    };
+
+    let gross_profit: f64 = wins.iter().sum();
+    let gross_loss: f64 = losses.iter().map(|v| v.abs()).sum();
+    let profit_factor = if gross_loss == 0.0 {
+        f64::INFINITY
+    } else {
+        gross_profit / gross_loss
+    };
+
+    let avg_hold = hold_bars.iter().sum::<f64>() / n as f64;
+
+    (win_rate, avg_win, avg_loss, profit_factor, avg_hold)
+}
+
+// ---------------------------------------------------------------------------
+// monthly_contribution
+// ---------------------------------------------------------------------------
+
+/// Group per-bar returns by month index and sum each month's contribution.
+///
+/// Returns `(months, contributions)` where `months` is sorted unique month
+/// indices and `contributions` is the corresponding total return per month.
+/// NaN returns are skipped.
+///
+/// # Panics
+/// Panics if `bar_returns.len() != month_index.len()`.
+pub fn monthly_contribution(bar_returns: &[f64], month_index: &[i64]) -> (Vec<i64>, Vec<f64>) {
+    let n = bar_returns.len();
+    assert_eq!(
+        n,
+        month_index.len(),
+        "bar_returns and month_index must have equal length"
+    );
+
+    let mut map: HashMap<i64, f64> = HashMap::new();
+    for i in 0..n {
+        if !bar_returns[i].is_nan() {
+            *map.entry(month_index[i]).or_insert(0.0) += bar_returns[i];
+        }
+    }
+
+    let mut months: Vec<i64> = map.keys().copied().collect();
+    months.sort_unstable();
+    let contributions: Vec<f64> = months.iter().map(|m| map[m]).collect();
+
+    (months, contributions)
+}
+
+// ---------------------------------------------------------------------------
+// signal_attribution
+// ---------------------------------------------------------------------------
+
+/// Attribute per-bar returns to each signal label.
+///
+/// Returns `(labels, contributions)` where `labels` is sorted unique signal
+/// labels and `contributions` is the corresponding total return per label.
+/// NaN returns are skipped.
+///
+/// # Panics
+/// Panics if `bar_returns.len() != signal_labels.len()`.
+pub fn signal_attribution(bar_returns: &[f64], signal_labels: &[i64]) -> (Vec<i64>, Vec<f64>) {
+    let n = bar_returns.len();
+    assert_eq!(
+        n,
+        signal_labels.len(),
+        "bar_returns and signal_labels must have equal length"
+    );
+
+    let mut map: HashMap<i64, f64> = HashMap::new();
+    for i in 0..n {
+        if !bar_returns[i].is_nan() {
+            *map.entry(signal_labels[i]).or_insert(0.0) += bar_returns[i];
+        }
+    }
+
+    let mut labels: Vec<i64> = map.keys().copied().collect();
+    labels.sort_unstable();
+    let contributions: Vec<f64> = labels.iter().map(|l| map[l]).collect();
+
+    (labels, contributions)
+}
+
+// ---------------------------------------------------------------------------
+// extract_trades
+// ---------------------------------------------------------------------------
+
+/// Extract trade-level PnL and hold durations from positions and strategy returns.
+///
+/// A trade is a maximal contiguous run of non-zero position values with the
+/// same sign/magnitude. Returns `(pnl, hold_durations)`.
+///
+/// # Panics
+/// Panics if `positions.len() != strategy_returns.len()`.
+pub fn extract_trades(positions: &[f64], strategy_returns: &[f64]) -> (Vec<f64>, Vec<f64>) {
+    let n = positions.len();
+    assert_eq!(
+        n,
+        strategy_returns.len(),
+        "positions and strategy_returns must have equal length"
+    );
+
+    let mut pnl = Vec::<f64>::new();
+    let mut hold = Vec::<f64>::new();
+
+    let mut i = 0usize;
+    while i < n {
+        if positions[i] == 0.0 {
+            i += 1;
+            continue;
+        }
+        let mut j = i + 1;
+        while j < n && positions[j] == positions[i] {
+            j += 1;
+        }
+        let mut trade_pnl = 0.0_f64;
+        for v in strategy_returns.iter().take(j).skip(i) {
+            trade_pnl += *v;
+        }
+        pnl.push(trade_pnl);
+        hold.push((j - i) as f64);
+        i = j;
+    }
+
+    (pnl, hold)
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // -- trade_stats ---------------------------------------------------------
+
+    #[test]
+    fn test_trade_stats_basic() {
+        let pnl = [100.0, -50.0, 200.0, -30.0, 150.0];
+        let hold = [5.0, 3.0, 7.0, 2.0, 6.0];
+        let (wr, aw, al, pf, ah) = trade_stats(&pnl, &hold);
+
+        // 3 wins out of 5
+        assert!((wr - 0.6).abs() < 1e-10);
+        // avg win = (100+200+150)/3
+        assert!((aw - 150.0).abs() < 1e-10);
+        // avg loss = (-50 + -30)/2 = -40
+        assert!((al - (-40.0)).abs() < 1e-10);
+        // profit_factor = 450 / 80
+        assert!((pf - 5.625).abs() < 1e-10);
+        // avg hold = (5+3+7+2+6)/5 = 4.6
+        assert!((ah - 4.6).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_trade_stats_all_wins() {
+        let pnl = [10.0, 20.0];
+        let hold = [1.0, 2.0];
+        let (wr, _aw, al, pf, _ah) = trade_stats(&pnl, &hold);
+        assert!((wr - 1.0).abs() < 1e-10);
+        assert!((al - 0.0).abs() < 1e-10);
+        assert!(pf.is_infinite());
+    }
+
+    #[test]
+    fn test_trade_stats_all_losses() {
+        let pnl = [-10.0, -20.0];
+        let hold = [1.0, 2.0];
+        let (wr, aw, _al, pf, _ah) = trade_stats(&pnl, &hold);
+        assert!((wr - 0.0).abs() < 1e-10);
+        assert!((aw - 0.0).abs() < 1e-10);
+        assert!((pf - 0.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "pnl must be non-empty")]
+    fn test_trade_stats_empty() {
+        trade_stats(&[], &[]);
+    }
+
+    // -- monthly_contribution ------------------------------------------------
+
+    #[test]
+    fn test_monthly_contribution_basic() {
+        let returns = [0.01, 0.02, -0.01, 0.03, -0.02];
+        let months = [0, 0, 1, 1, 2];
+        let (m, c) = monthly_contribution(&returns, &months);
+        assert_eq!(m, vec![0, 1, 2]);
+        assert!((c[0] - 0.03).abs() < 1e-10);
+        assert!((c[1] - 0.02).abs() < 1e-10);
+        assert!((c[2] - (-0.02)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_monthly_contribution_nan_skipped() {
+        let returns = [0.01, f64::NAN, 0.03];
+        let months = [0, 0, 1];
+        let (m, c) = monthly_contribution(&returns, &months);
+        assert_eq!(m, vec![0, 1]);
+        assert!((c[0] - 0.01).abs() < 1e-10);
+        assert!((c[1] - 0.03).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_monthly_contribution_empty() {
+        let (m, c) = monthly_contribution(&[], &[]);
+        assert!(m.is_empty());
+        assert!(c.is_empty());
+    }
+
+    // -- signal_attribution --------------------------------------------------
+
+    #[test]
+    fn test_signal_attribution_basic() {
+        let returns = [0.05, -0.02, 0.03, 0.01];
+        let labels = [1, -1, 2, 1];
+        let (l, c) = signal_attribution(&returns, &labels);
+        assert_eq!(l, vec![-1, 1, 2]);
+        assert!((c[0] - (-0.02)).abs() < 1e-10);
+        assert!((c[1] - 0.06).abs() < 1e-10); // 0.05 + 0.01
+        assert!((c[2] - 0.03).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_signal_attribution_nan_skipped() {
+        let returns = [0.05, f64::NAN];
+        let labels = [1, 2];
+        let (l, c) = signal_attribution(&returns, &labels);
+        assert_eq!(l, vec![1]);
+        assert!((c[0] - 0.05).abs() < 1e-10);
+    }
+
+    // -- extract_trades ------------------------------------------------------
+
+    #[test]
+    fn test_extract_trades_basic() {
+        // positions: flat, long, long, flat, short, short
+        let positions = [0.0, 1.0, 1.0, 0.0, -1.0, -1.0];
+        let strat_ret = [0.0, 0.01, 0.02, 0.0, -0.01, 0.03];
+        let (pnl, hold) = extract_trades(&positions, &strat_ret);
+        assert_eq!(pnl.len(), 2);
+        assert_eq!(hold.len(), 2);
+        // First trade: bars 1..3 => 0.01 + 0.02 = 0.03
+        assert!((pnl[0] - 0.03).abs() < 1e-10);
+        assert!((hold[0] - 2.0).abs() < 1e-10);
+        // Second trade: bars 4..6 => -0.01 + 0.03 = 0.02
+        assert!((pnl[1] - 0.02).abs() < 1e-10);
+        assert!((hold[1] - 2.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_extract_trades_all_flat() {
+        let positions = [0.0, 0.0, 0.0];
+        let strat_ret = [0.01, 0.02, 0.03];
+        let (pnl, hold) = extract_trades(&positions, &strat_ret);
+        assert!(pnl.is_empty());
+        assert!(hold.is_empty());
+    }
+
+    #[test]
+    fn test_extract_trades_empty() {
+        let (pnl, hold) = extract_trades(&[], &[]);
+        assert!(pnl.is_empty());
+        assert!(hold.is_empty());
+    }
+
+    #[test]
+    fn test_extract_trades_single_bar_trade() {
+        let positions = [0.0, 1.0, 0.0];
+        let strat_ret = [0.0, 0.05, 0.0];
+        let (pnl, hold) = extract_trades(&positions, &strat_ret);
+        assert_eq!(pnl.len(), 1);
+        assert!((pnl[0] - 0.05).abs() < 1e-10);
+        assert!((hold[0] - 1.0).abs() < 1e-10);
+    }
+}

crates/ferro_ta_core/src/batch.rs

@@ -0,0 +1,642 @@
+//! Pure-Rust batch operations — apply indicators across multiple series
+//! (columns) sequentially. The PyO3 wrapper can add Rayon parallelism on top.
+//!
+//! Input convention: `data[j]` is column *j* (one time-series). All columns
+//! must have the same length.
+
+use crate::{momentum, overlap, statistic, volatility};
+
+// ---------------------------------------------------------------------------
+// helpers
+// ---------------------------------------------------------------------------
+
+/// Validate that every column in `data` has the same length.  Returns `Ok(n)`
+/// where `n` is the common length, or `Err` with a message.
+fn validate_columns(data: &[Vec<f64>]) -> Result<usize, String> {
+    if data.is_empty() {
+        return Ok(0);
+    }
+    let n = data[0].len();
+    for (idx, col) in data.iter().enumerate() {
+        if col.len() != n {
+            return Err(format!(
+                "column 0 has length {n}, but column {idx} has length {}",
+                col.len()
+            ));
+        }
+    }
+    Ok(n)
+}
+
+fn validate_hlc_columns(
+    high: &[Vec<f64>],
+    low: &[Vec<f64>],
+    close: &[Vec<f64>],
+) -> Result<(usize, usize), String> {
+    let n_series = high.len();
+    if low.len() != n_series || close.len() != n_series {
+        return Err(format!(
+            "high has {} columns, low has {}, close has {} — must be equal",
+            n_series,
+            low.len(),
+            close.len()
+        ));
+    }
+    if n_series == 0 {
+        return Ok((0, 0));
+    }
+    let n = high[0].len();
+    for (idx, (h, (l, c))) in high
+        .iter()
+        .zip(low.iter().zip(close.iter()))
+        .enumerate()
+    {
+        if h.len() != n || l.len() != n || c.len() != n {
+            return Err(format!(
+                "column {idx}: high len={}, low len={}, close len={} — must all be {n}",
+                h.len(),
+                l.len(),
+                c.len()
+            ));
+        }
+    }
+    Ok((n, n_series))
+}
+
+// ---------------------------------------------------------------------------
+// rolling linear regression (self-contained so core has no PyO3 dep)
+// ---------------------------------------------------------------------------
+
+fn linreg(window: &[f64]) -> (f64, f64) {
+    let n = window.len() as f64;
+    let sum_x: f64 = (0..window.len()).map(|i| i as f64).sum();
+    let sum_y: f64 = window.iter().sum();
+    let sum_xy: f64 = window.iter().enumerate().map(|(i, &y)| i as f64 * y).sum();
+    let sum_x2: f64 = (0..window.len()).map(|i| (i as f64).powi(2)).sum();
+    let denom = n * sum_x2 - sum_x * sum_x;
+    let slope = if denom != 0.0 {
+        (n * sum_xy - sum_x * sum_y) / denom
+    } else {
+        0.0
+    };
+    let intercept = (sum_y - slope * sum_x) / n;
+    (slope, intercept)
+}
+
+fn rolling_linreg_apply<F>(prices: &[f64], timeperiod: usize, mut map: F) -> Vec<f64>
+where
+    F: FnMut(f64, f64) -> f64,
+{
+    let n = prices.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+
+    if prices.iter().any(|value| !value.is_finite()) {
+        for end in (timeperiod - 1)..n {
+            let window = &prices[(end + 1 - timeperiod)..=end];
+            let (slope, intercept) = linreg(window);
+            result[end] = map(slope, intercept);
+        }
+        return result;
+    }
+
+    let period = timeperiod as f64;
+    let last_x = (timeperiod - 1) as f64;
+    let sum_x = last_x * period / 2.0;
+    let sum_x2 = last_x * period * (2.0 * period - 1.0) / 6.0;
+    let denom = period * sum_x2 - sum_x * sum_x;
+
+    let mut sum_y = prices[..timeperiod].iter().sum::<f64>();
+    let mut sum_xy = prices[..timeperiod]
+        .iter()
+        .enumerate()
+        .map(|(idx, &value)| idx as f64 * value)
+        .sum::<f64>();
+
+    for end in (timeperiod - 1)..n {
+        let slope = if denom != 0.0 {
+            (period * sum_xy - sum_x * sum_y) / denom
+        } else {
+            0.0
+        };
+        let intercept = (sum_y - slope * sum_x) / period;
+        result[end] = map(slope, intercept);
+
+        if end + 1 < n {
+            let outgoing = prices[end + 1 - timeperiod];
+            let incoming = prices[end + 1];
+            let prev_sum_y = sum_y;
+
+            sum_y = prev_sum_y - outgoing + incoming;
+            sum_xy = sum_xy - (prev_sum_y - outgoing) + last_x * incoming;
+        }
+    }
+
+    result
+}
+
+// ---------------------------------------------------------------------------
+// CCI / WILLR helpers (no external dep)
+// ---------------------------------------------------------------------------
+
+fn compute_cci(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = high.len();
+    let typical_price: Vec<f64> = high
+        .iter()
+        .zip(low.iter())
+        .zip(close.iter())
+        .map(|((&h, &l), &c)| (h + l + c) / 3.0)
+        .collect();
+
+    let mut result = vec![f64::NAN; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+    for end in (timeperiod - 1)..n {
+        let window = &typical_price[(end + 1 - timeperiod)..=end];
+        let mean = window.iter().sum::<f64>() / timeperiod as f64;
+        let mad = window
+            .iter()
+            .map(|&value| (value - mean).abs())
+            .sum::<f64>()
+            / timeperiod as f64;
+        result[end] = if mad != 0.0 {
+            (typical_price[end] - mean) / (0.015 * mad)
+        } else {
+            0.0
+        };
+    }
+    result
+}
+
+fn compute_willr(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = high.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+
+    // Use simple sliding-window max/min
+    for end in (timeperiod - 1)..n {
+        let start = end + 1 - timeperiod;
+        let mut highest = f64::NEG_INFINITY;
+        let mut lowest = f64::INFINITY;
+        for i in start..=end {
+            if high[i] > highest {
+                highest = high[i];
+            }
+            if low[i] < lowest {
+                lowest = low[i];
+            }
+        }
+        let range = highest - lowest;
+        result[end] = if range != 0.0 {
+            -100.0 * (highest - close[end]) / range
+        } else {
+            -50.0
+        };
+    }
+
+    result
+}
+
+// ---------------------------------------------------------------------------
+// batch_sma
+// ---------------------------------------------------------------------------
+
+/// Apply SMA to each column. Returns one output column per input column.
+pub fn batch_sma(data: &[Vec<f64>], timeperiod: usize) -> Result<Vec<Vec<f64>>, String> {
+    if timeperiod == 0 {
+        return Err("timeperiod must be >= 1".into());
+    }
+    validate_columns(data)?;
+    Ok(data
+        .iter()
+        .map(|col| overlap::sma(col, timeperiod))
+        .collect())
+}
+
+// ---------------------------------------------------------------------------
+// batch_ema
+// ---------------------------------------------------------------------------
+
+/// Apply EMA to each column.
+pub fn batch_ema(data: &[Vec<f64>], timeperiod: usize) -> Result<Vec<Vec<f64>>, String> {
+    if timeperiod == 0 {
+        return Err("timeperiod must be >= 1".into());
+    }
+    validate_columns(data)?;
+    Ok(data
+        .iter()
+        .map(|col| overlap::ema(col, timeperiod))
+        .collect())
+}
+
+// ---------------------------------------------------------------------------
+// batch_rsi
+// ---------------------------------------------------------------------------
+
+/// Apply RSI to each column.
+pub fn batch_rsi(data: &[Vec<f64>], timeperiod: usize) -> Result<Vec<Vec<f64>>, String> {
+    if timeperiod == 0 {
+        return Err("timeperiod must be >= 1".into());
+    }
+    validate_columns(data)?;
+    Ok(data
+        .iter()
+        .map(|col| momentum::rsi(col, timeperiod))
+        .collect())
+}
+
+// ---------------------------------------------------------------------------
+// batch_atr
+// ---------------------------------------------------------------------------
+
+/// Apply ATR to each set of (high, low, close) columns.
+pub fn batch_atr(
+    high: &[Vec<f64>],
+    low: &[Vec<f64>],
+    close: &[Vec<f64>],
+    timeperiod: usize,
+) -> Result<Vec<Vec<f64>>, String> {
+    if timeperiod == 0 {
+        return Err("timeperiod must be >= 1".into());
+    }
+    validate_hlc_columns(high, low, close)?;
+    Ok((0..high.len())
+        .map(|i| volatility::atr(&high[i], &low[i], &close[i], timeperiod))
+        .collect())
+}
+
+// ---------------------------------------------------------------------------
+// batch_stoch
+// ---------------------------------------------------------------------------
+
+/// Apply Stochastic to each set of (high, low, close) columns.
+/// Returns `(slowk_columns, slowd_columns)`.
+pub fn batch_stoch(
+    high: &[Vec<f64>],
+    low: &[Vec<f64>],
+    close: &[Vec<f64>],
+    fastk_period: usize,
+    slowk_period: usize,
+    slowd_period: usize,
+) -> Result<(Vec<Vec<f64>>, Vec<Vec<f64>>), String> {
+    validate_hlc_columns(high, low, close)?;
+    let mut all_k = Vec::with_capacity(high.len());
+    let mut all_d = Vec::with_capacity(high.len());
+    for i in 0..high.len() {
+        let (k, d) = momentum::stoch(
+            &high[i],
+            &low[i],
+            &close[i],
+            fastk_period,
+            slowk_period,
+            slowd_period,
+        );
+        all_k.push(k);
+        all_d.push(d);
+    }
+    Ok((all_k, all_d))
+}
+
+// ---------------------------------------------------------------------------
+// batch_adx
+// ---------------------------------------------------------------------------
+
+/// Apply ADX to each set of (high, low, close) columns.
+pub fn batch_adx(
+    high: &[Vec<f64>],
+    low: &[Vec<f64>],
+    close: &[Vec<f64>],
+    timeperiod: usize,
+) -> Result<Vec<Vec<f64>>, String> {
+    if timeperiod == 0 {
+        return Err("timeperiod must be >= 1".into());
+    }
+    validate_hlc_columns(high, low, close)?;
+    Ok((0..high.len())
+        .map(|i| momentum::adx(&high[i], &low[i], &close[i], timeperiod))
+        .collect())
+}
+
+// ---------------------------------------------------------------------------
+// run_close_indicators
+// ---------------------------------------------------------------------------
+
+fn validate_indicator_requests(names: &[String], timeperiods: &[usize]) -> Result<(), String> {
+    if names.len() != timeperiods.len() {
+        return Err(format!(
+            "names length ({}) must equal timeperiods length ({})",
+            names.len(),
+            timeperiods.len()
+        ));
+    }
+    for (name, &tp) in names.iter().zip(timeperiods.iter()) {
+        if tp == 0 {
+            return Err(format!("{name}: timeperiod must be >= 1"));
+        }
+    }
+    Ok(())
+}
+
+fn compute_close_indicator(
+    name: &str,
+    close: &[f64],
+    timeperiod: usize,
+) -> Result<Vec<f64>, String> {
+    match name {
+        "SMA" => Ok(overlap::sma(close, timeperiod)),
+        "EMA" => Ok(overlap::ema(close, timeperiod)),
+        "RSI" => Ok(momentum::rsi(close, timeperiod)),
+        "STDDEV" => Ok(statistic::stddev(close, timeperiod, 1.0)),
+        "VAR" => Ok(statistic::stddev(close, timeperiod, 1.0)
+            .into_iter()
+            .map(|v| if v.is_nan() { v } else { v * v })
+            .collect()),
+        "LINEARREG" => {
+            let last_x = (timeperiod - 1) as f64;
+            Ok(rolling_linreg_apply(close, timeperiod, |slope, intercept| {
+                intercept + slope * last_x
+            }))
+        }
+        "LINEARREG_SLOPE" => Ok(rolling_linreg_apply(close, timeperiod, |slope, _| slope)),
+        "LINEARREG_INTERCEPT" => {
+            Ok(rolling_linreg_apply(close, timeperiod, |_, intercept| {
+                intercept
+            }))
+        }
+        "LINEARREG_ANGLE" => Ok(rolling_linreg_apply(close, timeperiod, |slope, _| {
+            slope.atan() * 180.0 / std::f64::consts::PI
+        })),
+        "TSF" => {
+            let forecast_x = timeperiod as f64;
+            Ok(rolling_linreg_apply(close, timeperiod, |slope, intercept| {
+                intercept + slope * forecast_x
+            }))
+        }
+        _ => Err(format!(
+            "unsupported close indicator for grouped execution: {name}"
+        )),
+    }
+}
+
+/// Run multiple close-only indicators on the same series.
+/// Returns `Vec<Result<Vec<f64>, String>>` — one result per (name, timeperiod) pair.
+pub fn run_close_indicators(
+    close: &[f64],
+    names: &[String],
+    timeperiods: &[usize],
+) -> Result<Vec<Vec<f64>>, String> {
+    validate_indicator_requests(names, timeperiods)?;
+    let mut results = Vec::with_capacity(names.len());
+    for (name, &tp) in names.iter().zip(timeperiods.iter()) {
+        results.push(compute_close_indicator(name, close, tp)?);
+    }
+    Ok(results)
+}
+
+// ---------------------------------------------------------------------------
+// run_hlc_indicators
+// ---------------------------------------------------------------------------
+
+fn compute_hlc_indicator(
+    name: &str,
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+) -> Result<Vec<f64>, String> {
+    match name {
+        "ATR" => Ok(volatility::atr(high, low, close, timeperiod)),
+        "NATR" => {
+            let atr_vals = volatility::atr(high, low, close, timeperiod);
+            Ok(atr_vals
+                .into_iter()
+                .zip(close.iter())
+                .map(|(a, &c)| {
+                    if a.is_nan() || c == 0.0 {
+                        f64::NAN
+                    } else {
+                        (a / c) * 100.0
+                    }
+                })
+                .collect())
+        }
+        "ADX" => Ok(momentum::adx(high, low, close, timeperiod)),
+        "ADXR" => Ok(momentum::adxr(high, low, close, timeperiod)),
+        "CCI" => Ok(compute_cci(high, low, close, timeperiod)),
+        "WILLR" => Ok(compute_willr(high, low, close, timeperiod)),
+        _ => Err(format!(
+            "unsupported HLC indicator for grouped execution: {name}"
+        )),
+    }
+}
+
+/// Run multiple HLC indicators on the same series.
+pub fn run_hlc_indicators(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    names: &[String],
+    timeperiods: &[usize],
+) -> Result<Vec<Vec<f64>>, String> {
+    validate_indicator_requests(names, timeperiods)?;
+    if high.len() != low.len() || high.len() != close.len() {
+        return Err("high, low, and close must have equal length".into());
+    }
+    let mut results = Vec::with_capacity(names.len());
+    for (name, &tp) in names.iter().zip(timeperiods.iter()) {
+        results.push(compute_hlc_indicator(name, high, low, close, tp)?);
+    }
+    Ok(results)
+}
+
+// ---------------------------------------------------------------------------
+// tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn close_data() -> Vec<f64> {
+        vec![
+            44.34, 44.09, 43.61, 44.33, 44.83, 45.10, 45.42, 45.84, 46.08, 45.89, 46.03, 45.61,
+            46.28, 46.28, 46.00, 46.03, 46.41, 46.22, 45.64,
+        ]
+    }
+
+    fn hlc_data() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
+        let close = close_data();
+        let high: Vec<f64> = close.iter().map(|c| c + 0.5).collect();
+        let low: Vec<f64> = close.iter().map(|c| c - 0.5).collect();
+        (high, low, close)
+    }
+
+    #[test]
+    fn test_batch_sma_basic() {
+        let col1 = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let col2 = vec![10.0, 20.0, 30.0, 40.0, 50.0];
+        let data = vec![col1, col2];
+        let result = batch_sma(&data, 3).unwrap();
+        assert_eq!(result.len(), 2);
+        assert!(result[0][0].is_nan());
+        assert!(result[0][1].is_nan());
+        assert!((result[0][2] - 2.0).abs() < 1e-10);
+        assert!((result[1][2] - 20.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_batch_sma_zero_period() {
+        let data = vec![vec![1.0, 2.0]];
+        assert!(batch_sma(&data, 0).is_err());
+    }
+
+    #[test]
+    fn test_batch_ema_basic() {
+        let data = vec![vec![1.0, 2.0, 3.0, 4.0, 5.0]];
+        let result = batch_ema(&data, 3).unwrap();
+        assert_eq!(result.len(), 1);
+        assert!(result[0][0].is_nan());
+    }
+
+    #[test]
+    fn test_batch_rsi_basic() {
+        let data = vec![close_data()];
+        let result = batch_rsi(&data, 14).unwrap();
+        assert_eq!(result.len(), 1);
+        // First 14 values should be NaN
+        for i in 0..14 {
+            assert!(result[0][i].is_nan(), "index {i} should be NaN");
+        }
+        // Value at index 14 should be a valid RSI
+        let rsi_val = result[0][14];
+        assert!(!rsi_val.is_nan());
+        assert!(rsi_val >= 0.0 && rsi_val <= 100.0);
+    }
+
+    #[test]
+    fn test_batch_atr_basic() {
+        let (h, l, c) = hlc_data();
+        let high = vec![h];
+        let low = vec![l];
+        let close = vec![c];
+        let result = batch_atr(&high, &low, &close, 14).unwrap();
+        assert_eq!(result.len(), 1);
+    }
+
+    #[test]
+    fn test_batch_stoch_basic() {
+        let (h, l, c) = hlc_data();
+        let high = vec![h];
+        let low = vec![l];
+        let close = vec![c];
+        let (k, d) = batch_stoch(&high, &low, &close, 5, 3, 3).unwrap();
+        assert_eq!(k.len(), 1);
+        assert_eq!(d.len(), 1);
+        assert_eq!(k[0].len(), d[0].len());
+    }
+
+    #[test]
+    fn test_batch_adx_basic() {
+        let (h, l, c) = hlc_data();
+        let high = vec![h];
+        let low = vec![l];
+        let close = vec![c];
+        let result = batch_adx(&high, &low, &close, 14).unwrap();
+        assert_eq!(result.len(), 1);
+    }
+
+    #[test]
+    fn test_run_close_indicators_basic() {
+        let close = close_data();
+        let names = vec!["SMA".to_string(), "EMA".to_string()];
+        let timeperiods = vec![5, 5];
+        let result = run_close_indicators(&close, &names, &timeperiods).unwrap();
+        assert_eq!(result.len(), 2);
+        assert_eq!(result[0].len(), close.len());
+        assert_eq!(result[1].len(), close.len());
+    }
+
+    #[test]
+    fn test_run_close_indicators_mismatched_lengths() {
+        let close = close_data();
+        let names = vec!["SMA".to_string()];
+        let timeperiods = vec![5, 10]; // different length
+        assert!(run_close_indicators(&close, &names, &timeperiods).is_err());
+    }
+
+    #[test]
+    fn test_run_close_indicators_linreg_variants() {
+        let close = close_data();
+        let names = vec![
+            "LINEARREG".to_string(),
+            "LINEARREG_SLOPE".to_string(),
+            "LINEARREG_INTERCEPT".to_string(),
+            "LINEARREG_ANGLE".to_string(),
+            "TSF".to_string(),
+        ];
+        let timeperiods = vec![5, 5, 5, 5, 5];
+        let result = run_close_indicators(&close, &names, &timeperiods).unwrap();
+        assert_eq!(result.len(), 5);
+        // First 4 values should be NaN for period=5
+        for series in &result {
+            for i in 0..4 {
+                assert!(series[i].is_nan());
+            }
+            assert!(!series[4].is_nan());
+        }
+    }
+
+    #[test]
+    fn test_run_hlc_indicators_basic() {
+        let (h, l, c) = hlc_data();
+        let names = vec!["ATR".to_string(), "CCI".to_string()];
+        let timeperiods = vec![14, 14];
+        let result = run_hlc_indicators(&h, &l, &c, &names, &timeperiods).unwrap();
+        assert_eq!(result.len(), 2);
+    }
+
+    #[test]
+    fn test_run_hlc_indicators_unsupported() {
+        let (h, l, c) = hlc_data();
+        let names = vec!["UNKNOWN".to_string()];
+        let timeperiods = vec![14];
+        assert!(run_hlc_indicators(&h, &l, &c, &names, &timeperiods).is_err());
+    }
+
+    #[test]
+    fn test_validate_hlc_mismatched_columns() {
+        let high = vec![vec![1.0, 2.0]];
+        let low = vec![vec![1.0, 2.0], vec![3.0, 4.0]]; // 2 cols vs 1
+        let close = vec![vec![1.0, 2.0]];
+        assert!(batch_atr(&high, &low, &close, 5).is_err());
+    }
+
+    #[test]
+    fn test_empty_data() {
+        let data: Vec<Vec<f64>> = vec![];
+        let result = batch_sma(&data, 3).unwrap();
+        assert!(result.is_empty());
+    }
+
+    #[test]
+    fn test_batch_multiple_columns() {
+        let data = vec![
+            vec![1.0, 2.0, 3.0, 4.0, 5.0],
+            vec![5.0, 4.0, 3.0, 2.0, 1.0],
+            vec![2.0, 4.0, 6.0, 8.0, 10.0],
+        ];
+        let result = batch_sma(&data, 3).unwrap();
+        assert_eq!(result.len(), 3);
+        // col 0: sma(3) at index 2 = (1+2+3)/3 = 2.0
+        assert!((result[0][2] - 2.0).abs() < 1e-10);
+        // col 1: sma(3) at index 2 = (5+4+3)/3 = 4.0
+        assert!((result[1][2] - 4.0).abs() < 1e-10);
+        // col 2: sma(3) at index 2 = (2+4+6)/3 = 4.0
+        assert!((result[2][2] - 4.0).abs() < 1e-10);
+    }
+}

crates/ferro_ta_core/src/chunked.rs

@@ -0,0 +1,123 @@
+//! Chunked / out-of-core execution helpers.
+//!
+//! - `trim_overlap`       — remove the first N elements from a slice
+//! - `stitch_chunks`      — concatenate trimmed chunk results
+//! - `make_chunk_ranges`  — compute (start, end) index pairs for chunked processing
+//! - `forward_fill_nan`   — forward-fill NaN values
+
+/// Remove the first `overlap` elements from a slice.
+pub fn trim_overlap(chunk_out: &[f64], overlap: usize) -> Vec<f64> {
+    if overlap > chunk_out.len() {
+        return vec![];
+    }
+    chunk_out[overlap..].to_vec()
+}
+
+/// Concatenate a list of slices into a single Vec.
+pub fn stitch_chunks(chunks: &[&[f64]]) -> Vec<f64> {
+    let mut out = Vec::new();
+    for &chunk in chunks {
+        out.extend_from_slice(chunk);
+    }
+    out
+}
+
+/// Compute (start, end) index pairs for chunked processing.
+///
+/// Returns a flat Vec of pairs: [start0, end0, start1, end1, ...].
+/// `chunk_size` is the desired output bars per chunk, `overlap` is the warm-up prefix.
+pub fn make_chunk_ranges(n: usize, chunk_size: usize, overlap: usize) -> Vec<i64> {
+    if chunk_size == 0 || n == 0 {
+        return vec![];
+    }
+    let mut ranges: Vec<i64> = Vec::new();
+    let mut start: usize = 0;
+    loop {
+        let end = (start + chunk_size + overlap).min(n);
+        ranges.push(start as i64);
+        ranges.push(end as i64);
+        if end >= n {
+            break;
+        }
+        start = end.saturating_sub(overlap);
+    }
+    ranges
+}
+
+/// Forward-fill NaN values in a 1-D array.
+/// Leading NaN values are preserved until the first non-NaN value appears.
+pub fn forward_fill_nan(values: &[f64]) -> Vec<f64> {
+    let mut out = Vec::with_capacity(values.len());
+    let mut last = f64::NAN;
+    for &value in values {
+        if value.is_nan() {
+            out.push(last);
+        } else {
+            last = value;
+            out.push(value);
+        }
+    }
+    out
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_trim_overlap() {
+        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let result = trim_overlap(&data, 2);
+        assert_eq!(result, vec![3.0, 4.0, 5.0]);
+    }
+
+    #[test]
+    fn test_trim_overlap_zero() {
+        let data = vec![1.0, 2.0, 3.0];
+        assert_eq!(trim_overlap(&data, 0), data);
+    }
+
+    #[test]
+    fn test_trim_overlap_exceeds() {
+        let data = vec![1.0, 2.0];
+        assert!(trim_overlap(&data, 5).is_empty());
+    }
+
+    #[test]
+    fn test_stitch_chunks() {
+        let a = vec![1.0, 2.0];
+        let b = vec![3.0, 4.0, 5.0];
+        let chunks: Vec<&[f64]> = vec![&a, &b];
+        let result = stitch_chunks(&chunks);
+        assert_eq!(result, vec![1.0, 2.0, 3.0, 4.0, 5.0]);
+    }
+
+    #[test]
+    fn test_make_chunk_ranges() {
+        let ranges = make_chunk_ranges(10, 4, 2);
+        // Expected: [0,6], [4,10]
+        assert_eq!(ranges.len() % 2, 0);
+        assert!(ranges.len() >= 4);
+        assert_eq!(ranges[0], 0);
+    }
+
+    #[test]
+    fn test_forward_fill_nan() {
+        let data = vec![f64::NAN, 1.0, f64::NAN, f64::NAN, 2.0, f64::NAN];
+        let result = forward_fill_nan(&data);
+        assert!(result[0].is_nan()); // leading NaN preserved
+        assert!((result[1] - 1.0).abs() < 1e-10);
+        assert!((result[2] - 1.0).abs() < 1e-10); // filled
+        assert!((result[3] - 1.0).abs() < 1e-10); // filled
+        assert!((result[4] - 2.0).abs() < 1e-10);
+        assert!((result[5] - 2.0).abs() < 1e-10); // filled
+    }
+
+    #[test]
+    fn test_empty() {
+        assert!(trim_overlap(&[], 0).is_empty());
+        assert!(stitch_chunks(&[]).is_empty());
+        assert!(make_chunk_ranges(0, 4, 2).is_empty());
+        assert!(forward_fill_nan(&[]).is_empty());
+    }
+}

crates/ferro_ta_core/src/commission.rs

@@ -0,0 +1,295 @@
+//! Commission, tax, and fee model for Indian and global markets.
+//!
+//! All `_rate` fields are fractions (0.001 = 0.1%).
+//! All per-unit fields (`flat_per_order`, `per_lot`) are in base currency units (e.g., INR).
+//! The model is self-contained: pass `trade_value`, `num_lots`, `is_buy` to get total cost.
+
+#[cfg(feature = "serde")]
+use serde::{Deserialize, Serialize};
+
+/// Advanced commission and tax model.
+///
+/// # Fields (all public for direct construction)
+/// - **Brokerage**: `flat_per_order`, `rate_of_value`, `per_lot`, `max_brokerage`
+/// - **STT**: `stt_rate`, `stt_on_buy`, `stt_on_sell`
+/// - **Levies**: `exchange_charges_rate`, `regulatory_charges_rate`, `gst_rate`, `stamp_duty_rate`
+/// - **Sizing**: `lot_size`
+///
+/// # Indian market notes
+/// - STT (Securities Transaction Tax) is applied on turnover (buy/sell legs vary by segment).
+/// - Exchange charges and regulatory body charges are on turnover.
+/// - GST (18%) applies on brokerage + exchange charges + regulatory body charges (not STT/stamp).
+/// - Stamp duty is on buy-side value only.
+#[derive(Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+pub struct CommissionModel {
+    // --- Brokerage ---------------------------------------------------------
+    /// Fixed fee per order (e.g., ₹20 flat fee per order). 0.0 = none.
+    pub flat_per_order: f64,
+    /// Proportional brokerage as fraction of `trade_value` (e.g., 0.001 = 0.1%). 0.0 = none.
+    pub rate_of_value: f64,
+    /// Fixed fee per lot (e.g., ₹2 per lot). 0.0 = none.
+    pub per_lot: f64,
+    /// Brokerage cap in currency units. 0.0 = no cap.
+    /// Effective brokerage = min(flat + rate × value + per_lot × lots, max_brokerage).
+    pub max_brokerage: f64,
+    /// Bid-ask spread model in basis points. Half-spread is paid on each leg (entry and exit),
+    /// so total roundtrip cost = spread_bps in bps. 0.0 = no spread cost.
+    pub spread_bps: f64,
+
+    // --- Securities Transaction Tax (STT) ----------------------------------
+    /// STT rate as fraction of trade value. 0.0 = no STT.
+    pub stt_rate: f64,
+    /// Apply STT on the buy leg.
+    pub stt_on_buy: bool,
+    /// Apply STT on the sell leg.
+    pub stt_on_sell: bool,
+
+    // --- Exchange & Regulatory Levies --------------------------------------
+    /// Exchange transaction charges rate (fraction of trade value).
+    pub exchange_charges_rate: f64,
+    /// Regulatory body turnover charges rate (fraction of trade value). Typically ~0.000001.
+    pub regulatory_charges_rate: f64,
+    /// Indirect tax (GST) rate applied on (brokerage + exchange_charges + regulatory_charges).
+    /// Typically 0.18 in India.
+    pub gst_rate: f64,
+    /// Stamp duty rate on buy side only (fraction of trade value).
+    pub stamp_duty_rate: f64,
+
+    // --- Instrument Sizing ------------------------------------------------
+    /// Lot size for the instrument.
+    /// Equities: 1.0. Index futures/options: contract lot size (e.g., 25, 50, 75).
+    /// Used for per_lot cost: cost += per_lot × ceil(quantity / lot_size).
+    pub lot_size: f64,
+
+    // --- Short Selling ----------------------------------------------------
+    /// Annualised short borrow rate as a fraction (e.g. 0.03 = 3% p.a.).
+    /// Applied per bar to short positions. 0.0 = no borrow cost.
+    pub short_borrow_rate_annual: f64,
+}
+
+impl Default for CommissionModel {
+    fn default() -> Self {
+        Self {
+            flat_per_order: 0.0,
+            rate_of_value: 0.0,
+            per_lot: 0.0,
+            max_brokerage: 0.0,
+            spread_bps: 0.0,
+            stt_rate: 0.0,
+            stt_on_buy: false,
+            stt_on_sell: false,
+            exchange_charges_rate: 0.0,
+            regulatory_charges_rate: 0.0,
+            gst_rate: 0.0,
+            stamp_duty_rate: 0.0,
+            lot_size: 1.0,
+            short_borrow_rate_annual: 0.0,
+        }
+    }
+}
+
+impl CommissionModel {
+    // ------------------------------------------------------------------
+    // Core computation
+    // ------------------------------------------------------------------
+
+    /// Compute total transaction cost in **absolute currency units**.
+    ///
+    /// # Parameters
+    /// - `trade_value`: price × quantity in base currency
+    /// - `num_lots`: number of lots transacted
+    /// - `is_buy`: true for buy (entry) leg, false for sell (exit) leg
+    pub fn total_cost(&self, trade_value: f64, num_lots: f64, is_buy: bool) -> f64 {
+        // Brokerage (optionally capped)
+        let raw_brokerage =
+            self.flat_per_order + self.rate_of_value * trade_value + self.per_lot * num_lots;
+        let brokerage = if self.max_brokerage > 0.0 {
+            raw_brokerage.min(self.max_brokerage)
+        } else {
+            raw_brokerage
+        };
+
+        // STT
+        let stt = if (is_buy && self.stt_on_buy) || (!is_buy && self.stt_on_sell) {
+            self.stt_rate * trade_value
+        } else {
+            0.0
+        };
+
+        let exchange = self.exchange_charges_rate * trade_value;
+        let regulatory = self.regulatory_charges_rate * trade_value;
+
+        // GST on brokerage + exchange + regulatory (NOT on STT or stamp duty)
+        let gst = self.gst_rate * (brokerage + exchange + regulatory);
+
+        // Stamp duty only on buy side
+        let stamp = if is_buy {
+            self.stamp_duty_rate * trade_value
+        } else {
+            0.0
+        };
+
+        // Bid-ask spread: half-spread paid on each leg
+        let spread_cost = self.spread_bps / 2.0 / 10_000.0 * trade_value;
+
+        brokerage + stt + exchange + regulatory + gst + stamp + spread_cost
+    }
+
+    /// Borrow cost per bar for a short position.
+    ///
+    /// # Parameters
+    /// - `trade_value`: abs(price × quantity)
+    /// - `periods_per_year`: 252 for daily, 52 for weekly, etc.
+    pub fn short_borrow_cost(&self, trade_value: f64, periods_per_year: f64) -> f64 {
+        if self.short_borrow_rate_annual <= 0.0 || periods_per_year <= 0.0 {
+            return 0.0;
+        }
+        self.short_borrow_rate_annual / periods_per_year * trade_value
+    }
+
+    /// Compute cost as a **fraction of `initial_capital`** for use in normalised equity loops.
+    ///
+    /// Returns 0.0 if `initial_capital` ≤ 0.
+    pub fn cost_fraction(
+        &self,
+        trade_value: f64,
+        num_lots: f64,
+        is_buy: bool,
+        initial_capital: f64,
+    ) -> f64 {
+        if initial_capital <= 0.0 {
+            return 0.0;
+        }
+        self.total_cost(trade_value, num_lots, is_buy) / initial_capital
+    }
+
+    // ------------------------------------------------------------------
+    // Built-in Presets
+    // ------------------------------------------------------------------
+
+    /// Zero commission — useful for clean research/comparison runs.
+    pub fn zero() -> Self {
+        Self::default()
+    }
+
+    /// Indian equity **delivery** (long-term hold).
+    ///
+    /// Brokerage: 0.1% (capped at ₹20), STT 0.1% both sides,
+    /// exchange charges, regulatory body charges, 18% GST, stamp duty.
+    pub fn equity_delivery_india() -> Self {
+        Self {
+            flat_per_order: 0.0,
+            rate_of_value: 0.001, // 0.1%
+            per_lot: 0.0,
+            max_brokerage: 20.0, // ₹20 cap
+            spread_bps: 0.0,
+            stt_rate: 0.001, // 0.1%
+            stt_on_buy: true,
+            stt_on_sell: true,
+            exchange_charges_rate: 0.0000297,
+            regulatory_charges_rate: 0.000001,
+            gst_rate: 0.18,
+            stamp_duty_rate: 0.00015,
+            lot_size: 1.0,
+            short_borrow_rate_annual: 0.0,
+        }
+    }
+
+    /// Indian equity **intraday** (same-day square-off).
+    ///
+    /// Brokerage: 0.03% (capped at ₹20), STT 0.025% sell side only,
+    /// exchange charges, regulatory body charges, 18% GST, stamp duty on buy.
+    pub fn equity_intraday_india() -> Self {
+        Self {
+            flat_per_order: 0.0,
+            rate_of_value: 0.0003, // 0.03%
+            per_lot: 0.0,
+            max_brokerage: 20.0,
+            spread_bps: 0.0,
+            stt_rate: 0.00025, // 0.025%
+            stt_on_buy: false,
+            stt_on_sell: true,
+            exchange_charges_rate: 0.0000297,
+            regulatory_charges_rate: 0.000001,
+            gst_rate: 0.18,
+            stamp_duty_rate: 0.000003,
+            lot_size: 1.0,
+            short_borrow_rate_annual: 0.0,
+        }
+    }
+
+    /// Indian **index futures** (indicative rates per current regulations).
+    ///
+    /// Flat ₹20 per order, STT 0.05% sell side only, exchange charges,
+    /// regulatory body charges, 18% GST, stamp duty on buy.
+    /// `lot_size` defaults to 25 — update as needed for the specific contract.
+    pub fn futures_india() -> Self {
+        Self {
+            flat_per_order: 20.0,
+            rate_of_value: 0.0,
+            per_lot: 0.0,
+            max_brokerage: 0.0,
+            spread_bps: 0.0,
+            stt_rate: 0.0005, // 0.05%
+            stt_on_buy: false,
+            stt_on_sell: true,
+            exchange_charges_rate: 0.0000019,
+            regulatory_charges_rate: 0.000001,
+            gst_rate: 0.18,
+            stamp_duty_rate: 0.00002,
+            lot_size: 25.0,
+            short_borrow_rate_annual: 0.0,
+        }
+    }
+
+    /// Indian **index options** (indicative rates per current regulations).
+    ///
+    /// Flat ₹20 per order, STT 0.15% on premium sell side only, exchange charges,
+    /// regulatory body charges, 18% GST, stamp duty on buy.
+    /// `lot_size` defaults to 25 — update as needed for the specific contract.
+    pub fn options_india() -> Self {
+        Self {
+            flat_per_order: 20.0,
+            rate_of_value: 0.0,
+            per_lot: 0.0,
+            max_brokerage: 0.0,
+            spread_bps: 0.0,
+            stt_rate: 0.0015, // 0.15% on premium
+            stt_on_buy: false,
+            stt_on_sell: true,
+            exchange_charges_rate: 0.0000053,
+            regulatory_charges_rate: 0.000001,
+            gst_rate: 0.18,
+            stamp_duty_rate: 0.000003,
+            lot_size: 25.0,
+            short_borrow_rate_annual: 0.0,
+        }
+    }
+
+    /// Simple proportional model — e.g., `proportional(0.001)` = 0.1% both sides.
+    ///
+    /// No taxes, no levies — suitable for non-Indian markets or simplified modelling.
+    pub fn proportional(rate: f64) -> Self {
+        Self {
+            rate_of_value: rate,
+            ..Default::default()
+        }
+    }
+
+    // ------------------------------------------------------------------
+    // JSON serialization (requires "serde" feature)
+    // ------------------------------------------------------------------
+
+    /// Serialize to a pretty-printed JSON string.
+    #[cfg(feature = "serde")]
+    pub fn to_json(&self) -> Result<String, serde_json::Error> {
+        serde_json::to_string_pretty(self)
+    }
+
+    /// Deserialize from a JSON string.
+    #[cfg(feature = "serde")]
+    pub fn from_json(s: &str) -> Result<Self, serde_json::Error> {
+        serde_json::from_str(s)
+    }
+}

crates/ferro_ta_core/src/crypto.rs

@@ -0,0 +1,91 @@
+//! Crypto and 24/7 market helpers.
+//!
+//! - `funding_cumulative_pnl` — cumulative PnL from periodic funding rate payments
+//! - `continuous_bar_labels`  — assign sequential integer labels based on fixed period size
+//! - `mark_session_boundaries` — return indices where a new UTC day begins
+
+/// Compute the cumulative PnL from funding rate payments.
+///
+/// `position_size` and `funding_rate` must have the same length.
+/// PnL at period i = -position_size[i] * funding_rate[i] (longs pay when rate > 0).
+pub fn funding_cumulative_pnl(position_size: &[f64], funding_rate: &[f64]) -> Vec<f64> {
+    let n = position_size.len();
+    let mut out = vec![0.0_f64; n];
+    let mut cumulative = 0.0_f64;
+    for i in 0..n {
+        cumulative += -position_size[i] * funding_rate[i];
+        out[i] = cumulative;
+    }
+    out
+}
+
+/// Assign a sequential integer label per bar based on a fixed-size period.
+///
+/// Bars 0..(period_bars-1) get label 0, bars period_bars..(2*period_bars-1) get label 1, etc.
+/// `period_bars` must be >= 1.
+pub fn continuous_bar_labels(n_bars: usize, period_bars: usize) -> Vec<i64> {
+    (0..n_bars).map(|i| (i / period_bars) as i64).collect()
+}
+
+/// Return bar indices where a new UTC day begins (based on nanosecond timestamps).
+///
+/// Bar 0 is always included as the first boundary.
+pub fn mark_session_boundaries(timestamps_ns: &[i64]) -> Vec<i64> {
+    let n = timestamps_ns.len();
+    if n == 0 {
+        return vec![];
+    }
+    const NS_PER_DAY: i64 = 86_400_000_000_000;
+    let mut out = vec![0i64]; // bar 0 is always a boundary
+    let mut prev_day = timestamps_ns[0].div_euclid(NS_PER_DAY);
+    for (i, &t) in timestamps_ns.iter().enumerate().skip(1) {
+        let day = t.div_euclid(NS_PER_DAY);
+        if day != prev_day {
+            out.push(i as i64);
+            prev_day = day;
+        }
+    }
+    out
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_funding_cumulative_pnl() {
+        let pos = vec![100.0, 100.0, -50.0];
+        let rate = vec![0.001, -0.002, 0.001];
+        let result = funding_cumulative_pnl(&pos, &rate);
+        assert!((result[0] - (-0.1)).abs() < 1e-10);
+        assert!((result[1] - 0.1).abs() < 1e-10); // -0.1 + 0.2 = 0.1
+        assert!((result[2] - 0.15).abs() < 1e-10); // 0.1 + 0.05 = 0.15
+    }
+
+    #[test]
+    fn test_continuous_bar_labels() {
+        let labels = continuous_bar_labels(7, 3);
+        assert_eq!(labels, vec![0, 0, 0, 1, 1, 1, 2]);
+    }
+
+    #[test]
+    fn test_mark_session_boundaries() {
+        let ns_per_day: i64 = 86_400_000_000_000;
+        let ts = vec![
+            0,                         // day 0
+            ns_per_day / 2,            // day 0
+            ns_per_day,                // day 1
+            ns_per_day + ns_per_day / 2, // day 1
+            ns_per_day * 2,            // day 2
+        ];
+        let result = mark_session_boundaries(&ts);
+        assert_eq!(result, vec![0, 2, 4]);
+    }
+
+    #[test]
+    fn test_empty() {
+        assert!(funding_cumulative_pnl(&[], &[]).is_empty());
+        assert!(continuous_bar_labels(0, 1).is_empty());
+        assert!(mark_session_boundaries(&[]).is_empty());
+    }
+}

crates/ferro_ta_core/src/currency.rs

@@ -0,0 +1,173 @@
+//! Currency metadata and Indian number formatting.
+
+/// Immutable currency descriptor.
+///
+/// Carries the currency code, symbol, decimal places, and whether to use
+/// Indian lakh/crore grouping (1,23,45,678.00) instead of standard
+/// Western grouping (1,234,567.89).
+#[derive(Clone, Debug, PartialEq, Eq)]
+pub struct Currency {
+    /// IETF currency code, e.g. "INR", "USD".
+    pub code: &'static str,
+    /// Display symbol, e.g. "₹", "$".
+    pub symbol: &'static str,
+    /// Number of decimal places for formatting.
+    pub decimal_places: u8,
+    /// Use Indian lakh/crore digit grouping (true only for INR).
+    pub lakh_grouping: bool,
+}
+
+impl Currency {
+    pub const INR: Currency = Currency {
+        code: "INR",
+        symbol: "₹",
+        decimal_places: 2,
+        lakh_grouping: true,
+    };
+    pub const USD: Currency = Currency {
+        code: "USD",
+        symbol: "$",
+        decimal_places: 2,
+        lakh_grouping: false,
+    };
+    pub const EUR: Currency = Currency {
+        code: "EUR",
+        symbol: "€",
+        decimal_places: 2,
+        lakh_grouping: false,
+    };
+    pub const GBP: Currency = Currency {
+        code: "GBP",
+        symbol: "£",
+        decimal_places: 2,
+        lakh_grouping: false,
+    };
+    pub const JPY: Currency = Currency {
+        code: "JPY",
+        symbol: "¥",
+        decimal_places: 0,
+        lakh_grouping: false,
+    };
+    pub const USDT: Currency = Currency {
+        code: "USDT",
+        symbol: "₮",
+        decimal_places: 2,
+        lakh_grouping: false,
+    };
+
+    /// Look up a currency by IETF code (case-insensitive).
+    /// Returns `None` if the code is not recognised.
+    pub fn from_code(code: &str) -> Option<&'static Currency> {
+        match code.to_ascii_uppercase().as_str() {
+            "INR" => Some(&Currency::INR),
+            "USD" => Some(&Currency::USD),
+            "EUR" => Some(&Currency::EUR),
+            "GBP" => Some(&Currency::GBP),
+            "JPY" => Some(&Currency::JPY),
+            "USDT" => Some(&Currency::USDT),
+            _ => None,
+        }
+    }
+
+    /// Format `amount` according to this currency's style.
+    ///
+    /// - INR uses Indian lakh/crore grouping: `₹1,23,45,678.00`
+    /// - Others use standard Western grouping: `$1,234,567.89`
+    pub fn format(&self, amount: f64) -> String {
+        let neg = amount < 0.0;
+        let abs = amount.abs();
+        let integer_part = abs.floor() as u64;
+        let frac_part = abs - abs.floor();
+
+        let grouped = if self.lakh_grouping {
+            format_lakh(integer_part)
+        } else {
+            format_standard(integer_part)
+        };
+
+        let dp = self.decimal_places as usize;
+        let decimal_str = if dp > 0 {
+            let frac = (frac_part * 10f64.powi(dp as i32)).round() as u64;
+            format!(".{:0>width$}", frac, width = dp)
+        } else {
+            String::new()
+        };
+
+        let sign = if neg { "-" } else { "" };
+        format!("{}{}{}{}", sign, self.symbol, grouped, decimal_str)
+    }
+}
+
+/// Indian lakh/crore grouping: last 3 digits, then groups of 2 from the right.
+/// e.g. 12345678 → "1,23,45,678"
+fn format_lakh(n: u64) -> String {
+    let s = n.to_string();
+    if s.len() <= 3 {
+        return s;
+    }
+    let (rest, last3) = s.split_at(s.len() - 3);
+    let mut out = String::new();
+    let chars: Vec<char> = rest.chars().collect();
+    let first_len = chars.len() % 2;
+    if first_len > 0 {
+        out.push_str(&chars[..first_len].iter().collect::<String>());
+    }
+    let mut i = first_len;
+    while i < chars.len() {
+        if !out.is_empty() {
+            out.push(',');
+        }
+        out.push_str(&chars[i..i + 2].iter().collect::<String>());
+        i += 2;
+    }
+    if !out.is_empty() {
+        out.push(',');
+    }
+    out.push_str(last3);
+    out
+}
+
+/// Standard Western grouping: groups of 3 digits from the right.
+/// e.g. 1234567 → "1,234,567"
+fn format_standard(n: u64) -> String {
+    let s = n.to_string();
+    let mut out = String::new();
+    for (i, c) in s.chars().rev().enumerate() {
+        if i > 0 && i % 3 == 0 {
+            out.push(',');
+        }
+        out.push(c);
+    }
+    out.chars().rev().collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_inr_format() {
+        assert_eq!(Currency::INR.format(123456.78), "₹1,23,456.78");
+        assert_eq!(Currency::INR.format(10000000.0), "₹1,00,00,000.00");
+        assert_eq!(Currency::INR.format(100.0), "₹100.00");
+        assert_eq!(Currency::INR.format(-5000.0), "-₹5,000.00");
+    }
+
+    #[test]
+    fn test_usd_format() {
+        assert_eq!(Currency::USD.format(1234567.89), "$1,234,567.89");
+        assert_eq!(Currency::USD.format(0.5), "$0.50");
+    }
+
+    #[test]
+    fn test_jpy_format() {
+        assert_eq!(Currency::JPY.format(1000000.0), "¥1,000,000");
+    }
+
+    #[test]
+    fn test_from_code() {
+        assert_eq!(Currency::from_code("inr"), Some(&Currency::INR));
+        assert_eq!(Currency::from_code("USD"), Some(&Currency::USD));
+        assert_eq!(Currency::from_code("UNKNOWN"), None);
+    }
+}

crates/ferro_ta_core/src/cycle.rs

@@ -0,0 +1,370 @@
+//! Cycle indicators — Hilbert Transform-based cycle analysis (Ehlers).
+//!
+//! Based on John Ehlers' Discrete Hilbert Transform as implemented in TA-Lib.
+//! Reference: "Cybernetic Analysis for Stocks and Futures" by J.F. Ehlers
+//!
+//! All HT functions share a 63-bar lookback period.
+
+use std::f64::consts::PI;
+
+/// Number of leading bars that are set to NaN / zero.
+pub const HT_LOOKBACK: usize = 63;
+
+/// Shared output from the core Hilbert Transform computation.
+pub struct HtCore {
+    pub trendline: Vec<f64>,
+    pub dc_period: Vec<f64>,
+    pub dc_phase: Vec<f64>,
+    pub inphase: Vec<f64>,
+    pub quadrature: Vec<f64>,
+    pub trend_mode: Vec<i32>,
+}
+
+/// Run the full Hilbert Transform pipeline on a slice of close prices.
+pub fn compute_ht_core(prices: &[f64]) -> HtCore {
+    let n = prices.len();
+
+    let mut trendline = vec![f64::NAN; n];
+    let mut dc_period = vec![f64::NAN; n];
+    let mut dc_phase = vec![f64::NAN; n];
+    let mut inphase = vec![f64::NAN; n];
+    let mut quadrature = vec![f64::NAN; n];
+    let mut trend_mode = vec![0i32; n];
+
+    if n <= HT_LOOKBACK {
+        return HtCore {
+            trendline,
+            dc_period,
+            dc_phase,
+            inphase,
+            quadrature,
+            trend_mode,
+        };
+    }
+
+    // Step 1: Smooth the price series (4-bar weighted average)
+    let mut smooth = vec![0.0f64; n];
+    for i in 0..n {
+        smooth[i] = if i >= 3 {
+            (4.0 * prices[i] + 3.0 * prices[i - 1] + 2.0 * prices[i - 2] + prices[i - 3]) / 10.0
+        } else {
+            prices[i]
+        };
+    }
+
+    // Step 2: Full Hilbert Transform pipeline
+    let mut detrender = vec![0.0f64; n];
+    let mut q1 = vec![0.0f64; n];
+    let mut i1 = vec![0.0f64; n];
+    let mut ji = vec![0.0f64; n];
+    let mut jq = vec![0.0f64; n];
+    let mut i2 = vec![0.0f64; n];
+    let mut q2 = vec![0.0f64; n];
+    let mut re = vec![0.0f64; n];
+    let mut im = vec![0.0f64; n];
+    let mut period = vec![0.0f64; n];
+    let mut smooth_period = vec![0.0f64; n];
+    let mut phase = vec![0.0f64; n];
+
+    for i in 6..n {
+        let prev_period = period[i - 1];
+        // Alpha coefficient for HT filters depends on the current period estimate
+        let alpha = 0.075 * prev_period + 0.54;
+
+        // Discrete Hilbert Transform of smooth price (detrender)
+        detrender[i] = (0.0962 * smooth[i] + 0.5769 * smooth[i - 2]
+            - 0.5769 * smooth[i - 4]
+            - 0.0962 * smooth[i - 6])
+            * alpha;
+
+        // Q1: HT of detrender
+        if i >= 12 {
+            q1[i] = (0.0962 * detrender[i] + 0.5769 * detrender[i - 2]
+                - 0.5769 * detrender[i - 4]
+                - 0.0962 * detrender[i - 6])
+                * alpha;
+        }
+
+        // I1: delayed detrender
+        if i >= 9 {
+            i1[i] = detrender[i - 3];
+        }
+
+        // jI: HT of I1
+        if i >= 15 {
+            ji[i] = (0.0962 * i1[i] + 0.5769 * i1[i - 2] - 0.5769 * i1[i - 4] - 0.0962 * i1[i - 6])
+                * alpha;
+        }
+
+        // jQ: HT of Q1
+        if i >= 18 {
+            jq[i] = (0.0962 * q1[i] + 0.5769 * q1[i - 2] - 0.5769 * q1[i - 4] - 0.0962 * q1[i - 6])
+                * alpha;
+        }
+
+        // Phase components
+        let i2_raw = i1[i] - jq[i];
+        let q2_raw = q1[i] + ji[i];
+
+        // EMA smoothing of I2 and Q2
+        let i2_prev = i2[i - 1];
+        let q2_prev = q2[i - 1];
+        i2[i] = 0.2 * i2_raw + 0.8 * i2_prev;
+        q2[i] = 0.2 * q2_raw + 0.8 * q2_prev;
+
+        // Cross-product for period estimation
+        let re_raw = i2[i] * i2_prev + q2[i] * q2_prev;
+        let im_raw = i2[i] * q2_prev - q2[i] * i2_prev;
+
+        // EMA smoothing of Re and Im
+        re[i] = 0.2 * re_raw + 0.8 * re[i - 1];
+        im[i] = 0.2 * im_raw + 0.8 * im[i - 1];
+
+        // Compute period from cross-product of consecutive phasors.
+        let mut p = if re[i] != 0.0 && im[i] != 0.0 && re[i] > 0.0 {
+            2.0 * PI / (im[i] / re[i]).atan()
+        } else {
+            prev_period
+        };
+
+        // Clamp period relative to previous
+        if prev_period > 0.0 {
+            if p > 1.5 * prev_period {
+                p = 1.5 * prev_period;
+            }
+            if p < 0.67 * prev_period {
+                p = 0.67 * prev_period;
+            }
+        }
+        // Hard clamp to [6, 50] bars
+        p = p.clamp(6.0, 50.0);
+
+        // EMA smooth the period
+        period[i] = 0.2 * p + 0.8 * prev_period;
+
+        // Smooth the smoothed period once more
+        smooth_period[i] = 0.33 * period[i] + 0.67 * smooth_period[i - 1];
+
+        // Phase from I1 and Q1
+        phase[i] = if i1[i] != 0.0 {
+            q1[i].atan2(i1[i]) * 180.0 / PI
+        } else if q1[i] > 0.0 {
+            90.0
+        } else if q1[i] < 0.0 {
+            -90.0
+        } else {
+            0.0
+        };
+
+        // Write outputs once past lookback
+        if i >= HT_LOOKBACK {
+            dc_period[i] = smooth_period[i];
+            dc_phase[i] = phase[i];
+            inphase[i] = i1[i];
+            quadrature[i] = q1[i];
+
+            // Trend mode: cycle when SmoothPeriod >= 20, trend when < 20
+            trend_mode[i] = if smooth_period[i] < 20.0 { 1 } else { 0 };
+        }
+    }
+
+    // Trendline: average over the current dominant cycle period
+    for i in HT_LOOKBACK..n {
+        let sp = smooth_period[i];
+        let dc = (sp.round() as usize).max(1).min(i + 1);
+        let sum: f64 = (0..dc).map(|j| smooth[i - j]).sum();
+        trendline[i] = sum / dc as f64;
+    }
+
+    HtCore {
+        trendline,
+        dc_period,
+        dc_phase,
+        inphase,
+        quadrature,
+        trend_mode,
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Public indicator functions
+// ---------------------------------------------------------------------------
+
+/// Hilbert Transform Instantaneous Trendline (Ehlers).
+/// Smooths price over the dominant cycle period.
+pub fn ht_trendline(close: &[f64]) -> Vec<f64> {
+    compute_ht_core(close).trendline
+}
+
+/// Hilbert Transform Dominant Cycle Period in bars.
+pub fn ht_dcperiod(close: &[f64]) -> Vec<f64> {
+    compute_ht_core(close).dc_period
+}
+
+/// Hilbert Transform Dominant Cycle Phase in degrees.
+pub fn ht_dcphase(close: &[f64]) -> Vec<f64> {
+    compute_ht_core(close).dc_phase
+}
+
+/// Hilbert Transform Phasor components. Returns `(inphase, quadrature)`.
+pub fn ht_phasor(close: &[f64]) -> (Vec<f64>, Vec<f64>) {
+    let core = compute_ht_core(close);
+    (core.inphase, core.quadrature)
+}
+
+/// Hilbert Transform SineWave. Returns `(sine, leadsine)` where leadsine
+/// leads sine by 45 degrees.
+pub fn ht_sine(close: &[f64]) -> (Vec<f64>, Vec<f64>) {
+    let n = close.len();
+    let core = compute_ht_core(close);
+
+    let mut sine = vec![f64::NAN; n];
+    let mut lead_sine = vec![f64::NAN; n];
+
+    for i in HT_LOOKBACK..n {
+        if !core.dc_phase[i].is_nan() {
+            let phase_rad = core.dc_phase[i] * PI / 180.0;
+            sine[i] = phase_rad.sin();
+            lead_sine[i] = (phase_rad + PI / 4.0).sin(); // 45-degree lead
+        }
+    }
+
+    (sine, lead_sine)
+}
+
+/// Hilbert Transform Trend vs Cycle Mode: 1 = trending, 0 = cycling.
+pub fn ht_trendmode(close: &[f64]) -> Vec<i32> {
+    compute_ht_core(close).trend_mode
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Generate a simple sine wave for testing cycle detection.
+    fn sine_wave(n: usize, period: f64) -> Vec<f64> {
+        (0..n)
+            .map(|i| 100.0 + 10.0 * (2.0 * PI * i as f64 / period).sin())
+            .collect()
+    }
+
+    /// Flat price series for baseline testing.
+    fn flat_prices(n: usize) -> Vec<f64> {
+        vec![100.0; n]
+    }
+
+    #[test]
+    fn test_ht_trendline_length_and_lookback() {
+        let close = sine_wave(200, 20.0);
+        let result = ht_trendline(&close);
+        assert_eq!(result.len(), close.len());
+        // First HT_LOOKBACK values must be NaN
+        for v in &result[..HT_LOOKBACK] {
+            assert!(v.is_nan(), "expected NaN in lookback region");
+        }
+        // Values after lookback must be finite
+        for v in &result[HT_LOOKBACK..] {
+            assert!(v.is_finite(), "expected finite value after lookback");
+        }
+    }
+
+    #[test]
+    fn test_ht_dcperiod_length_and_lookback() {
+        let close = sine_wave(200, 20.0);
+        let result = ht_dcperiod(&close);
+        assert_eq!(result.len(), close.len());
+        for v in &result[..HT_LOOKBACK] {
+            assert!(v.is_nan());
+        }
+        // After lookback, period should be positive and finite
+        for v in &result[HT_LOOKBACK..] {
+            assert!(v.is_finite());
+            assert!(*v >= 6.0 && *v <= 50.0, "period {} out of [6,50]", v);
+        }
+    }
+
+    #[test]
+    fn test_ht_dcphase_length_and_lookback() {
+        let close = sine_wave(200, 20.0);
+        let result = ht_dcphase(&close);
+        assert_eq!(result.len(), close.len());
+        for v in &result[..HT_LOOKBACK] {
+            assert!(v.is_nan());
+        }
+        for v in &result[HT_LOOKBACK..] {
+            assert!(v.is_finite());
+        }
+    }
+
+    #[test]
+    fn test_ht_phasor_dual_output() {
+        let close = sine_wave(200, 20.0);
+        let (inp, quad) = ht_phasor(&close);
+        assert_eq!(inp.len(), close.len());
+        assert_eq!(quad.len(), close.len());
+        for v in &inp[..HT_LOOKBACK] {
+            assert!(v.is_nan());
+        }
+        for v in &quad[..HT_LOOKBACK] {
+            assert!(v.is_nan());
+        }
+    }
+
+    #[test]
+    fn test_ht_sine_dual_output() {
+        let close = sine_wave(200, 20.0);
+        let (s, ls) = ht_sine(&close);
+        assert_eq!(s.len(), close.len());
+        assert_eq!(ls.len(), close.len());
+        for v in &s[..HT_LOOKBACK] {
+            assert!(v.is_nan());
+        }
+        // Sine values should be in [-1, 1]
+        for v in &s[HT_LOOKBACK..] {
+            assert!(v.is_finite());
+            assert!(*v >= -1.0 && *v <= 1.0, "sine {} out of [-1,1]", v);
+        }
+        for v in &ls[HT_LOOKBACK..] {
+            assert!(v.is_finite());
+            assert!(*v >= -1.0 && *v <= 1.0, "leadsine {} out of [-1,1]", v);
+        }
+    }
+
+    #[test]
+    fn test_ht_trendmode_values() {
+        let close = sine_wave(200, 20.0);
+        let result = ht_trendmode(&close);
+        assert_eq!(result.len(), close.len());
+        // All values must be 0 or 1
+        for v in &result {
+            assert!(*v == 0 || *v == 1, "trend_mode {} not 0 or 1", v);
+        }
+    }
+
+    #[test]
+    fn test_short_input_all_nan() {
+        let close = vec![100.0; HT_LOOKBACK]; // exactly HT_LOOKBACK, not enough
+        let tl = ht_trendline(&close);
+        assert!(tl.iter().all(|v| v.is_nan()));
+        let dp = ht_dcperiod(&close);
+        assert!(dp.iter().all(|v| v.is_nan()));
+    }
+
+    #[test]
+    fn test_flat_prices_trendline_equals_price() {
+        let close = flat_prices(200);
+        let tl = ht_trendline(&close);
+        // For a flat price, trendline after lookback should be very close to the price
+        for v in &tl[HT_LOOKBACK..] {
+            assert!(
+                (v - 100.0).abs() < 1e-6,
+                "trendline {} diverged from flat price",
+                v
+            );
+        }
+    }
+}

crates/ferro_ta_core/src/extended.rs

@@ -0,0 +1,977 @@
+//! Extended indicators — pure Rust implementations (no PyO3, no numpy).
+//!
+//! These indicators are not part of TA-Lib and provide additional technical
+//! analysis capabilities. All functions operate on `&[f64]` slices and return
+//! `Vec<f64>` (or tuples thereof).
+
+#![allow(clippy::too_many_arguments)]
+
+use crate::math;
+use crate::overlap;
+// Note: we use a local compute_atr helper (seeds from bar 0) rather than
+// crate::volatility::atr (which seeds from bar 1, TA-Lib style).
+
+// ---------------------------------------------------------------------------
+// Internal helpers
+// ---------------------------------------------------------------------------
+
+/// Compute ATR array using Wilder smoothing (same algorithm as in the PyO3
+/// extended module — seeds from bar 0, not bar 1 like TA-Lib's `volatility::atr`).
+fn compute_atr(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = high.len();
+    let mut result = vec![f64::NAN; n];
+    if n <= timeperiod {
+        return result;
+    }
+    // Seed: SMA of first `timeperiod` true range values
+    let mut seed_sum = high[0] - low[0]; // first TR has no prev_close
+    for i in 1..timeperiod {
+        let hl = high[i] - low[i];
+        let hc = (high[i] - close[i - 1]).abs();
+        let lc = (low[i] - close[i - 1]).abs();
+        seed_sum += hl.max(hc).max(lc);
+    }
+    let mut atr = seed_sum / timeperiod as f64;
+    result[timeperiod - 1] = atr;
+    let pf = (timeperiod - 1) as f64;
+    for i in timeperiod..n {
+        let hl = high[i] - low[i];
+        let hc = (high[i] - close[i - 1]).abs();
+        let lc = (low[i] - close[i - 1]).abs();
+        let tr = hl.max(hc).max(lc);
+        atr = (atr * pf + tr) / timeperiod as f64;
+        result[i] = atr;
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// VWAP
+// ---------------------------------------------------------------------------
+
+/// Volume Weighted Average Price (cumulative or rolling).
+///
+/// # Arguments
+/// * `high`, `low`, `close`, `volume` — equal-length price/volume slices.
+/// * `timeperiod` — 0 for cumulative VWAP from bar 0; >= 1 for a rolling window.
+///
+/// # Returns
+/// A `Vec<f64>` of VWAP values. For rolling mode the first `timeperiod - 1`
+/// entries are `NaN`.
+pub fn vwap(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    volume: &[f64],
+    timeperiod: usize,
+) -> Vec<f64> {
+    let n = high.len();
+    let mut result = vec![f64::NAN; n];
+
+    if timeperiod == 0 {
+        let mut cum_tpv = 0.0_f64;
+        let mut cum_vol = 0.0_f64;
+        for i in 0..n {
+            let tp = (high[i] + low[i] + close[i]) / 3.0;
+            cum_tpv += tp * volume[i];
+            cum_vol += volume[i];
+            result[i] = if cum_vol != 0.0 {
+                cum_tpv / cum_vol
+            } else {
+                f64::NAN
+            };
+        }
+    } else {
+        // Pre-compute cumulative sums for O(n) rolling window
+        let mut cum_tpv_arr = vec![0.0_f64; n];
+        let mut cum_vol_arr = vec![0.0_f64; n];
+        for i in 0..n {
+            let tp = (high[i] + low[i] + close[i]) / 3.0;
+            let tpv = tp * volume[i];
+            cum_tpv_arr[i] = tpv + if i > 0 { cum_tpv_arr[i - 1] } else { 0.0 };
+            cum_vol_arr[i] = volume[i] + if i > 0 { cum_vol_arr[i - 1] } else { 0.0 };
+        }
+        for i in (timeperiod - 1)..n {
+            let prev_tpv = if i >= timeperiod {
+                cum_tpv_arr[i - timeperiod]
+            } else {
+                0.0
+            };
+            let prev_vol = if i >= timeperiod {
+                cum_vol_arr[i - timeperiod]
+            } else {
+                0.0
+            };
+            let w_tpv = cum_tpv_arr[i] - prev_tpv;
+            let w_vol = cum_vol_arr[i] - prev_vol;
+            result[i] = if w_vol != 0.0 {
+                w_tpv / w_vol
+            } else {
+                f64::NAN
+            };
+        }
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// VWMA
+// ---------------------------------------------------------------------------
+
+/// Volume Weighted Moving Average.
+///
+/// `VWMA = sum(close * volume, n) / sum(volume, n)`
+///
+/// # Arguments
+/// * `close` — price series.
+/// * `volume` — volume series (same length as `close`).
+/// * `timeperiod` — rolling window size (>= 1).
+///
+/// # Returns
+/// A `Vec<f64>` with `NaN` for the first `timeperiod - 1` entries.
+pub fn vwma(close: &[f64], volume: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = close.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod < 1 || n < timeperiod {
+        return result;
+    }
+
+    let mut cum_cv = vec![0.0_f64; n];
+    let mut cum_v = vec![0.0_f64; n];
+    for i in 0..n {
+        cum_cv[i] = close[i] * volume[i] + if i > 0 { cum_cv[i - 1] } else { 0.0 };
+        cum_v[i] = volume[i] + if i > 0 { cum_v[i - 1] } else { 0.0 };
+    }
+
+    for i in (timeperiod - 1)..n {
+        let prev_cv = if i >= timeperiod {
+            cum_cv[i - timeperiod]
+        } else {
+            0.0
+        };
+        let prev_v = if i >= timeperiod {
+            cum_v[i - timeperiod]
+        } else {
+            0.0
+        };
+        let w_cv = cum_cv[i] - prev_cv;
+        let w_v = cum_v[i] - prev_v;
+        result[i] = if w_v != 0.0 { w_cv / w_v } else { f64::NAN };
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// SUPERTREND
+// ---------------------------------------------------------------------------
+
+/// ATR-based Supertrend indicator.
+///
+/// # Returns
+/// `(supertrend_line, direction)` where direction values are:
+/// * `1` = uptrend
+/// * `-1` = downtrend
+/// * `0` = warmup (first `timeperiod` bars)
+pub fn supertrend(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+    multiplier: f64,
+) -> (Vec<f64>, Vec<i8>) {
+    let n = high.len();
+    let mut supertrend_out = vec![f64::NAN; n];
+    let mut direction = vec![0_i8; n];
+
+    if timeperiod < 1 || n <= timeperiod {
+        return (supertrend_out, direction);
+    }
+
+    let atr = compute_atr(high, low, close, timeperiod);
+
+    let mut upper_band = vec![f64::NAN; n];
+    let mut lower_band = vec![f64::NAN; n];
+
+    let first_valid = timeperiod - 1;
+    if first_valid >= n || atr[first_valid].is_nan() {
+        return (supertrend_out, direction);
+    }
+
+    // Initialize band state at first valid ATR bar (compute basic bands inline)
+    {
+        let hl2 = (high[first_valid] + low[first_valid]) / 2.0;
+        upper_band[first_valid] = hl2 + multiplier * atr[first_valid];
+        lower_band[first_valid] = hl2 - multiplier * atr[first_valid];
+    }
+
+    for i in (first_valid + 1)..n {
+        if atr[i].is_nan() {
+            continue;
+        }
+
+        // Compute basic bands as scalars — no Vec allocation needed
+        let hl2 = (high[i] + low[i]) / 2.0;
+        let upper_basic = hl2 + multiplier * atr[i];
+        let lower_basic = hl2 - multiplier * atr[i];
+
+        // Adjust lower band
+        lower_band[i] = if lower_basic > lower_band[i - 1] || close[i - 1] < lower_band[i - 1]
+        {
+            lower_basic
+        } else {
+            lower_band[i - 1]
+        };
+
+        // Adjust upper band
+        upper_band[i] = if upper_basic < upper_band[i - 1] || close[i - 1] > upper_band[i - 1]
+        {
+            upper_basic
+        } else {
+            upper_band[i - 1]
+        };
+
+        // Direction and output only from index timeperiod (warmup = 0, NaN)
+        if i >= timeperiod {
+            let prev_dir = direction[i - 1];
+            direction[i] = if prev_dir == 0 {
+                if close[i] > upper_band[i] {
+                    1
+                } else {
+                    -1
+                }
+            } else if prev_dir == -1 {
+                if close[i] > upper_band[i] {
+                    1
+                } else {
+                    -1
+                }
+            } else if close[i] < lower_band[i] {
+                -1
+            } else {
+                1
+            };
+            supertrend_out[i] = if direction[i] == 1 {
+                lower_band[i]
+            } else {
+                upper_band[i]
+            };
+        }
+    }
+
+    (supertrend_out, direction)
+}
+
+// ---------------------------------------------------------------------------
+// DONCHIAN
+// ---------------------------------------------------------------------------
+
+/// Donchian Channels — rolling highest high / lowest low.
+///
+/// # Returns
+/// `(upper, middle, lower)` arrays.
+pub fn donchian(
+    high: &[f64],
+    low: &[f64],
+    timeperiod: usize,
+) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
+    let n = high.len();
+    let mut upper = vec![f64::NAN; n];
+    let mut lower = vec![f64::NAN; n];
+    let mut middle = vec![f64::NAN; n];
+
+    if timeperiod < 1 || n < timeperiod {
+        return (upper, middle, lower);
+    }
+
+    let hh = math::sliding_max(high, timeperiod);
+    let ll = math::sliding_min(low, timeperiod);
+
+    for i in 0..n {
+        if !hh[i].is_nan() {
+            upper[i] = hh[i];
+            lower[i] = ll[i];
+            middle[i] = (upper[i] + lower[i]) / 2.0;
+        }
+    }
+
+    (upper, middle, lower)
+}
+
+// ---------------------------------------------------------------------------
+// CHOPPINESS_INDEX
+// ---------------------------------------------------------------------------
+
+/// Choppiness Index — measures market choppiness vs trending.
+///
+/// Values near 100 indicate a choppy market; near 0 indicates trending.
+/// The first `timeperiod` values are `NaN`.
+pub fn choppiness_index(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+) -> Vec<f64> {
+    let n = high.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod < 1 || n <= timeperiod {
+        return result;
+    }
+
+    // ATR(1) = True Range per bar
+    let mut tr = vec![0.0_f64; n];
+    tr[0] = high[0] - low[0];
+    for i in 1..n {
+        let hl = high[i] - low[i];
+        let hc = (high[i] - close[i - 1]).abs();
+        let lc = (low[i] - close[i - 1]).abs();
+        tr[i] = hl.max(hc).max(lc);
+    }
+
+    // Cumulative TR for rolling sum
+    let mut cum_tr = vec![0.0_f64; n];
+    cum_tr[0] = tr[0];
+    for i in 1..n {
+        cum_tr[i] = cum_tr[i - 1] + tr[i];
+    }
+
+    let log_n = (timeperiod as f64).log10();
+
+    let hh = math::sliding_max(high, timeperiod);
+    let ll = math::sliding_min(low, timeperiod);
+
+    for i in (timeperiod)..n {
+        let prev_cum = cum_tr[i - timeperiod];
+        let sum_tr = cum_tr[i] - prev_cum;
+        let hl_range = hh[i] - ll[i];
+        if hl_range > 0.0 && log_n > 0.0 {
+            result[i] = 100.0 * (sum_tr / hl_range).log10() / log_n;
+        }
+    }
+
+    result
+}
+
+// ---------------------------------------------------------------------------
+// KELTNER_CHANNELS
+// ---------------------------------------------------------------------------
+
+/// Keltner Channels — EMA +/- (multiplier x ATR).
+///
+/// # Returns
+/// `(upper, middle, lower)` arrays.
+pub fn keltner_channels(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+    atr_period: usize,
+    multiplier: f64,
+) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
+    let n = high.len();
+    if timeperiod < 1 || atr_period < 1 || n < timeperiod || n < atr_period {
+        let nan = vec![f64::NAN; n];
+        return (nan.clone(), nan.clone(), nan);
+    }
+
+    let middle = overlap::ema(close, timeperiod);
+    let atr = compute_atr(high, low, close, atr_period);
+
+    let mut upper = vec![f64::NAN; n];
+    let mut lower = vec![f64::NAN; n];
+    for i in 0..n {
+        if !middle[i].is_nan() && !atr[i].is_nan() {
+            let band = multiplier * atr[i];
+            upper[i] = middle[i] + band;
+            lower[i] = middle[i] - band;
+        }
+    }
+
+    (upper, middle, lower)
+}
+
+// ---------------------------------------------------------------------------
+// HULL_MA
+// ---------------------------------------------------------------------------
+
+/// Hull Moving Average (HMA).
+///
+/// `HMA(n) = WMA(2 * WMA(n/2) - WMA(n), sqrt(n))`
+pub fn hull_ma(close: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = close.len();
+    if timeperiod < 1 || n < timeperiod {
+        return vec![f64::NAN; n];
+    }
+
+    let half = (timeperiod / 2).max(1);
+    let sqrt_p = ((timeperiod as f64).sqrt().round() as usize).max(1);
+
+    let wma_full = overlap::wma(close, timeperiod);
+    let wma_half = overlap::wma(close, half);
+
+    // raw = 2 * wma_half - wma_full
+    let mut raw = vec![f64::NAN; n];
+    for i in 0..n {
+        if !wma_full[i].is_nan() && !wma_half[i].is_nan() {
+            raw[i] = 2.0 * wma_half[i] - wma_full[i];
+        }
+    }
+
+    // Find first valid index in raw
+    let first_valid = raw.iter().position(|x| !x.is_nan()).unwrap_or(n);
+    let mut hull = vec![f64::NAN; n];
+    if first_valid < n {
+        let raw_valid = &raw[first_valid..];
+        let hma_slice = overlap::wma(raw_valid, sqrt_p);
+        for (k, &v) in hma_slice.iter().enumerate() {
+            hull[first_valid + k] = v;
+        }
+    }
+
+    hull
+}
+
+// ---------------------------------------------------------------------------
+// CHANDELIER_EXIT
+// ---------------------------------------------------------------------------
+
+/// Chandelier Exit — ATR-based trailing stop levels.
+///
+/// # Returns
+/// `(long_exit, short_exit)` arrays.
+pub fn chandelier_exit(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    timeperiod: usize,
+    multiplier: f64,
+) -> (Vec<f64>, Vec<f64>) {
+    let n = high.len();
+    if timeperiod < 1 || n < timeperiod {
+        return (vec![f64::NAN; n], vec![f64::NAN; n]);
+    }
+
+    let atr = compute_atr(high, low, close, timeperiod);
+
+    let highest_high = math::sliding_max(high, timeperiod);
+    let lowest_low = math::sliding_min(low, timeperiod);
+
+    let mut long_exit = vec![f64::NAN; n];
+    let mut short_exit = vec![f64::NAN; n];
+    for i in 0..n {
+        if !highest_high[i].is_nan() && !atr[i].is_nan() {
+            long_exit[i] = highest_high[i] - multiplier * atr[i];
+            short_exit[i] = lowest_low[i] + multiplier * atr[i];
+        }
+    }
+
+    (long_exit, short_exit)
+}
+
+// ---------------------------------------------------------------------------
+// ICHIMOKU
+// ---------------------------------------------------------------------------
+
+/// Ichimoku Cloud (Ichimoku Kinko Hyo).
+///
+/// # Returns
+/// `(tenkan, kijun, senkou_a, senkou_b, chikou)` arrays.
+pub fn ichimoku(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    tenkan_period: usize,
+    kijun_period: usize,
+    senkou_b_period: usize,
+    displacement: usize,
+) -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
+    let n = high.len();
+    let nan = || vec![f64::NAN; n];
+
+    if tenkan_period < 1 || kijun_period < 1 || senkou_b_period < 1 {
+        return (nan(), nan(), nan(), nan(), nan());
+    }
+
+    // Helper: rolling (H+L)/2 via shared sliding_max / sliding_min
+    let midpoint_rolling = |period: usize| -> Vec<f64> {
+        let hh = math::sliding_max(high, period);
+        let ll = math::sliding_min(low, period);
+        let mut result = vec![f64::NAN; n];
+        for i in 0..n {
+            if !hh[i].is_nan() {
+                result[i] = (hh[i] + ll[i]) / 2.0;
+            }
+        }
+        result
+    };
+
+    let tenkan = midpoint_rolling(tenkan_period);
+    let kijun = midpoint_rolling(kijun_period);
+    let raw_b = midpoint_rolling(senkou_b_period);
+
+    // Senkou A: (tenkan + kijun) / 2 shifted back `displacement` bars
+    let mut senkou_a = vec![f64::NAN; n];
+    if n > displacement {
+        for i in displacement..n {
+            if !tenkan[i].is_nan() && !kijun[i].is_nan() {
+                senkou_a[i - displacement] = (tenkan[i] + kijun[i]) / 2.0;
+            }
+        }
+    }
+
+    // Senkou B: raw_b shifted back `displacement` bars
+    let mut senkou_b = vec![f64::NAN; n];
+    if n > displacement {
+        senkou_b[..n - displacement].copy_from_slice(&raw_b[displacement..]);
+    }
+
+    // Chikou: close shifted forward `displacement` bars
+    let mut chikou = vec![f64::NAN; n];
+    if n > displacement {
+        chikou[displacement..].copy_from_slice(&close[..n - displacement]);
+    }
+
+    (tenkan, kijun, senkou_a, senkou_b, chikou)
+}
+
+// ---------------------------------------------------------------------------
+// PIVOT_POINTS
+// ---------------------------------------------------------------------------
+
+/// Pivot Points — support / resistance levels computed from the previous bar.
+///
+/// # Arguments
+/// * `method` — `"classic"`, `"fibonacci"`, or `"camarilla"`. Returns all-NaN
+///   vectors for unknown methods.
+///
+/// # Returns
+/// `(pivot, r1, s1, r2, s2)` arrays. Index 0 is always `NaN` (no previous bar).
+pub fn pivot_points(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    method: &str,
+) -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
+    let n = high.len();
+    let mut pivot = vec![f64::NAN; n];
+    let mut r1 = vec![f64::NAN; n];
+    let mut s1 = vec![f64::NAN; n];
+    let mut r2 = vec![f64::NAN; n];
+    let mut s2 = vec![f64::NAN; n];
+
+    let method_lower = method.to_lowercase();
+    if !matches!(method_lower.as_str(), "classic" | "fibonacci" | "camarilla") {
+        // Unknown method — return all NaN
+        return (pivot, r1, s1, r2, s2);
+    }
+
+    for i in 1..n {
+        let ph = high[i - 1];
+        let pl = low[i - 1];
+        let pc = close[i - 1];
+        let hl = ph - pl;
+        let p = (ph + pl + pc) / 3.0;
+        pivot[i] = p;
+        match method_lower.as_str() {
+            "classic" => {
+                r1[i] = 2.0 * p - pl;
+                s1[i] = 2.0 * p - ph;
+                r2[i] = p + hl;
+                s2[i] = p - hl;
+            }
+            "fibonacci" => {
+                r1[i] = p + 0.382 * hl;
+                s1[i] = p - 0.382 * hl;
+                r2[i] = p + 0.618 * hl;
+                s2[i] = p - 0.618 * hl;
+            }
+            "camarilla" => {
+                r1[i] = pc + 1.1 * hl / 12.0;
+                s1[i] = pc - 1.1 * hl / 12.0;
+                r2[i] = pc + 1.1 * hl / 6.0;
+                s2[i] = pc - 1.1 * hl / 6.0;
+            }
+            _ => unreachable!(),
+        }
+    }
+
+    (pivot, r1, s1, r2, s2)
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // Shared test data: 10-bar OHLCV
+    fn sample_ohlcv() -> (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>) {
+        let high = vec![11.0, 12.0, 13.0, 14.0, 15.0, 14.5, 15.5, 16.0, 15.0, 14.0];
+        let low = vec![9.0, 10.0, 11.0, 12.0, 13.0, 12.5, 13.5, 14.0, 13.0, 12.0];
+        let close = vec![10.0, 11.0, 12.0, 13.0, 14.0, 13.5, 14.5, 15.0, 14.0, 13.0];
+        let volume = vec![
+            100.0, 150.0, 200.0, 250.0, 300.0, 200.0, 350.0, 400.0, 180.0, 220.0,
+        ];
+        (high, low, close, volume)
+    }
+
+    // -----------------------------------------------------------------------
+    // VWAP tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn vwap_cumulative_basic() {
+        let (h, l, c, v) = sample_ohlcv();
+        let result = vwap(&h, &l, &c, &v, 0);
+        assert_eq!(result.len(), h.len());
+        // First bar: tp = (11+9+10)/3 = 10.0, tpv = 1000.0, vol = 100.0 => 10.0
+        assert!((result[0] - 10.0).abs() < 1e-10);
+        // All values should be non-NaN for cumulative
+        for val in &result {
+            assert!(!val.is_nan());
+        }
+    }
+
+    #[test]
+    fn vwap_empty_input() {
+        let result = vwap(&[], &[], &[], &[], 0);
+        assert!(result.is_empty());
+    }
+
+    #[test]
+    fn vwap_rolling_basic() {
+        let (h, l, c, v) = sample_ohlcv();
+        let result = vwap(&h, &l, &c, &v, 3);
+        assert_eq!(result.len(), h.len());
+        // First 2 values should be NaN
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        // From index 2 onward should be valid
+        assert!(!result[2].is_nan());
+    }
+
+    // -----------------------------------------------------------------------
+    // VWMA tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn vwma_basic() {
+        let (_, _, c, v) = sample_ohlcv();
+        let result = vwma(&c, &v, 3);
+        assert_eq!(result.len(), c.len());
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        // Index 2: sum(c*v, 0..3) / sum(v, 0..3) = (1000+1650+2400)/(100+150+200) = 5050/450
+        let expected = (10.0 * 100.0 + 11.0 * 150.0 + 12.0 * 200.0) / (100.0 + 150.0 + 200.0);
+        assert!((result[2] - expected).abs() < 1e-10);
+    }
+
+    #[test]
+    fn vwma_empty_input() {
+        let result = vwma(&[], &[], 3);
+        assert!(result.is_empty());
+    }
+
+    #[test]
+    fn vwma_period_larger_than_data() {
+        let result = vwma(&[1.0, 2.0], &[100.0, 200.0], 5);
+        assert_eq!(result.len(), 2);
+        assert!(result.iter().all(|v| v.is_nan()));
+    }
+
+    // -----------------------------------------------------------------------
+    // SUPERTREND tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn supertrend_basic() {
+        let (h, l, c, _) = sample_ohlcv();
+        let (st, dir) = supertrend(&h, &l, &c, 3, 2.0);
+        assert_eq!(st.len(), h.len());
+        assert_eq!(dir.len(), h.len());
+        // First 3 bars should be warmup (direction = 0, st = NaN)
+        for i in 0..3 {
+            assert_eq!(dir[i], 0);
+            assert!(st[i].is_nan());
+        }
+        // From bar 3 onward, direction should be 1 or -1
+        for i in 3..h.len() {
+            assert!(dir[i] == 1 || dir[i] == -1);
+            assert!(!st[i].is_nan());
+        }
+    }
+
+    #[test]
+    fn supertrend_empty_input() {
+        let (st, dir) = supertrend(&[], &[], &[], 3, 2.0);
+        assert!(st.is_empty());
+        assert!(dir.is_empty());
+    }
+
+    #[test]
+    fn supertrend_insufficient_data() {
+        let (st, dir) = supertrend(&[1.0, 2.0], &[0.5, 1.5], &[1.5, 1.8], 5, 2.0);
+        assert!(st.iter().all(|v| v.is_nan()));
+        assert!(dir.iter().all(|&d| d == 0));
+    }
+
+    // -----------------------------------------------------------------------
+    // DONCHIAN tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn donchian_basic() {
+        let (h, l, _, _) = sample_ohlcv();
+        let (upper, middle, lower) = donchian(&h, &l, 3);
+        assert_eq!(upper.len(), h.len());
+        // First 2 are NaN
+        assert!(upper[0].is_nan());
+        assert!(upper[1].is_nan());
+        // Index 2: max(11,12,13)=13, min(9,10,11)=9
+        assert!((upper[2] - 13.0).abs() < 1e-10);
+        assert!((lower[2] - 9.0).abs() < 1e-10);
+        assert!((middle[2] - 11.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn donchian_empty_input() {
+        let (u, m, l) = donchian(&[], &[], 3);
+        assert!(u.is_empty());
+        assert!(m.is_empty());
+        assert!(l.is_empty());
+    }
+
+    #[test]
+    fn donchian_period_1() {
+        let h = vec![5.0, 3.0, 7.0];
+        let l = vec![2.0, 1.0, 4.0];
+        let (upper, middle, lower) = donchian(&h, &l, 1);
+        // Every bar is its own window
+        assert!((upper[0] - 5.0).abs() < 1e-10);
+        assert!((lower[0] - 2.0).abs() < 1e-10);
+        assert!((middle[0] - 3.5).abs() < 1e-10);
+    }
+
+    // -----------------------------------------------------------------------
+    // CHOPPINESS_INDEX tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn choppiness_index_basic() {
+        let (h, l, c, _) = sample_ohlcv();
+        let result = choppiness_index(&h, &l, &c, 3);
+        assert_eq!(result.len(), h.len());
+        // First 3 values should be NaN (timeperiod=3, i+1 > 3 starts at i=3)
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        assert!(result[2].is_nan());
+        // Index 3 should have a valid value (i+1=4 > 3)
+        assert!(!result[3].is_nan());
+        // CI should be between 0 and 100
+        for val in result.iter().filter(|v| !v.is_nan()) {
+            assert!(*val >= 0.0 && *val <= 100.0);
+        }
+    }
+
+    #[test]
+    fn choppiness_index_empty_input() {
+        let result = choppiness_index(&[], &[], &[], 3);
+        assert!(result.is_empty());
+    }
+
+    // -----------------------------------------------------------------------
+    // KELTNER_CHANNELS tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn keltner_channels_basic() {
+        let (h, l, c, _) = sample_ohlcv();
+        let (upper, middle, lower) = keltner_channels(&h, &l, &c, 3, 3, 1.5);
+        assert_eq!(upper.len(), h.len());
+        // Where both EMA and ATR are valid, upper > middle > lower
+        for i in 0..h.len() {
+            if !upper[i].is_nan() && !lower[i].is_nan() {
+                assert!(upper[i] > middle[i]);
+                assert!(lower[i] < middle[i]);
+            }
+        }
+    }
+
+    #[test]
+    fn keltner_channels_empty_input() {
+        let (u, m, l) = keltner_channels(&[], &[], &[], 3, 3, 1.5);
+        assert!(u.is_empty());
+        assert!(m.is_empty());
+        assert!(l.is_empty());
+    }
+
+    // -----------------------------------------------------------------------
+    // HULL_MA tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn hull_ma_basic() {
+        let prices: Vec<f64> = (1..=20).map(|i| i as f64).collect();
+        let result = hull_ma(&prices, 4);
+        assert_eq!(result.len(), prices.len());
+        // Should have some NaN warmup, then valid values
+        let valid_count = result.iter().filter(|v| !v.is_nan()).count();
+        assert!(valid_count > 0);
+    }
+
+    #[test]
+    fn hull_ma_empty_input() {
+        let result = hull_ma(&[], 4);
+        assert!(result.is_empty());
+    }
+
+    #[test]
+    fn hull_ma_period_larger_than_data() {
+        let result = hull_ma(&[1.0, 2.0], 10);
+        assert!(result.iter().all(|v| v.is_nan()));
+    }
+
+    // -----------------------------------------------------------------------
+    // CHANDELIER_EXIT tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn chandelier_exit_basic() {
+        let (h, l, c, _) = sample_ohlcv();
+        let (long_exit, short_exit) = chandelier_exit(&h, &l, &c, 3, 2.0);
+        assert_eq!(long_exit.len(), h.len());
+        assert_eq!(short_exit.len(), h.len());
+        // Where valid, long_exit should be below highest high
+        for i in 0..h.len() {
+            if !long_exit[i].is_nan() {
+                // long_exit = highest_high - multiplier * atr, should be < max high
+                assert!(long_exit[i] < 20.0); // sanity
+            }
+        }
+    }
+
+    #[test]
+    fn chandelier_exit_empty_input() {
+        let (le, se) = chandelier_exit(&[], &[], &[], 3, 2.0);
+        assert!(le.is_empty());
+        assert!(se.is_empty());
+    }
+
+    // -----------------------------------------------------------------------
+    // ICHIMOKU tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn ichimoku_basic() {
+        // Use a larger dataset for ichimoku
+        let n = 60;
+        let high: Vec<f64> = (0..n).map(|i| 100.0 + i as f64 + 1.0).collect();
+        let low: Vec<f64> = (0..n).map(|i| 100.0 + i as f64 - 1.0).collect();
+        let close: Vec<f64> = (0..n).map(|i| 100.0 + i as f64).collect();
+
+        let (tenkan, kijun, senkou_a, senkou_b, chikou) =
+            ichimoku(&high, &low, &close, 9, 26, 52, 26);
+
+        assert_eq!(tenkan.len(), n);
+        assert_eq!(kijun.len(), n);
+        assert_eq!(senkou_a.len(), n);
+        assert_eq!(senkou_b.len(), n);
+        assert_eq!(chikou.len(), n);
+
+        // Tenkan: period 9, first valid at index 8
+        assert!(tenkan[7].is_nan());
+        assert!(!tenkan[8].is_nan());
+
+        // Kijun: period 26, first valid at index 25
+        assert!(kijun[24].is_nan());
+        assert!(!kijun[25].is_nan());
+
+        // Chikou: close shifted forward by 26 bars
+        assert!(chikou[25].is_nan());
+        assert!(!chikou[26].is_nan());
+        assert!((chikou[26] - close[0]).abs() < 1e-10);
+    }
+
+    #[test]
+    fn ichimoku_empty_input() {
+        let (t, k, sa, sb, ch) = ichimoku(&[], &[], &[], 9, 26, 52, 26);
+        assert!(t.is_empty());
+        assert!(k.is_empty());
+        assert!(sa.is_empty());
+        assert!(sb.is_empty());
+        assert!(ch.is_empty());
+    }
+
+    // -----------------------------------------------------------------------
+    // PIVOT_POINTS tests
+    // -----------------------------------------------------------------------
+
+    #[test]
+    fn pivot_points_classic() {
+        let h = vec![10.0, 12.0, 11.0];
+        let l = vec![8.0, 9.0, 8.5];
+        let c = vec![9.0, 11.0, 10.0];
+        let (pivot, r1, s1, r2, s2) = pivot_points(&h, &l, &c, "classic");
+        assert_eq!(pivot.len(), 3);
+        // Index 0 is NaN
+        assert!(pivot[0].is_nan());
+        // Index 1: prev bar H=10, L=8, C=9 => P=(10+8+9)/3=9.0
+        assert!((pivot[1] - 9.0).abs() < 1e-10);
+        // R1 = 2*P - L = 18 - 8 = 10
+        assert!((r1[1] - 10.0).abs() < 1e-10);
+        // S1 = 2*P - H = 18 - 10 = 8
+        assert!((s1[1] - 8.0).abs() < 1e-10);
+        // R2 = P + (H-L) = 9 + 2 = 11
+        assert!((r2[1] - 11.0).abs() < 1e-10);
+        // S2 = P - (H-L) = 9 - 2 = 7
+        assert!((s2[1] - 7.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn pivot_points_fibonacci() {
+        let h = vec![10.0, 12.0];
+        let l = vec![8.0, 9.0];
+        let c = vec![9.0, 11.0];
+        let (pivot, r1, s1, _, _) = pivot_points(&h, &l, &c, "fibonacci");
+        // Index 1: P = (10+8+9)/3 = 9.0, HL = 2
+        assert!((pivot[1] - 9.0).abs() < 1e-10);
+        assert!((r1[1] - (9.0 + 0.382 * 2.0)).abs() < 1e-10);
+        assert!((s1[1] - (9.0 - 0.382 * 2.0)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn pivot_points_camarilla() {
+        let h = vec![10.0, 12.0];
+        let l = vec![8.0, 9.0];
+        let c = vec![9.0, 11.0];
+        let (pivot, r1, s1, _, _) = pivot_points(&h, &l, &c, "camarilla");
+        assert!((pivot[1] - 9.0).abs() < 1e-10);
+        // R1 = C + 1.1 * HL / 12 = 9 + 1.1*2/12
+        assert!((r1[1] - (9.0 + 1.1 * 2.0 / 12.0)).abs() < 1e-10);
+        assert!((s1[1] - (9.0 - 1.1 * 2.0 / 12.0)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn pivot_points_unknown_method() {
+        let h = vec![10.0, 12.0];
+        let l = vec![8.0, 9.0];
+        let c = vec![9.0, 11.0];
+        let (pivot, r1, s1, r2, s2) = pivot_points(&h, &l, &c, "unknown");
+        assert!(pivot.iter().all(|v| v.is_nan()));
+        assert!(r1.iter().all(|v| v.is_nan()));
+        assert!(s1.iter().all(|v| v.is_nan()));
+        assert!(r2.iter().all(|v| v.is_nan()));
+        assert!(s2.iter().all(|v| v.is_nan()));
+    }
+
+    #[test]
+    fn pivot_points_empty_input() {
+        let (p, r1, s1, r2, s2) = pivot_points(&[], &[], &[], "classic");
+        assert!(p.is_empty());
+        assert!(r1.is_empty());
+        assert!(s1.is_empty());
+        assert!(r2.is_empty());
+        assert!(s2.is_empty());
+    }
+}

crates/ferro_ta_core/src/lib.rs

@@ -26,11 +26,30 @@ assert!((sma[2] - 2.0).abs() < 1e-10);
 ```
 */
 
+pub mod aggregation;
+pub mod alerts;
+pub mod attribution;
+pub mod backtest;
+pub mod batch;
+pub mod chunked;
+pub mod commission;
+pub mod crypto;
+pub mod currency;
+pub mod cycle;
+pub mod extended;
 pub mod futures;
 pub mod math;
+pub mod math_ops;
 pub mod momentum;
 pub mod options;
 pub mod overlap;
+pub mod pattern;
+pub mod portfolio;
+pub mod price_transform;
+pub mod regime;
+pub mod resampling;
+pub mod signals;
 pub mod statistic;
+pub mod streaming;
 pub mod volatility;
 pub mod volume;

crates/ferro_ta_core/src/math.rs

@@ -2,7 +2,14 @@
 
 use std::collections::VecDeque;
 
-/// Rolling sum over `timeperiod` bars.
+/// Compute the rolling sum over `timeperiod` bars.
+///
+/// Returns a `Vec<f64>` of length `n`. The first `timeperiod - 1` values
+/// are `NaN`. Uses an incremental algorithm (add new, subtract old) for O(n).
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn sum(real: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = real.len();
     let mut result = vec![f64::NAN; n];
@@ -18,20 +25,38 @@ pub fn sum(real: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Rolling maximum over `timeperiod` bars — O(n) via monotonic deque.
+/// Compute the rolling maximum over `timeperiod` bars.
+///
+/// Delegates to [`sliding_max`] for O(n) performance via a monotonic deque.
+/// The first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn max(real: &[f64], timeperiod: usize) -> Vec<f64> {
     sliding_max(real, timeperiod)
 }
 
-/// Rolling minimum over `timeperiod` bars — O(n) via monotonic deque.
+/// Compute the rolling minimum over `timeperiod` bars.
+///
+/// Delegates to [`sliding_min`] for O(n) performance via a monotonic deque.
+/// The first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn min(real: &[f64], timeperiod: usize) -> Vec<f64> {
     sliding_min(real, timeperiod)
 }
 
-/// Sliding maximum over `timeperiod` bars — O(n) via monotonic deque.
+/// Compute the sliding maximum over `timeperiod` bars in O(n) time.
 ///
-/// Equivalent to `max` but uses a monotonic deque for O(n) total time.
-/// Leading `timeperiod - 1` values are NaN.
+/// Uses a monotonic decreasing deque so each element is pushed/popped at
+/// most once. The first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn sliding_max(real: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = real.len();
     let mut result = vec![f64::NAN; n];
@@ -56,10 +81,14 @@ pub fn sliding_max(real: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Sliding minimum over `timeperiod` bars — O(n) via monotonic deque.
+/// Compute the sliding minimum over `timeperiod` bars in O(n) time.
 ///
-/// Equivalent to `min` but uses a monotonic deque for O(n) total time.
-/// Leading `timeperiod - 1` values are NaN.
+/// Uses a monotonic increasing deque so each element is pushed/popped at
+/// most once. The first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn sliding_min(real: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = real.len();
     let mut result = vec![f64::NAN; n];

crates/ferro_ta_core/src/math_ops.rs

@@ -0,0 +1,154 @@
+//! Rolling math operators — O(n) sliding window implementations.
+//!
+//! - `rolling_sum`      — rolling sum over `timeperiod` bars (prefix-sum based)
+//! - `rolling_max`      — rolling maximum (O(n) monotonic deque)
+//! - `rolling_min`      — rolling minimum (O(n) monotonic deque)
+//! - `rolling_maxindex` — index of rolling maximum
+//! - `rolling_minindex` — index of rolling minimum
+
+use std::collections::VecDeque;
+
+/// Rolling sum over `timeperiod` bars using a prefix-sum array.
+/// Leading `timeperiod - 1` values are NaN.
+pub fn rolling_sum(real: &[f64], timeperiod: usize) -> Vec<f64> {
+    let n = real.len();
+    let mut result = vec![f64::NAN; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+    let mut cs = vec![0.0f64; n + 1];
+    for i in 0..n {
+        cs[i + 1] = cs[i] + real[i];
+    }
+    for i in (timeperiod - 1)..n {
+        result[i] = cs[i + 1] - cs[i + 1 - timeperiod];
+    }
+    result
+}
+
+/// Rolling maximum over `timeperiod` bars (O(n) monotonic deque).
+/// Delegates to `math::sliding_max`.
+pub fn rolling_max(real: &[f64], timeperiod: usize) -> Vec<f64> {
+    crate::math::sliding_max(real, timeperiod)
+}
+
+/// Rolling minimum over `timeperiod` bars (O(n) monotonic deque).
+/// Delegates to `math::sliding_min`.
+pub fn rolling_min(real: &[f64], timeperiod: usize) -> Vec<f64> {
+    crate::math::sliding_min(real, timeperiod)
+}
+
+/// Index of rolling maximum over `timeperiod` bars.
+/// Returns 0-based index. During warmup the value is `-1`.
+pub fn rolling_maxindex(real: &[f64], timeperiod: usize) -> Vec<i64> {
+    let n = real.len();
+    let mut result = vec![-1i64; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+    let mut dq: VecDeque<usize> = VecDeque::new();
+    for i in 0..n {
+        while dq.front().map(|&j| j + timeperiod <= i).unwrap_or(false) {
+            dq.pop_front();
+        }
+        while dq.back().map(|&j| real[j] <= real[i]).unwrap_or(false) {
+            dq.pop_back();
+        }
+        dq.push_back(i);
+        if i + 1 >= timeperiod {
+            result[i] = *dq.front().unwrap() as i64;
+        }
+    }
+    result
+}
+
+/// Index of rolling minimum over `timeperiod` bars.
+/// Returns 0-based index. During warmup the value is `-1`.
+pub fn rolling_minindex(real: &[f64], timeperiod: usize) -> Vec<i64> {
+    let n = real.len();
+    let mut result = vec![-1i64; n];
+    if timeperiod == 0 || n < timeperiod {
+        return result;
+    }
+    let mut dq: VecDeque<usize> = VecDeque::new();
+    for i in 0..n {
+        while dq.front().map(|&j| j + timeperiod <= i).unwrap_or(false) {
+            dq.pop_front();
+        }
+        while dq.back().map(|&j| real[j] >= real[i]).unwrap_or(false) {
+            dq.pop_back();
+        }
+        dq.push_back(i);
+        if i + 1 >= timeperiod {
+            result[i] = *dq.front().unwrap() as i64;
+        }
+    }
+    result
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_rolling_sum() {
+        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let result = rolling_sum(&data, 3);
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        assert!((result[2] - 6.0).abs() < 1e-10); // 1+2+3
+        assert!((result[3] - 9.0).abs() < 1e-10); // 2+3+4
+        assert!((result[4] - 12.0).abs() < 1e-10); // 3+4+5
+    }
+
+    #[test]
+    fn test_rolling_max() {
+        let data = vec![1.0, 3.0, 2.0, 5.0, 4.0];
+        let result = rolling_max(&data, 3);
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        assert!((result[2] - 3.0).abs() < 1e-10);
+        assert!((result[3] - 5.0).abs() < 1e-10);
+        assert!((result[4] - 5.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_rolling_min() {
+        let data = vec![5.0, 3.0, 4.0, 1.0, 2.0];
+        let result = rolling_min(&data, 3);
+        assert!(result[0].is_nan());
+        assert!(result[1].is_nan());
+        assert!((result[2] - 3.0).abs() < 1e-10);
+        assert!((result[3] - 1.0).abs() < 1e-10);
+        assert!((result[4] - 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_rolling_maxindex() {
+        let data = vec![1.0, 3.0, 2.0, 5.0, 4.0];
+        let result = rolling_maxindex(&data, 3);
+        assert_eq!(result[0], -1);
+        assert_eq!(result[1], -1);
+        assert_eq!(result[2], 1); // max(1,3,2) at index 1
+        assert_eq!(result[3], 3); // max(3,2,5) at index 3
+        assert_eq!(result[4], 3); // max(2,5,4) at index 3
+    }
+
+    #[test]
+    fn test_rolling_minindex() {
+        let data = vec![5.0, 3.0, 4.0, 1.0, 2.0];
+        let result = rolling_minindex(&data, 3);
+        assert_eq!(result[0], -1);
+        assert_eq!(result[1], -1);
+        assert_eq!(result[2], 1); // min(5,3,4) at index 1
+        assert_eq!(result[3], 3); // min(3,4,1) at index 3
+        assert_eq!(result[4], 3); // min(4,1,2) at index 3
+    }
+
+    #[test]
+    fn test_short_input() {
+        let data = vec![1.0, 2.0];
+        let result = rolling_sum(&data, 5);
+        assert!(result.iter().all(|v| v.is_nan()));
+    }
+}

crates/ferro_ta_core/src/momentum.rs

@@ -1,11 +1,15 @@
 //! Momentum indicators.
 
-use crate::math::{sliding_max, sliding_min};
 
-/// Relative Strength Index — TA-Lib compatible Wilder smoothing.
+/// Compute the Relative Strength Index (RSI).
 ///
-/// Seeds avg_gain/avg_loss with SMA of first `timeperiod` changes.
-/// Uses branchless gain/loss split: `gain = diff.max(0.0)`, `loss = (-diff).max(0.0)`.
+/// Returns values in the range `[0, 100]`. Uses Wilder's smoothing method
+/// (TA-Lib compatible), seeding avg_gain/avg_loss with the SMA of the first
+/// `timeperiod` price changes. The first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - Lookback period (typically 14).
 pub fn rsi(close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = close.len();
     let mut result = vec![f64::NAN; n];
@@ -46,7 +50,14 @@ pub fn rsi(close: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Momentum — `close[i] - close[i - timeperiod]`.
+/// Compute the Momentum indicator: `close[i] - close[i - timeperiod]`.
+///
+/// Returns a `Vec<f64>` of length `n`. The first `timeperiod` values are `NaN`.
+/// Positive values indicate upward price movement over the lookback window.
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - Number of bars to look back (must be >= 1).
 pub fn mom(close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = close.len();
     let mut result = vec![f64::NAN; n];
@@ -59,14 +70,21 @@ pub fn mom(close: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Stochastic Oscillator — TA-Lib compatible.
+/// Compute the Stochastic Oscillator (TA-Lib compatible).
 ///
-/// Returns `(slowk, slowd)`.
-///  - Fast %K[i] = 100 * (close[i] - min(low, fastk_period)) / (max(high, fastk_period) - min(low, fastk_period))
-///  - Slow %K = SMA(fast %K, slowk_period)
-///  - Slow %D = SMA(slow %K, slowd_period)
+/// Returns `(slow_k, slow_d)`, both in the range `[0, 100]`.
+///  - Fast %K = 100 * (close - lowest low) / (highest high - lowest low)
+///  - Slow %K = SMA(fast %K, `slowk_period`)
+///  - Slow %D = SMA(slow %K, `slowd_period`)
 ///
-/// Uses O(n) sliding max/min via monotonic deques.
+/// Uses O(n) sliding max/min via monotonic deques. Both outputs are
+/// `NaN`-padded until slow %D becomes valid (TA-Lib convention).
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `fastk_period` - Lookback for highest high / lowest low.
+/// * `slowk_period` - SMA period applied to fast %K.
+/// * `slowd_period` - SMA period applied to slow %K.
 pub fn stoch(
     high: &[f64],
     low: &[f64],
@@ -84,29 +102,38 @@ pub fn stoch(
         return nan_pair();
     }
 
-    let max_h = sliding_max(high, fastk_period);
-    let min_l = sliding_min(low, fastk_period);
-
     let mut slowk = vec![f64::NAN; n];
     let mut slowd = vec![f64::NAN; n];
 
-    // Fast %K is valid from index fastk_period-1 onward.
+    // Fused pass: compute fast %K inline with sliding max/min.
+    // For typical small windows (5-14), inline scan beats VecDeque overhead.
     let fastk_start = fastk_period - 1;
-    let mut fastk_valid = vec![0.0; n - fastk_start];
+    let fk_len = n - fastk_start;
+    let mut fastk_valid = vec![0.0_f64; fk_len];
+
     for i in fastk_start..n {
-        let range = max_h[i] - min_l[i];
+        // Inline sliding max(high) and min(low) over [i - fastk_period + 1 .. i].
+        let win_start = i + 1 - fastk_period;
+        let mut hh = high[win_start];
+        let mut ll = low[win_start];
+        for j in (win_start + 1)..=i {
+            let h = high[j];
+            let l = low[j];
+            if h > hh { hh = h; }
+            if l < ll { ll = l; }
+        }
+        let range = hh - ll;
         fastk_valid[i - fastk_start] = if range != 0.0 {
-            100.0 * (close[i] - min_l[i]) / range
+            100.0 * (close[i] - ll) / range
         } else {
             0.0
         };
     }
 
-    // Slow %K = SMA(fastk_valid, slowk_period); write directly into `slowk` offset by `fastk_start`.
+    // Slow %K = SMA(fastk_valid, slowk_period).
     crate::overlap::sma_into(&fastk_valid, slowk_period, &mut slowk, fastk_start);
 
     // Slow %D = SMA(slowk, slowd_period).
-    // The valid part of slowk starts at `fastk_start + slowk_period - 1`.
     let slowk_valid_start = fastk_start + slowk_period - 1;
     let slowd_valid_start = slowk_valid_start + slowd_period - 1;
 
@@ -249,50 +276,164 @@ fn adx_inner(high: &[f64], low: &[f64], close: &[f64], period: usize) -> AdxInne
     (b_pdm, b_mdm, b_pdi, b_mdi, b_dx, b_adx)
 }
 
-/// Plus Directional Movement (Wilder smoothed). Output length = n (bar 0 is NaN).
-pub fn plus_dm(high: &[f64], low: &[f64], timeperiod: usize) -> Vec<f64> {
+/// Compute all six ADX-family outputs in a single pass.
+///
+/// Returns `(plus_dm, minus_dm, plus_di, minus_di, dx, adx)`.
+/// Use this when you need multiple ADX-family outputs to avoid redundant
+/// computation. All values are in `[0, 100]` except DM which is unbounded.
+/// Warmup: DI/DX valid from index `timeperiod`; ADX from `2 * timeperiod - 1`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period (typically 14).
+pub fn adx_all(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> AdxInnerOutput {
+    adx_inner(high, low, close, timeperiod)
+}
+
+/// Internal helper for plus_dm and minus_dm that doesn't allocate dummy close prices.
+/// Returns (plus_dm, minus_dm) smoothed with Wilder's method.
+fn dm_only_inner(high: &[f64], low: &[f64], period: usize) -> (Vec<f64>, Vec<f64>) {
     let n = high.len();
-    let closes = vec![0.0_f64; n];
-    let (pdm, _, _, _, _, _) = adx_inner(high, low, &closes, timeperiod);
+    let mut b_pdm = vec![f64::NAN; n];
+    let mut b_mdm = vec![f64::NAN; n];
+
+    if n < period || period < 1 || n < 2 {
+        return (b_pdm, b_mdm);
+    }
+
+    let m = n - 1;
+    let mut pdm = vec![0.0_f64; m];
+    let mut mdm = vec![0.0_f64; m];
+
+    for i in 0..m {
+        let j = i + 1;
+        let h_diff = high[j] - high[i];
+        let l_diff = low[i] - low[j];
+        pdm[i] = if h_diff > l_diff && h_diff > 0.0 {
+            h_diff
+        } else {
+            0.0
+        };
+        mdm[i] = if l_diff > h_diff && l_diff > 0.0 {
+            l_diff
+        } else {
+            0.0
+        };
+    }
+
+    if m < period {
+        return (b_pdm, b_mdm);
+    }
+
+    let mut pdm_s = pdm[..period].iter().sum::<f64>();
+    let mut mdm_s = mdm[..period].iter().sum::<f64>();
+
+    b_pdm[period] = pdm_s;
+    b_mdm[period] = mdm_s;
+
+    let decay = (period - 1) as f64 / period as f64;
+    for i in period..m {
+        pdm_s = pdm_s * decay + pdm[i];
+        mdm_s = mdm_s * decay + mdm[i];
+        b_pdm[i + 1] = pdm_s;
+        b_mdm[i + 1] = mdm_s;
+    }
+
+    (b_pdm, b_mdm)
+}
+
+/// Compute the Plus Directional Movement (+DM), Wilder smoothed.
+///
+/// Measures upward price movement. Returns a `Vec<f64>` of length `n`;
+/// the first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` - High and low price series (same length).
+/// * `timeperiod` - Wilder smoothing period.
+pub fn plus_dm(high: &[f64], low: &[f64], timeperiod: usize) -> Vec<f64> {
+    let (pdm, _) = dm_only_inner(high, low, timeperiod);
     pdm
 }
 
-/// Minus Directional Movement (Wilder smoothed). Output length = n (bar 0 is NaN).
+/// Compute the Minus Directional Movement (-DM), Wilder smoothed.
+///
+/// Measures downward price movement. Returns a `Vec<f64>` of length `n`;
+/// the first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` - High and low price series (same length).
+/// * `timeperiod` - Wilder smoothing period.
 pub fn minus_dm(high: &[f64], low: &[f64], timeperiod: usize) -> Vec<f64> {
-    let n = high.len();
-    let closes = vec![0.0_f64; n];
-    let (_, mdm, _, _, _, _) = adx_inner(high, low, &closes, timeperiod);
+    let (_, mdm) = dm_only_inner(high, low, timeperiod);
     mdm
 }
 
-/// Plus Directional Indicator (Wilder smoothed). Output length = n.
+/// Compute the Plus Directional Indicator (+DI), Wilder smoothed.
+///
+/// `+DI = 100 * smoothed(+DM) / smoothed(TR)`. Returns values in `[0, 100]`.
+/// The first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period.
 pub fn plus_di(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let (_, _, pdi, _, _, _) = adx_inner(high, low, close, timeperiod);
     pdi
 }
 
-/// Minus Directional Indicator (Wilder smoothed). Output length = n.
+/// Compute the Minus Directional Indicator (-DI), Wilder smoothed.
+///
+/// `-DI = 100 * smoothed(-DM) / smoothed(TR)`. Returns values in `[0, 100]`.
+/// The first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period.
 pub fn minus_di(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let (_, _, _, mdi, _, _) = adx_inner(high, low, close, timeperiod);
     mdi
 }
 
-/// Directional Movement Index: 100 * |+DI − −DI| / (+DI + −DI).
+/// Compute the Directional Movement Index (DX).
+///
+/// `DX = 100 * |+DI - -DI| / (+DI + -DI)`. Returns values in `[0, 100]`.
+/// The first `timeperiod` values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period.
 pub fn dx(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let (_, _, _, _, dx_vals, _) = adx_inner(high, low, close, timeperiod);
     dx_vals
 }
 
-/// Average Directional Movement Index (Wilder smoothing of DX).
+/// Compute the Average Directional Movement Index (ADX).
+///
+/// ADX is Wilder's smoothing of DX, measuring trend strength regardless of
+/// direction. Returns values in `[0, 100]`. The first `2 * timeperiod - 1`
+/// values are `NaN` (DX warmup + ADX smoothing warmup).
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period (typically 14).
 pub fn adx(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let (_, _, _, _, _, adx_vals) = adx_inner(high, low, close, timeperiod);
     adx_vals
 }
 
-/// ADX Rating: (ADX[i] + ADX[i − timeperiod]) / 2.
+/// Compute the ADX Rating (ADXR).
+///
+/// `ADXR[i] = (ADX[i] + ADX[i - timeperiod]) / 2`. Smooths ADX further
+/// by averaging current ADX with its value `timeperiod` bars ago.
+/// Returns values in `[0, 100]`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Wilder smoothing period (typically 14).
 pub fn adxr(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = high.len();
-    let adx_vals = adx(high, low, close, timeperiod);
+    // Reuse adx_all to compute ADX once, then derive ADXR from it
+    let (_, _, _, _, _, adx_vals) = adx_inner(high, low, close, timeperiod);
     let mut result = vec![f64::NAN; n];
     for i in timeperiod..n {
         if !adx_vals[i].is_nan() && !adx_vals[i - timeperiod].is_nan() {

crates/ferro_ta_core/src/overlap.rs

@@ -3,7 +3,14 @@
 //! All functions return a `Vec<f64>` of the same length as the input.
 //! Leading values are `f64::NAN` for the warm-up period.
 
-/// Simple Moving Average over `timeperiod` bars.
+/// Compute the Simple Moving Average (SMA) over a rolling window.
+///
+/// Returns a `Vec<f64>` of the same length as `close`. The first
+/// `timeperiod - 1` values are `NaN` (warmup period).
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 ///
 /// # Edge Cases
 /// Returns all-NaN when `timeperiod < 1` or `close.len() < timeperiod`.
@@ -14,8 +21,17 @@ pub fn sma(close: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Simple Moving Average written directly into `dest` starting at `dest_offset`.
-/// Leaves values before `dest_offset + timeperiod - 1` untouched (e.g. they can be NaN).
+/// Write a Simple Moving Average directly into a pre-allocated buffer.
+///
+/// Values before `dest_offset + timeperiod - 1` are left untouched.
+/// This avoids an intermediate allocation when composing indicators
+/// (e.g., Stochastic slow %K and slow %D).
+///
+/// # Arguments
+/// * `src` - Input price series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
+/// * `dest` - Output buffer (must be at least `dest_offset + src.len()` long).
+/// * `dest_offset` - Starting index in `dest` to write results.
 pub fn sma_into(src: &[f64], timeperiod: usize, dest: &mut [f64], dest_offset: usize) {
     let n = src.len();
     if timeperiod < 1 || n < timeperiod {
@@ -66,7 +82,15 @@ pub fn sma_into(src: &[f64], timeperiod: usize, dest: &mut [f64], dest_offset: u
     }
 }
 
-/// Exponential Moving Average — seeded with SMA of first `timeperiod` bars.
+/// Compute the Exponential Moving Average (EMA).
+///
+/// The EMA is seeded with the SMA of the first `timeperiod` bars and uses
+/// a smoothing factor of `k = 2 / (timeperiod + 1)`. Returns a `Vec<f64>`
+/// of the same length as `close`; the first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - Lookback period (must be >= 1).
 pub fn ema(close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = close.len();
     let mut result = vec![f64::NAN; n];
@@ -82,10 +106,15 @@ pub fn ema(close: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Weighted Moving Average — O(n) incremental algorithm using running weighted sum.
+/// Compute the Weighted Moving Average (WMA).
 ///
-/// Recurrence: `T[i] = T[i-1] + n*close[i] - S[i-1]`
-/// where `S[i]` is the rolling sum over `timeperiod` bars.
+/// Assigns linearly increasing weights (1, 2, ..., timeperiod) to the window.
+/// Uses an O(n) incremental recurrence to avoid recomputing weights each bar.
+/// Returns a `Vec<f64>` of length `n`; the first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
 pub fn wma(close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = close.len();
     let mut result = vec![f64::NAN; n];
@@ -157,10 +186,43 @@ pub fn wma(close: &[f64], timeperiod: usize) -> Vec<f64> {
     result
 }
 
-/// Bollinger Bands — returns `(upper, middle, lower)`.
+/// Compute Bollinger Bands, returning `(upper, middle, lower)`.
 ///
-/// Middle is SMA; bands are `± nbdev * stddev`.
-/// Uses O(n) sliding `sum` and `sum_sq` windows for mean and variance.
+/// The middle band is the SMA; upper and lower bands are offset by
+/// `nbdevup` and `nbdevdn` standard deviations respectively. Uses
+/// Welford's rolling algorithm for numerically stable variance in O(n).
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `timeperiod` - SMA / standard deviation window (must be >= 1).
+/// * `nbdevup` - Number of standard deviations above the mean for the upper band.
+/// * `nbdevdn` - Number of standard deviations below the mean for the lower band.
+///
+/// # Returns
+/// `(upper, middle, lower)` -- each `Vec<f64>` of length `n`. The first
+/// `timeperiod - 1` values in each vector are `NaN`.
+///
+/// ## Welford's rolling algorithm
+///
+/// We maintain `mean` and `m2` (sum of squared deviations from the current
+/// mean) across a sliding window of size `N`.  When a new value `x_new`
+/// replaces an old value `x_old` (window size stays constant):
+///
+/// ```text
+/// delta     = x_new - x_old
+/// old_mean  = mean
+/// mean     += delta / N
+/// m2       += delta * ((x_new - mean) + (x_old - old_mean))
+///
+/// variance  = m2 / N               // population variance
+/// stddev    = sqrt(variance)
+/// ```
+///
+/// The initial window is seeded using the standard (non-rolling) Welford
+/// incremental algorithm.
+///
+/// This avoids the catastrophic cancellation inherent in the naïve
+/// `Σx²/N − mean²` formula when values are large but close together.
 pub fn bbands(
     close: &[f64],
     timeperiod: usize,
@@ -177,90 +239,135 @@ pub fn bbands(
     let mut lower = vec![f64::NAN; n];
     let p = timeperiod as f64;
 
-    // Seed sliding sums for the first window.
-    #[cfg(feature = "simd")]
-    let (mut sum, mut sum_sq) = {
-        use wide::f64x4;
-        let p_data = &close[..timeperiod];
-        let mut sum_simd = f64x4::splat(0.0);
-        let mut sq_simd = f64x4::splat(0.0);
-        let mut chunks = p_data.chunks_exact(4);
-        for chunk in &mut chunks {
-            let vals = f64x4::new([chunk[0], chunk[1], chunk[2], chunk[3]]);
-            sum_simd += vals;
-            sq_simd += vals * vals;
-        }
-        let s_arr = sum_simd.to_array();
-        let sq_arr = sq_simd.to_array();
-        let mut sum = s_arr[0] + s_arr[1] + s_arr[2] + s_arr[3];
-        let mut sum_sq = sq_arr[0] + sq_arr[1] + sq_arr[2] + sq_arr[3];
-        for &v in chunks.remainder() {
-            sum += v;
-            sum_sq += v * v;
-        }
-        (sum, sum_sq)
-    };
+    // --- Seed: build initial mean and m2 for the first window using
+    //     Welford's incremental (non-rolling) algorithm. ---
+    let mut mean = 0.0_f64;
+    let mut m2 = 0.0_f64;
+    for (k, &x) in close[..timeperiod].iter().enumerate() {
+        let count = (k + 1) as f64;
+        let delta = x - mean;
+        mean += delta / count;
+        let delta2 = x - mean;
+        m2 += delta * delta2;
+    }
 
-    #[cfg(not(feature = "simd"))]
-    let (mut sum, mut sum_sq) = {
-        let s: f64 = close[..timeperiod].iter().sum();
-        let sq: f64 = close[..timeperiod].iter().map(|&x| x * x).sum();
-        (s, sq)
-    };
-
-    let mean = sum / p;
-    let var = (sum_sq / p - mean * mean).max(0.0);
+    let var = (m2 / p).max(0.0);
     let std = var.sqrt();
     middle[timeperiod - 1] = mean;
     upper[timeperiod - 1] = mean + nbdevup * std;
     lower[timeperiod - 1] = mean - nbdevdn * std;
 
+    // --- Rolling phase: slide the window one element at a time,
+    //     removing the oldest value and adding the newest. ---
+
+    /// Inline helper: replace `x_old` with `x_new` in the Welford accumulator
+    /// (constant window size `p`), then write band values into the output slots.
+    ///
+    /// Combined rolling Welford update (window size stays constant at N):
+    ///
+    /// ```text
+    /// delta     = x_new - x_old
+    /// old_mean  = mean
+    /// mean     += delta / N
+    /// m2       += delta * ((x_new - mean) + (x_old - old_mean))
+    /// ```
+    ///
+    /// This is algebraically equivalent to removing `x_old` and adding `x_new`
+    /// in two separate Welford steps, but avoids the intermediate N-1 state.
+    #[inline(always)]
+    #[allow(clippy::too_many_arguments)]
+    fn welford_step(
+        x_old: f64,
+        x_new: f64,
+        mean: &mut f64,
+        m2: &mut f64,
+        p: f64,
+        nbdevup: f64,
+        nbdevdn: f64,
+        upper: &mut f64,
+        middle: &mut f64,
+        lower: &mut f64,
+    ) {
+        let delta = x_new - x_old;
+        let old_mean = *mean;
+        *mean += delta / p;
+        // Update m2 using both old and new deviations.
+        *m2 += delta * ((x_new - *mean) + (x_old - old_mean));
+
+        // Clamp m2 to zero to guard against floating-point drift.
+        if *m2 < 0.0 {
+            *m2 = 0.0;
+        }
+
+        let var = *m2 / p;
+        let std = var.sqrt();
+        *middle = *mean;
+        *upper = *mean + nbdevup * std;
+        *lower = *mean - nbdevdn * std;
+    }
+
+    // Process two iterations at a time (loop unrolling) for throughput.
     let mut i = timeperiod;
     while i + 1 < n {
-        let old0 = close[i - timeperiod];
-        sum += close[i] - old0;
-        sum_sq += close[i] * close[i] - old0 * old0;
-        let mean = sum / p;
-        let var = (sum_sq / p - mean * mean).max(0.0);
-        let std = var.sqrt();
-        middle[i] = mean;
-        upper[i] = mean + nbdevup * std;
-        lower[i] = mean - nbdevdn * std;
-
-        let old1 = close[i + 1 - timeperiod];
-        sum += close[i + 1] - old1;
-        sum_sq += close[i + 1] * close[i + 1] - old1 * old1;
-        let mean1 = sum / p;
-        let var1 = (sum_sq / p - mean1 * mean1).max(0.0);
-        let std1 = var1.sqrt();
-        middle[i + 1] = mean1;
-        upper[i + 1] = mean1 + nbdevup * std1;
-        lower[i + 1] = mean1 - nbdevdn * std1;
-
+        welford_step(
+            close[i - timeperiod],
+            close[i],
+            &mut mean,
+            &mut m2,
+            p,
+            nbdevup,
+            nbdevdn,
+            &mut upper[i],
+            &mut middle[i],
+            &mut lower[i],
+        );
+        welford_step(
+            close[i + 1 - timeperiod],
+            close[i + 1],
+            &mut mean,
+            &mut m2,
+            p,
+            nbdevup,
+            nbdevdn,
+            &mut upper[i + 1],
+            &mut middle[i + 1],
+            &mut lower[i + 1],
+        );
         i += 2;
     }
     if i < n {
-        let old = close[i - timeperiod];
-        sum += close[i] - old;
-        sum_sq += close[i] * close[i] - old * old;
-        let mean = sum / p;
-        let var = (sum_sq / p - mean * mean).max(0.0);
-        let std = var.sqrt();
-        middle[i] = mean;
-        upper[i] = mean + nbdevup * std;
-        lower[i] = mean - nbdevdn * std;
+        welford_step(
+            close[i - timeperiod],
+            close[i],
+            &mut mean,
+            &mut m2,
+            p,
+            nbdevup,
+            nbdevdn,
+            &mut upper[i],
+            &mut middle[i],
+            &mut lower[i],
+        );
     }
+
     (upper, middle, lower)
 }
 
-/// MACD — EMA(fastperiod) minus EMA(slowperiod), signal = EMA(macd, signalperiod).
+/// Compute the Moving Average Convergence/Divergence (MACD).
 ///
-/// Returns `(macd_line, signal_line, histogram)`, each of length `n`.
-/// Leading values are `NaN` during warmup.
-/// `fastperiod` must be less than `slowperiod`.
+/// `MACD = EMA(close, fastperiod) - EMA(close, slowperiod)`.
+/// The signal line is `EMA(macd, signalperiod)` and the histogram is
+/// `macd - signal`. TA-Lib compatible: leading values are `NaN` up to
+/// the point where all three outputs are valid.
 ///
-/// Fast and slow EMAs are computed in a **single combined loop** to minimise
-/// memory round-trips, then the signal EMA is computed in a second pass.
+/// # Arguments
+/// * `close` - Price series.
+/// * `fastperiod` - Fast EMA period (must be < `slowperiod`).
+/// * `slowperiod` - Slow EMA period.
+/// * `signalperiod` - Signal line EMA period.
+///
+/// # Returns
+/// `(macd_line, signal_line, histogram)` -- each `Vec<f64>` of length `n`.
 pub fn macd(
     close: &[f64],
     fastperiod: usize,
@@ -383,6 +490,74 @@ mod tests {
         assert!((lower[2] - 2.0).abs() < 1e-10);
     }
 
+    #[test]
+    fn bbands_varying_prices() {
+        // Verify against hand-computed values for a small window.
+        let prices = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let (upper, middle, lower) = bbands(&prices, 3, 2.0, 2.0);
+
+        // First two values should be NaN (warmup).
+        assert!(middle[0].is_nan());
+        assert!(middle[1].is_nan());
+
+        // Window [1,2,3]: mean = 2.0, pop_var = 2/3, std = sqrt(2/3)
+        let expected_mean = 2.0;
+        let expected_std = (2.0_f64 / 3.0).sqrt();
+        assert!((middle[2] - expected_mean).abs() < 1e-10);
+        assert!((upper[2] - (expected_mean + 2.0 * expected_std)).abs() < 1e-10);
+        assert!((lower[2] - (expected_mean - 2.0 * expected_std)).abs() < 1e-10);
+
+        // Window [2,3,4]: mean = 3.0, pop_var = 2/3, std = sqrt(2/3)
+        assert!((middle[3] - 3.0).abs() < 1e-10);
+        assert!((upper[3] - (3.0 + 2.0 * expected_std)).abs() < 1e-10);
+
+        // Window [3,4,5]: mean = 4.0, pop_var = 2/3, std = sqrt(2/3)
+        assert!((middle[4] - 4.0).abs() < 1e-10);
+        assert!((upper[4] - (4.0 + 2.0 * expected_std)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn bbands_numerical_stability() {
+        // Large offset with tiny variation — this is where the naïve sum_sq
+        // formula suffers from catastrophic cancellation.
+        let base = 1e12;
+        let prices: Vec<f64> = (0..100).map(|i| base + (i as f64) * 0.01).collect();
+        let (upper, middle, lower) = bbands(&prices, 20, 2.0, 2.0);
+
+        // Check that middle band matches SMA.
+        for i in 19..100 {
+            let window = &prices[i - 19..=i];
+            let expected_mean: f64 = window.iter().sum::<f64>() / 20.0;
+            assert!(
+                (middle[i] - expected_mean).abs() < 1e-4,
+                "mean mismatch at {i}: got {} expected {}",
+                middle[i],
+                expected_mean,
+            );
+            // Bands should be above/below middle.
+            assert!(upper[i] >= middle[i]);
+            assert!(lower[i] <= middle[i]);
+        }
+    }
+
+    #[test]
+    fn bbands_edge_cases() {
+        // timeperiod == 1: every bar should have std = 0, bands == price.
+        let prices = vec![10.0, 20.0, 30.0];
+        let (upper, middle, lower) = bbands(&prices, 1, 2.0, 2.0);
+        for i in 0..3 {
+            assert!((middle[i] - prices[i]).abs() < 1e-10);
+            assert!((upper[i] - prices[i]).abs() < 1e-10);
+            assert!((lower[i] - prices[i]).abs() < 1e-10);
+        }
+
+        // Input shorter than timeperiod: all NaN.
+        let (u, m, l) = bbands(&[1.0, 2.0], 5, 2.0, 2.0);
+        assert!(u.iter().all(|v| v.is_nan()));
+        assert!(m.iter().all(|v| v.is_nan()));
+        assert!(l.iter().all(|v| v.is_nan()));
+    }
+
     #[test]
     fn macd_basic() {
         // 40 bars of linearly increasing prices — MACD line should converge

crates/ferro_ta_core/src/portfolio.rs

@@ -0,0 +1,627 @@
+//! Pure Rust portfolio analytics — no PyO3, no numpy, no ndarray.
+//!
+//! Functions:
+//! - `portfolio_volatility`  — sqrt(w' Σ w)
+//! - `beta_full`             — Cov/Var OLS beta
+//! - `rolling_beta`          — rolling beta with NaN warmup
+//! - `drawdown_series`       — per-bar drawdown + max drawdown
+//! - `correlation_matrix`    — pairwise Pearson correlation
+//! - `relative_strength`     — cumulative return ratio
+//! - `spread`                — A - hedge * B
+//! - `ratio`                 — A / B (NaN for zero)
+//! - `zscore_series`         — rolling z-score, NaN warmup
+//! - `compose_weighted`      — weighted sum per row
+
+// ---------------------------------------------------------------------------
+// portfolio_volatility
+// ---------------------------------------------------------------------------
+
+/// Compute portfolio volatility: sqrt(w' Σ w).
+///
+/// `cov_matrix` is an n×n covariance matrix stored as a slice of row-Vecs.
+/// `weights` has length n.
+///
+/// Panics if dimensions are inconsistent.
+pub fn portfolio_volatility(cov_matrix: &[Vec<f64>], weights: &[f64]) -> f64 {
+    let n = weights.len();
+    assert!(
+        cov_matrix.len() == n,
+        "cov_matrix must have {} rows, got {}",
+        n,
+        cov_matrix.len()
+    );
+    let mut variance = 0.0_f64;
+    for i in 0..n {
+        assert!(
+            cov_matrix[i].len() == n,
+            "cov_matrix row {} must have length {}, got {}",
+            i,
+            n,
+            cov_matrix[i].len()
+        );
+        let mut row_sum = 0.0_f64;
+        for j in 0..n {
+            row_sum += weights[j] * cov_matrix[i][j];
+        }
+        variance += weights[i] * row_sum;
+    }
+    variance.max(0.0).sqrt()
+}
+
+// ---------------------------------------------------------------------------
+// beta_full
+// ---------------------------------------------------------------------------
+
+/// Compute the full-sample OLS beta of `asset_returns` vs `benchmark_returns`.
+///
+/// Beta = Cov(asset, bench) / Var(bench).
+///
+/// Panics if lengths differ or are < 2, or if benchmark has zero variance.
+pub fn beta_full(asset_returns: &[f64], benchmark_returns: &[f64]) -> f64 {
+    let n = asset_returns.len();
+    assert!(
+        n >= 2 && benchmark_returns.len() == n,
+        "asset_returns and benchmark_returns must have equal length >= 2"
+    );
+    let mean_a: f64 = asset_returns.iter().sum::<f64>() / n as f64;
+    let mean_b: f64 = benchmark_returns.iter().sum::<f64>() / n as f64;
+    let mut cov = 0.0_f64;
+    let mut var_b = 0.0_f64;
+    for i in 0..n {
+        let da = asset_returns[i] - mean_a;
+        let db = benchmark_returns[i] - mean_b;
+        cov += da * db;
+        var_b += db * db;
+    }
+    assert!(var_b != 0.0, "benchmark_returns has zero variance; cannot compute beta");
+    cov / var_b
+}
+
+// ---------------------------------------------------------------------------
+// rolling_beta
+// ---------------------------------------------------------------------------
+
+/// Compute rolling beta of `asset` vs `benchmark` over a sliding `window`.
+///
+/// Returns a Vec of the same length as the inputs. The first `window - 1`
+/// entries are NaN (warmup period). `window` must be >= 2.
+pub fn rolling_beta(asset: &[f64], benchmark: &[f64], window: usize) -> Vec<f64> {
+    assert!(window >= 2, "window must be >= 2");
+    let n = asset.len();
+    assert!(
+        n > 0 && benchmark.len() == n,
+        "asset and benchmark must be non-empty and equal length"
+    );
+    let mut result = vec![f64::NAN; n];
+    for i in (window - 1)..n {
+        let start = i + 1 - window;
+        let a_win = &asset[start..=i];
+        let b_win = &benchmark[start..=i];
+        let mean_a: f64 = a_win.iter().sum::<f64>() / window as f64;
+        let mean_b: f64 = b_win.iter().sum::<f64>() / window as f64;
+        let mut cov = 0.0_f64;
+        let mut var_b = 0.0_f64;
+        for k in 0..window {
+            let da = a_win[k] - mean_a;
+            let db = b_win[k] - mean_b;
+            cov += da * db;
+            var_b += db * db;
+        }
+        result[i] = if var_b == 0.0 { f64::NAN } else { cov / var_b };
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// drawdown_series
+// ---------------------------------------------------------------------------
+
+/// Compute the drawdown series and maximum drawdown for an equity/price series.
+///
+/// Drawdown at bar i = (equity[i] - running_max) / running_max  (always <= 0).
+///
+/// Returns `(dd_array, max_dd)` where `max_dd` is the most negative drawdown.
+///
+/// Panics if `equity` is empty.
+pub fn drawdown_series(equity: &[f64]) -> (Vec<f64>, f64) {
+    let n = equity.len();
+    assert!(n > 0, "equity must be non-empty");
+    let mut dd = vec![0.0_f64; n];
+    let mut peak = equity[0];
+    let mut max_dd = 0.0_f64;
+    for i in 0..n {
+        if equity[i] > peak {
+            peak = equity[i];
+        }
+        let d = if peak == 0.0 {
+            0.0
+        } else {
+            (equity[i] - peak) / peak
+        };
+        dd[i] = d;
+        if d < max_dd {
+            max_dd = d;
+        }
+    }
+    (dd, max_dd)
+}
+
+// ---------------------------------------------------------------------------
+// correlation_matrix
+// ---------------------------------------------------------------------------
+
+/// Compute the pairwise Pearson correlation matrix.
+///
+/// `data` is a slice of column vectors — `data[j]` is the return series for
+/// asset j, so `data[j][i]` is the return of asset j at bar i. All columns
+/// must have the same length (>= 2).
+///
+/// Returns an n_assets × n_assets matrix stored as `Vec<Vec<f64>>`.
+pub fn correlation_matrix(data: &[Vec<f64>]) -> Vec<Vec<f64>> {
+    let n_assets = data.len();
+    assert!(n_assets > 0, "data must contain at least one asset column");
+    let n_bars = data[0].len();
+    assert!(n_bars >= 2, "data must have at least 2 rows (bars)");
+    for j in 1..n_assets {
+        assert!(
+            data[j].len() == n_bars,
+            "all columns must have equal length; column 0 has {} but column {} has {}",
+            n_bars,
+            j,
+            data[j].len()
+        );
+    }
+
+    // Means
+    let mut means = vec![0.0_f64; n_assets];
+    for j in 0..n_assets {
+        means[j] = data[j].iter().sum::<f64>() / n_bars as f64;
+    }
+
+    // Standard deviations (population)
+    let mut stds = vec![0.0_f64; n_assets];
+    for j in 0..n_assets {
+        let var: f64 = data[j].iter().map(|&v| (v - means[j]).powi(2)).sum::<f64>() / n_bars as f64;
+        stds[j] = var.sqrt();
+    }
+
+    // Build correlation matrix (exploit symmetry: compute each pair once)
+    let mut result = vec![vec![0.0_f64; n_assets]; n_assets];
+    for j1 in 0..n_assets {
+        result[j1][j1] = 1.0;
+        for j2 in (j1 + 1)..n_assets {
+            let mut cov = 0.0_f64;
+            for i in 0..n_bars {
+                cov += (data[j1][i] - means[j1]) * (data[j2][i] - means[j2]);
+            }
+            cov /= n_bars as f64;
+            let denom = stds[j1] * stds[j2];
+            let corr = if denom == 0.0 { f64::NAN } else { cov / denom };
+            result[j1][j2] = corr;
+            result[j2][j1] = corr;
+        }
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// relative_strength
+// ---------------------------------------------------------------------------
+
+/// Compute relative strength of an asset vs a benchmark.
+///
+/// result[i] = cumprod(1 + asset_returns[0..=i]) / cumprod(1 + benchmark_returns[0..=i])
+///
+/// Panics if lengths differ or are zero.
+pub fn relative_strength(asset_returns: &[f64], benchmark_returns: &[f64]) -> Vec<f64> {
+    let n = asset_returns.len();
+    assert!(
+        n > 0 && benchmark_returns.len() == n,
+        "asset_returns and benchmark_returns must be non-empty and equal length"
+    );
+    let mut result = vec![0.0_f64; n];
+    let mut cum_a = 1.0_f64;
+    let mut cum_b = 1.0_f64;
+    for i in 0..n {
+        cum_a *= 1.0 + asset_returns[i];
+        cum_b *= 1.0 + benchmark_returns[i];
+        result[i] = if cum_b == 0.0 { f64::NAN } else { cum_a / cum_b };
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// spread
+// ---------------------------------------------------------------------------
+
+/// Compute the spread between two series: a - hedge * b.
+///
+/// Panics if lengths differ or are zero.
+pub fn spread(a: &[f64], b: &[f64], hedge: f64) -> Vec<f64> {
+    let n = a.len();
+    assert!(
+        n > 0 && b.len() == n,
+        "a and b must be non-empty and equal length"
+    );
+    a.iter().zip(b.iter()).map(|(&x, &y)| x - hedge * y).collect()
+}
+
+// ---------------------------------------------------------------------------
+// ratio
+// ---------------------------------------------------------------------------
+
+/// Compute the ratio between two series: a / b.
+///
+/// Where b is 0, returns NaN.
+///
+/// Panics if lengths differ or are zero.
+pub fn ratio(a: &[f64], b: &[f64]) -> Vec<f64> {
+    let n = a.len();
+    assert!(
+        n > 0 && b.len() == n,
+        "a and b must be non-empty and equal length"
+    );
+    a.iter()
+        .zip(b.iter())
+        .map(|(&x, &y)| if y == 0.0 { f64::NAN } else { x / y })
+        .collect()
+}
+
+// ---------------------------------------------------------------------------
+// zscore_series
+// ---------------------------------------------------------------------------
+
+/// Compute the rolling Z-score of a 1-D series.
+///
+/// Z[i] = (x[i] - mean(window)) / std(window)
+///
+/// The first `window - 1` entries are NaN. `window` must be >= 2.
+///
+/// Panics if `x` is empty or `window < 2`.
+pub fn zscore_series(x: &[f64], window: usize) -> Vec<f64> {
+    assert!(window >= 2, "window must be >= 2");
+    let n = x.len();
+    assert!(n > 0, "x must be non-empty");
+    let mut result = vec![f64::NAN; n];
+    for i in (window - 1)..n {
+        let win = &x[i + 1 - window..=i];
+        let mean: f64 = win.iter().sum::<f64>() / window as f64;
+        let var: f64 = win.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / window as f64;
+        let std = var.sqrt();
+        result[i] = if std == 0.0 {
+            f64::NAN
+        } else {
+            (x[i] - mean) / std
+        };
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// compose_weighted
+// ---------------------------------------------------------------------------
+
+/// Weighted combination of multiple signal columns.
+///
+/// `data` is a slice of column vectors — `data[j]` is one signal column.
+/// `weights` has one entry per column.
+///
+/// Returns a Vec of length n_bars where each entry is the weighted sum across
+/// columns for that bar.
+///
+/// Panics if weights length != number of columns, or columns have unequal lengths.
+pub fn compose_weighted(data: &[Vec<f64>], weights: &[f64]) -> Vec<f64> {
+    let n_sigs = data.len();
+    assert!(
+        weights.len() == n_sigs,
+        "weights length ({}) must equal number of signal columns ({})",
+        weights.len(),
+        n_sigs
+    );
+    if n_sigs == 0 {
+        return vec![];
+    }
+    let n_bars = data[0].len();
+    for j in 1..n_sigs {
+        assert!(
+            data[j].len() == n_bars,
+            "all columns must have equal length"
+        );
+    }
+    let mut result = vec![0.0_f64; n_bars];
+    for i in 0..n_bars {
+        let mut s = 0.0_f64;
+        for j in 0..n_sigs {
+            s += data[j][i] * weights[j];
+        }
+        result[i] = s;
+    }
+    result
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    const EPS: f64 = 1e-10;
+
+    fn approx_eq(a: f64, b: f64) -> bool {
+        (a - b).abs() < EPS
+    }
+
+    // -- portfolio_volatility -------------------------------------------------
+
+    #[test]
+    fn test_portfolio_volatility_identity_cov() {
+        // Identity covariance, equal weights => sqrt(sum(w_i^2))
+        let cov = vec![
+            vec![1.0, 0.0],
+            vec![0.0, 1.0],
+        ];
+        let w = vec![0.5, 0.5];
+        let vol = portfolio_volatility(&cov, &w);
+        // w' I w = 0.25 + 0.25 = 0.5, sqrt = 0.7071...
+        assert!(approx_eq(vol, (0.5_f64).sqrt()));
+    }
+
+    #[test]
+    fn test_portfolio_volatility_single_asset() {
+        let cov = vec![vec![0.04]];
+        let w = vec![1.0];
+        assert!(approx_eq(portfolio_volatility(&cov, &w), 0.2));
+    }
+
+    #[test]
+    fn test_portfolio_volatility_correlated() {
+        // Fully correlated: cov = [[0.04, 0.04], [0.04, 0.04]]
+        let cov = vec![
+            vec![0.04, 0.04],
+            vec![0.04, 0.04],
+        ];
+        let w = vec![0.5, 0.5];
+        // w' Σ w = 0.04, sqrt = 0.2
+        let vol = portfolio_volatility(&cov, &w);
+        assert!(approx_eq(vol, 0.2));
+    }
+
+    // -- beta_full ------------------------------------------------------------
+
+    #[test]
+    fn test_beta_full_same_series() {
+        let r = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        assert!(approx_eq(beta_full(&r, &r), 1.0));
+    }
+
+    #[test]
+    fn test_beta_full_double() {
+        let bench = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        let asset: Vec<f64> = bench.iter().map(|x| x * 2.0).collect();
+        assert!(approx_eq(beta_full(&asset, &bench), 2.0));
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_beta_full_zero_variance() {
+        let a = vec![0.01, 0.02];
+        let b = vec![0.05, 0.05]; // zero variance
+        beta_full(&a, &b);
+    }
+
+    // -- rolling_beta ---------------------------------------------------------
+
+    #[test]
+    fn test_rolling_beta_warmup_nan() {
+        let a = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        let b = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        let rb = rolling_beta(&a, &b, 3);
+        assert_eq!(rb.len(), 5);
+        assert!(rb[0].is_nan());
+        assert!(rb[1].is_nan());
+        // From index 2 onward, beta of identical series = 1.0
+        assert!(approx_eq(rb[2], 1.0));
+        assert!(approx_eq(rb[3], 1.0));
+        assert!(approx_eq(rb[4], 1.0));
+    }
+
+    #[test]
+    fn test_rolling_beta_double() {
+        let bench = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        let asset: Vec<f64> = bench.iter().map(|x| x * 3.0).collect();
+        let rb = rolling_beta(&asset, &bench, 3);
+        for i in 2..5 {
+            assert!(approx_eq(rb[i], 3.0));
+        }
+    }
+
+    // -- drawdown_series ------------------------------------------------------
+
+    #[test]
+    fn test_drawdown_series_monotonic_up() {
+        let eq = vec![100.0, 110.0, 120.0, 130.0];
+        let (dd, max_dd) = drawdown_series(&eq);
+        for &d in &dd {
+            assert!(approx_eq(d, 0.0));
+        }
+        assert!(approx_eq(max_dd, 0.0));
+    }
+
+    #[test]
+    fn test_drawdown_series_with_dip() {
+        let eq = vec![100.0, 120.0, 90.0, 110.0];
+        let (dd, max_dd) = drawdown_series(&eq);
+        assert!(approx_eq(dd[0], 0.0));
+        assert!(approx_eq(dd[1], 0.0));
+        // dd[2] = (90 - 120) / 120 = -0.25
+        assert!(approx_eq(dd[2], -0.25));
+        // dd[3] = (110 - 120) / 120 = -1/12
+        assert!((dd[3] - (-1.0 / 12.0)).abs() < EPS);
+        assert!(approx_eq(max_dd, -0.25));
+    }
+
+    // -- correlation_matrix ---------------------------------------------------
+
+    #[test]
+    fn test_correlation_matrix_identical() {
+        let col = vec![0.01, -0.02, 0.03, -0.01, 0.02];
+        let data = vec![col.clone(), col.clone()];
+        let cm = correlation_matrix(&data);
+        assert_eq!(cm.len(), 2);
+        assert!(approx_eq(cm[0][0], 1.0));
+        assert!(approx_eq(cm[1][1], 1.0));
+        assert!(approx_eq(cm[0][1], 1.0));
+        assert!(approx_eq(cm[1][0], 1.0));
+    }
+
+    #[test]
+    fn test_correlation_matrix_negatively_correlated() {
+        let col_a = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let col_b: Vec<f64> = col_a.iter().map(|x| -x).collect();
+        let data = vec![col_a, col_b];
+        let cm = correlation_matrix(&data);
+        assert!(approx_eq(cm[0][1], -1.0));
+        assert!(approx_eq(cm[1][0], -1.0));
+    }
+
+    #[test]
+    fn test_correlation_matrix_single_asset() {
+        let data = vec![vec![1.0, 2.0, 3.0]];
+        let cm = correlation_matrix(&data);
+        assert_eq!(cm.len(), 1);
+        assert!(approx_eq(cm[0][0], 1.0));
+    }
+
+    // -- relative_strength ----------------------------------------------------
+
+    #[test]
+    fn test_relative_strength_equal() {
+        let r = vec![0.01, -0.02, 0.03];
+        let rs = relative_strength(&r, &r);
+        for &v in &rs {
+            assert!(approx_eq(v, 1.0));
+        }
+    }
+
+    #[test]
+    fn test_relative_strength_outperformance() {
+        let a = vec![0.10, 0.10];
+        let b = vec![0.05, 0.05];
+        let rs = relative_strength(&a, &b);
+        // rs[0] = 1.10 / 1.05
+        assert!((rs[0] - 1.10 / 1.05).abs() < EPS);
+        // rs[1] = 1.21 / 1.1025
+        assert!((rs[1] - 1.21 / 1.1025).abs() < EPS);
+    }
+
+    // -- spread ---------------------------------------------------------------
+
+    #[test]
+    fn test_spread_basic() {
+        let a = vec![10.0, 20.0, 30.0];
+        let b = vec![5.0, 10.0, 15.0];
+        let s = spread(&a, &b, 2.0);
+        assert!(approx_eq(s[0], 0.0));
+        assert!(approx_eq(s[1], 0.0));
+        assert!(approx_eq(s[2], 0.0));
+    }
+
+    #[test]
+    fn test_spread_hedge_one() {
+        let a = vec![10.0, 20.0];
+        let b = vec![3.0, 7.0];
+        let s = spread(&a, &b, 1.0);
+        assert!(approx_eq(s[0], 7.0));
+        assert!(approx_eq(s[1], 13.0));
+    }
+
+    // -- ratio ----------------------------------------------------------------
+
+    #[test]
+    fn test_ratio_basic() {
+        let a = vec![10.0, 20.0, 30.0];
+        let b = vec![5.0, 10.0, 15.0];
+        let r = ratio(&a, &b);
+        assert!(approx_eq(r[0], 2.0));
+        assert!(approx_eq(r[1], 2.0));
+        assert!(approx_eq(r[2], 2.0));
+    }
+
+    #[test]
+    fn test_ratio_zero_denominator() {
+        let a = vec![10.0, 20.0];
+        let b = vec![0.0, 5.0];
+        let r = ratio(&a, &b);
+        assert!(r[0].is_nan());
+        assert!(approx_eq(r[1], 4.0));
+    }
+
+    // -- zscore_series --------------------------------------------------------
+
+    #[test]
+    fn test_zscore_warmup_nan() {
+        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let z = zscore_series(&x, 3);
+        assert!(z[0].is_nan());
+        assert!(z[1].is_nan());
+        assert!(!z[2].is_nan());
+        assert!(!z[3].is_nan());
+        assert!(!z[4].is_nan());
+    }
+
+    #[test]
+    fn test_zscore_constant_window() {
+        // All same values in window => std = 0 => NaN
+        let x = vec![5.0, 5.0, 5.0, 5.0];
+        let z = zscore_series(&x, 3);
+        assert!(z[2].is_nan());
+        assert!(z[3].is_nan());
+    }
+
+    #[test]
+    fn test_zscore_known_value() {
+        // Window [1, 2, 3]: mean=2, pop_std = sqrt(2/3) ~0.8165
+        // z = (3 - 2) / sqrt(2/3) = sqrt(3/2) ~ 1.2247
+        let x = vec![1.0, 2.0, 3.0];
+        let z = zscore_series(&x, 3);
+        let expected = (3.0_f64 / 2.0).sqrt();
+        assert!((z[2] - expected).abs() < EPS);
+    }
+
+    // -- compose_weighted -----------------------------------------------------
+
+    #[test]
+    fn test_compose_weighted_basic() {
+        let data = vec![
+            vec![1.0, 2.0, 3.0],
+            vec![4.0, 5.0, 6.0],
+        ];
+        let weights = vec![0.3, 0.7];
+        let cw = compose_weighted(&data, &weights);
+        // bar 0: 1*0.3 + 4*0.7 = 3.1
+        assert!(approx_eq(cw[0], 3.1));
+        // bar 1: 2*0.3 + 5*0.7 = 4.1
+        assert!(approx_eq(cw[1], 4.1));
+        // bar 2: 3*0.3 + 6*0.7 = 5.1
+        assert!(approx_eq(cw[2], 5.1));
+    }
+
+    #[test]
+    fn test_compose_weighted_single_column() {
+        let data = vec![vec![10.0, 20.0]];
+        let weights = vec![2.0];
+        let cw = compose_weighted(&data, &weights);
+        assert!(approx_eq(cw[0], 20.0));
+        assert!(approx_eq(cw[1], 40.0));
+    }
+
+    #[test]
+    fn test_compose_weighted_empty() {
+        let data: Vec<Vec<f64>> = vec![];
+        let weights: Vec<f64> = vec![];
+        let cw = compose_weighted(&data, &weights);
+        assert!(cw.is_empty());
+    }
+}

crates/ferro_ta_core/src/price_transform.rs

@@ -0,0 +1,89 @@
+//! Price transformations — synthesize OHLC arrays into single price arrays.
+
+/// Average Price: (open + high + low + close) / 4.
+pub fn avgprice(open: &[f64], high: &[f64], low: &[f64], close: &[f64]) -> Vec<f64> {
+    open.iter()
+        .zip(high.iter())
+        .zip(low.iter())
+        .zip(close.iter())
+        .map(|(((&o, &h), &l), &c)| (o + h + l + c) / 4.0)
+        .collect()
+}
+
+/// Median Price: (high + low) / 2.
+pub fn medprice(high: &[f64], low: &[f64]) -> Vec<f64> {
+    high.iter()
+        .zip(low.iter())
+        .map(|(&h, &l)| (h + l) / 2.0)
+        .collect()
+}
+
+/// Typical Price: (high + low + close) / 3.
+pub fn typprice(high: &[f64], low: &[f64], close: &[f64]) -> Vec<f64> {
+    high.iter()
+        .zip(low.iter())
+        .zip(close.iter())
+        .map(|((&h, &l), &c)| (h + l + c) / 3.0)
+        .collect()
+}
+
+/// Weighted Close Price: (high + low + close * 2) / 4.
+pub fn wclprice(high: &[f64], low: &[f64], close: &[f64]) -> Vec<f64> {
+    high.iter()
+        .zip(low.iter())
+        .zip(close.iter())
+        .map(|((&h, &l), &c)| (h + l + c * 2.0) / 4.0)
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_avgprice() {
+        let o = vec![1.0, 2.0, 3.0];
+        let h = vec![4.0, 5.0, 6.0];
+        let l = vec![0.5, 1.5, 2.5];
+        let c = vec![2.5, 3.5, 4.5];
+        let result = avgprice(&o, &h, &l, &c);
+        assert_eq!(result.len(), 3);
+        assert!((result[0] - 2.0).abs() < 1e-10); // (1+4+0.5+2.5)/4 = 2.0
+    }
+
+    #[test]
+    fn test_medprice() {
+        let h = vec![10.0, 20.0];
+        let l = vec![6.0, 12.0];
+        let result = medprice(&h, &l);
+        assert!((result[0] - 8.0).abs() < 1e-10);
+        assert!((result[1] - 16.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_typprice() {
+        let h = vec![10.0];
+        let l = vec![6.0];
+        let c = vec![8.0];
+        let result = typprice(&h, &l, &c);
+        assert!((result[0] - 8.0).abs() < 1e-10); // (10+6+8)/3 = 8.0
+    }
+
+    #[test]
+    fn test_wclprice() {
+        let h = vec![10.0];
+        let l = vec![6.0];
+        let c = vec![8.0];
+        let result = wclprice(&h, &l, &c);
+        assert!((result[0] - 8.0).abs() < 1e-10); // (10+6+16)/4 = 8.0
+    }
+
+    #[test]
+    fn test_empty_inputs() {
+        let empty: Vec<f64> = vec![];
+        assert!(avgprice(&empty, &empty, &empty, &empty).is_empty());
+        assert!(medprice(&empty, &empty).is_empty());
+        assert!(typprice(&empty, &empty, &empty).is_empty());
+        assert!(wclprice(&empty, &empty, &empty).is_empty());
+    }
+}

crates/ferro_ta_core/src/regime.rs

@@ -0,0 +1,171 @@
+//! Regime detection and structural breaks.
+//!
+//! - `regime_adx`             — label trend (1) vs range (0) using ADX threshold
+//! - `regime_combined`        — combine ADX + ATR-ratio for robust regime labelling
+//! - `detect_breaks_cusum`    — CUSUM-based structural break detection
+//! - `rolling_variance_break` — variance ratio break detection
+
+/// Label each bar as trend (1) or range (0) based on ADX level.
+///
+/// Returns `Vec<i8>`: `1` = trend (ADX > threshold), `0` = range, `-1` = NaN/warmup.
+pub fn regime_adx(adx: &[f64], threshold: f64) -> Vec<i8> {
+    adx.iter()
+        .map(|&v| {
+            if v.is_nan() {
+                -1i8
+            } else if v > threshold {
+                1i8
+            } else {
+                0i8
+            }
+        })
+        .collect()
+}
+
+/// Label each bar as trend (1) or range (0) using ADX + ATR-ratio rule.
+///
+/// A bar is trending when: `adx[i] > adx_threshold` AND `atr[i] / close[i] > atr_pct_threshold`.
+///
+/// Returns `Vec<i8>`: `1` = trend, `0` = range, `-1` = NaN.
+pub fn regime_combined(
+    adx: &[f64],
+    atr: &[f64],
+    close: &[f64],
+    adx_threshold: f64,
+    atr_pct_threshold: f64,
+) -> Vec<i8> {
+    let n = adx.len();
+    (0..n)
+        .map(|i| {
+            let av = adx[i];
+            let rv = atr[i];
+            let cv = close[i];
+            if av.is_nan() || rv.is_nan() || cv.is_nan() || cv == 0.0 {
+                -1i8
+            } else if av > adx_threshold && (rv / cv) > atr_pct_threshold {
+                1i8
+            } else {
+                0i8
+            }
+        })
+        .collect()
+}
+
+/// Detect structural breaks using a CUSUM (cumulative sum) approach.
+///
+/// `window` must be >= 2. Returns `Vec<i8>`: `1` at break bars, `0` elsewhere.
+pub fn detect_breaks_cusum(
+    series: &[f64],
+    window: usize,
+    threshold: f64,
+    slack: f64,
+) -> Vec<i8> {
+    let n = series.len();
+    let mut out = vec![0i8; n];
+    if n < window || window < 2 {
+        return out;
+    }
+    let mut cusum_pos = 0.0_f64;
+    let mut cusum_neg = 0.0_f64;
+    for i in window..n {
+        let slice = &series[(i - window)..i];
+        let mean: f64 = slice.iter().sum::<f64>() / window as f64;
+        let var: f64 =
+            slice.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>() / (window - 1) as f64;
+        let std = var.sqrt();
+        if std == 0.0 || std.is_nan() || series[i].is_nan() {
+            continue;
+        }
+        let z = (series[i] - mean) / std;
+        cusum_pos = (cusum_pos + z - slack).max(0.0);
+        cusum_neg = (cusum_neg - z - slack).max(0.0);
+        if cusum_pos > threshold || cusum_neg > threshold {
+            out[i] = 1;
+            cusum_pos = 0.0;
+            cusum_neg = 0.0;
+        }
+    }
+    out
+}
+
+/// Detect volatility regime breaks using rolling variance ratio.
+///
+/// `short_window` must be >= 2, `long_window` must be > `short_window`.
+/// Returns `Vec<i8>`: `1` at break bars, `0` elsewhere.
+pub fn rolling_variance_break(
+    series: &[f64],
+    short_window: usize,
+    long_window: usize,
+    threshold: f64,
+) -> Vec<i8> {
+    let n = series.len();
+    let mut out = vec![0i8; n];
+    if n < long_window || short_window < 2 || long_window <= short_window {
+        return out;
+    }
+
+    let variance = |slice: &[f64]| -> f64 {
+        let k = slice.len();
+        let mean: f64 = slice.iter().sum::<f64>() / k as f64;
+        slice.iter().map(|&v| (v - mean) * (v - mean)).sum::<f64>() / (k - 1) as f64
+    };
+
+    for i in long_window..n {
+        let long_slice = &series[(i - long_window)..i];
+        let short_slice = &series[(i - short_window)..i];
+        let long_var = variance(long_slice);
+        let short_var = variance(short_slice);
+        if long_var == 0.0 || long_var.is_nan() || short_var.is_nan() {
+            continue;
+        }
+        if short_var / long_var > threshold {
+            out[i] = 1;
+        }
+    }
+    out
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_regime_adx_basic() {
+        let adx = vec![f64::NAN, 20.0, 30.0, 10.0, 50.0];
+        let result = regime_adx(&adx, 25.0);
+        assert_eq!(result, vec![-1, 0, 1, 0, 1]);
+    }
+
+    #[test]
+    fn test_regime_combined() {
+        let adx = vec![30.0, 30.0, 10.0];
+        let atr = vec![1.0, 0.001, 1.0];
+        let close = vec![100.0, 100.0, 100.0];
+        let result = regime_combined(&adx, &atr, &close, 25.0, 0.005);
+        assert_eq!(result[0], 1); // ADX>25 and ATR/close=0.01>0.005
+        assert_eq!(result[1], 0); // ATR/close=0.00001 < 0.005
+        assert_eq!(result[2], 0); // ADX<25
+    }
+
+    #[test]
+    fn test_detect_breaks_cusum_short_input() {
+        let series = vec![1.0, 2.0];
+        let result = detect_breaks_cusum(&series, 5, 3.0, 0.5);
+        assert!(result.iter().all(|&v| v == 0));
+    }
+
+    #[test]
+    fn test_rolling_variance_break_short_input() {
+        let series = vec![1.0, 2.0, 3.0];
+        let result = rolling_variance_break(&series, 2, 5, 2.0);
+        assert!(result.iter().all(|&v| v == 0));
+    }
+
+    #[test]
+    fn test_empty() {
+        assert!(regime_adx(&[], 25.0).is_empty());
+        assert!(regime_combined(&[], &[], &[], 25.0, 0.005).is_empty());
+        assert!(detect_breaks_cusum(&[], 2, 3.0, 0.5).is_empty());
+        assert!(rolling_variance_break(&[], 2, 5, 2.0).is_empty());
+    }
+}

crates/ferro_ta_core/src/resampling.rs

@@ -0,0 +1,278 @@
+//! Resampling — OHLCV resampling and multi-timeframe helpers, pure Rust.
+//!
+//! # Functions
+//! - `volume_bars` — Aggregate OHLCV bars into bars of fixed volume size.
+//! - `ohlcv_agg`   — Aggregate OHLCV bars given contiguous integer group labels.
+
+/// OHLCV 5-tuple return type alias.
+type Ohlcv5 = (Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>, Vec<f64>);
+
+// ---------------------------------------------------------------------------
+// volume_bars
+// ---------------------------------------------------------------------------
+
+/// Aggregate OHLCV data into volume bars of a fixed volume threshold.
+///
+/// Each output bar accumulates input bars until `volume_threshold` units of
+/// volume have been consumed.  The resulting bar has:
+///   - open  = first open of the group
+///   - high  = max high of the group
+///   - low   = min low of the group
+///   - close = last close of the group
+///   - volume = sum of volumes (approximately `volume_threshold`)
+///
+/// Returns `(open, high, low, close, volume)`.
+///
+/// # Panics
+/// Panics if arrays are empty, have unequal lengths, or `volume_threshold <= 0`.
+pub fn volume_bars(
+    open: &[f64],
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    volume: &[f64],
+    volume_threshold: f64,
+) -> Ohlcv5 {
+    assert!(volume_threshold > 0.0, "volume_threshold must be > 0");
+    let n = open.len();
+    assert!(n > 0, "input arrays must be non-empty");
+    assert!(
+        high.len() == n && low.len() == n && close.len() == n && volume.len() == n,
+        "all input arrays must have equal length"
+    );
+
+    let mut out_open: Vec<f64> = Vec::new();
+    let mut out_high: Vec<f64> = Vec::new();
+    let mut out_low: Vec<f64> = Vec::new();
+    let mut out_close: Vec<f64> = Vec::new();
+    let mut out_vol: Vec<f64> = Vec::new();
+
+    let mut bar_open = open[0];
+    let mut bar_high = high[0];
+    let mut bar_low = low[0];
+    let mut bar_close = close[0];
+    let mut bar_vol = volume[0];
+
+    for i in 1..n {
+        bar_high = bar_high.max(high[i]);
+        bar_low = bar_low.min(low[i]);
+        bar_close = close[i];
+        bar_vol += volume[i];
+
+        if bar_vol >= volume_threshold {
+            out_open.push(bar_open);
+            out_high.push(bar_high);
+            out_low.push(bar_low);
+            out_close.push(bar_close);
+            out_vol.push(bar_vol);
+            // Start new bar
+            if i + 1 < n {
+                bar_open = open[i + 1];
+                bar_high = high[i + 1];
+                bar_low = low[i + 1];
+                bar_close = close[i + 1];
+                bar_vol = volume[i + 1];
+            }
+        }
+    }
+    // Push any remaining partial bar
+    if bar_vol > 0.0 && out_vol.last().is_none_or(|&last| last != bar_vol) {
+        out_open.push(bar_open);
+        out_high.push(bar_high);
+        out_low.push(bar_low);
+        out_close.push(bar_close);
+        out_vol.push(bar_vol);
+    }
+
+    (out_open, out_high, out_low, out_close, out_vol)
+}
+
+// ---------------------------------------------------------------------------
+// ohlcv_agg
+// ---------------------------------------------------------------------------
+
+/// Aggregate OHLCV bars by integer group labels.
+///
+/// Groups consecutive bars with the same label and computes:
+///   - open  = first open of the group
+///   - high  = max high of the group
+///   - low   = min low of the group
+///   - close = last close of the group
+///   - volume = sum of volumes
+///
+/// `labels` must be non-decreasing (groups are contiguous).
+///
+/// Returns `(open, high, low, close, volume)`.
+///
+/// # Panics
+/// Panics if arrays are empty or have unequal lengths.
+pub fn ohlcv_agg(
+    open: &[f64],
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    volume: &[f64],
+    labels: &[i64],
+) -> Ohlcv5 {
+    let n = open.len();
+    assert!(n > 0, "input arrays must be non-empty");
+    assert!(
+        high.len() == n
+            && low.len() == n
+            && close.len() == n
+            && volume.len() == n
+            && labels.len() == n,
+        "all input arrays must have equal length"
+    );
+
+    let mut out_open: Vec<f64> = Vec::new();
+    let mut out_high: Vec<f64> = Vec::new();
+    let mut out_low: Vec<f64> = Vec::new();
+    let mut out_close: Vec<f64> = Vec::new();
+    let mut out_vol: Vec<f64> = Vec::new();
+
+    let mut cur_label = labels[0];
+    let mut bar_open = open[0];
+    let mut bar_high = high[0];
+    let mut bar_low = low[0];
+    let mut bar_close = close[0];
+    let mut bar_vol = volume[0];
+
+    for i in 1..n {
+        if labels[i] != cur_label {
+            out_open.push(bar_open);
+            out_high.push(bar_high);
+            out_low.push(bar_low);
+            out_close.push(bar_close);
+            out_vol.push(bar_vol);
+            cur_label = labels[i];
+            bar_open = open[i];
+            bar_high = high[i];
+            bar_low = low[i];
+            bar_close = close[i];
+            bar_vol = volume[i];
+        } else {
+            bar_high = bar_high.max(high[i]);
+            bar_low = bar_low.min(low[i]);
+            bar_close = close[i];
+            bar_vol += volume[i];
+        }
+    }
+    out_open.push(bar_open);
+    out_high.push(bar_high);
+    out_low.push(bar_low);
+    out_close.push(bar_close);
+    out_vol.push(bar_vol);
+
+    (out_open, out_high, out_low, out_close, out_vol)
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // -- volume_bars ---------------------------------------------------------
+
+    #[test]
+    fn test_volume_bars_basic() {
+        let o = [100.0, 101.0, 102.0, 103.0, 104.0];
+        let h = [105.0, 106.0, 107.0, 108.0, 109.0];
+        let l = [95.0, 96.0, 97.0, 98.0, 99.0];
+        let c = [101.0, 102.0, 103.0, 104.0, 105.0];
+        let v = [50.0, 60.0, 40.0, 70.0, 30.0];
+        // threshold 100: first bar covers indices 0..2 (vol=110>=100)
+        let (ro, rh, rl, rc, rv) = volume_bars(&o, &h, &l, &c, &v, 100.0);
+        assert!(rv.len() >= 2);
+        // First bar: vol = 50+60 = 110
+        assert!((rv[0] - 110.0).abs() < 1e-10);
+        assert!((ro[0] - 100.0).abs() < 1e-10);
+        assert!((rh[0] - 106.0).abs() < 1e-10);
+        assert!((rl[0] - 95.0).abs() < 1e-10);
+        assert!((rc[0] - 102.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_volume_bars_single_element() {
+        let (ro, rh, rl, rc, rv) =
+            volume_bars(&[10.0], &[12.0], &[8.0], &[11.0], &[50.0], 100.0);
+        assert_eq!(rv.len(), 1);
+        assert!((rv[0] - 50.0).abs() < 1e-10);
+        assert!((ro[0] - 10.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "volume_threshold must be > 0")]
+    fn test_volume_bars_zero_threshold() {
+        volume_bars(&[1.0], &[1.0], &[1.0], &[1.0], &[1.0], 0.0);
+    }
+
+    #[test]
+    #[should_panic(expected = "input arrays must be non-empty")]
+    fn test_volume_bars_empty() {
+        volume_bars(&[], &[], &[], &[], &[], 100.0);
+    }
+
+    // -- ohlcv_agg -----------------------------------------------------------
+
+    #[test]
+    fn test_ohlcv_agg_basic() {
+        let o = [100.0, 101.0, 102.0, 103.0];
+        let h = [105.0, 106.0, 108.0, 109.0];
+        let l = [95.0, 96.0, 97.0, 98.0];
+        let c = [101.0, 102.0, 103.0, 104.0];
+        let v = [10.0, 20.0, 30.0, 40.0];
+        let labels: [i64; 4] = [0, 0, 1, 1];
+        let (ro, rh, rl, rc, rv) = ohlcv_agg(&o, &h, &l, &c, &v, &labels);
+        assert_eq!(ro.len(), 2);
+        // Group 0: open=100, high=max(105,106)=106, low=min(95,96)=95, close=102, vol=30
+        assert!((ro[0] - 100.0).abs() < 1e-10);
+        assert!((rh[0] - 106.0).abs() < 1e-10);
+        assert!((rl[0] - 95.0).abs() < 1e-10);
+        assert!((rc[0] - 102.0).abs() < 1e-10);
+        assert!((rv[0] - 30.0).abs() < 1e-10);
+        // Group 1: open=102, high=max(108,109)=109, low=min(97,98)=97, close=104, vol=70
+        assert!((ro[1] - 102.0).abs() < 1e-10);
+        assert!((rh[1] - 109.0).abs() < 1e-10);
+        assert!((rl[1] - 97.0).abs() < 1e-10);
+        assert!((rc[1] - 104.0).abs() < 1e-10);
+        assert!((rv[1] - 70.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_ohlcv_agg_single_group() {
+        let o = [100.0, 101.0];
+        let h = [105.0, 106.0];
+        let l = [95.0, 96.0];
+        let c = [101.0, 102.0];
+        let v = [10.0, 20.0];
+        let labels: [i64; 2] = [0, 0];
+        let (ro, rh, rl, rc, rv) = ohlcv_agg(&o, &h, &l, &c, &v, &labels);
+        assert_eq!(ro.len(), 1);
+        assert!((rv[0] - 30.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_ohlcv_agg_each_bar_own_group() {
+        let o = [100.0, 101.0, 102.0];
+        let h = [105.0, 106.0, 107.0];
+        let l = [95.0, 96.0, 97.0];
+        let c = [101.0, 102.0, 103.0];
+        let v = [10.0, 20.0, 30.0];
+        let labels: [i64; 3] = [0, 1, 2];
+        let (ro, _rh, _rl, _rc, rv) = ohlcv_agg(&o, &h, &l, &c, &v, &labels);
+        assert_eq!(ro.len(), 3);
+        assert!((rv[0] - 10.0).abs() < 1e-10);
+        assert!((rv[1] - 20.0).abs() < 1e-10);
+        assert!((rv[2] - 30.0).abs() < 1e-10);
+    }
+
+    #[test]
+    #[should_panic(expected = "input arrays must be non-empty")]
+    fn test_ohlcv_agg_empty() {
+        ohlcv_agg(&[], &[], &[], &[], &[], &[]);
+    }
+}

crates/ferro_ta_core/src/signals.rs

@@ -0,0 +1,131 @@
+//! Signal processing helpers.
+//!
+//! - `rank_values`      — fractional rank of a slice (1-based, ties averaged)
+//! - `compose_rank`     — rank-based composite scores for a 2-D signal matrix
+//! - `top_n_indices`    — indices of the N largest values
+//! - `bottom_n_indices` — indices of the N smallest values
+
+/// Compute fractional rank of each element (1-based, ascending).
+/// Ties receive the average of their rank positions.
+pub fn rank_values(x: &[f64]) -> Vec<f64> {
+    let n = x.len();
+    let mut order: Vec<usize> = (0..n).collect();
+    order.sort_by(|&a, &b| x[a].partial_cmp(&x[b]).unwrap_or(std::cmp::Ordering::Equal));
+
+    let mut ranks = vec![0.0_f64; n];
+    let mut i = 0;
+    while i < n {
+        let val = x[order[i]];
+        let mut j = i + 1;
+        while j < n && x[order[j]] == val {
+            j += 1;
+        }
+        let avg_rank = (i + 1 + j) as f64 / 2.0;
+        for k in i..j {
+            ranks[order[k]] = avg_rank;
+        }
+        i = j;
+    }
+    ranks
+}
+
+/// Compute rank-based composite scores for a 2-D signal matrix.
+///
+/// Each column is ranked independently, and the per-row ranks are summed.
+/// `signals` is a slice of columns, each column being a `&[f64]` of the same length.
+pub fn compose_rank(signals: &[&[f64]]) -> Vec<f64> {
+    if signals.is_empty() {
+        return vec![];
+    }
+    let n_bars = signals[0].len();
+    let mut scores = vec![0.0_f64; n_bars];
+    for &column in signals {
+        let ranks = rank_values(column);
+        for (bar_idx, rank) in ranks.into_iter().enumerate() {
+            scores[bar_idx] += rank;
+        }
+    }
+    scores
+}
+
+/// Return the indices of the N largest values in `x` (descending by value).
+pub fn top_n_indices(x: &[f64], n: usize) -> Vec<i64> {
+    let len = x.len();
+    let k = n.min(len);
+    let mut order: Vec<usize> = (0..len).collect();
+    order.sort_by(|&a, &b| x[b].partial_cmp(&x[a]).unwrap_or(std::cmp::Ordering::Equal));
+    order[..k].iter().map(|&i| i as i64).collect()
+}
+
+/// Return the indices of the N smallest values in `x` (ascending by value).
+pub fn bottom_n_indices(x: &[f64], n: usize) -> Vec<i64> {
+    let len = x.len();
+    let k = n.min(len);
+    let mut order: Vec<usize> = (0..len).collect();
+    order.sort_by(|&a, &b| x[a].partial_cmp(&x[b]).unwrap_or(std::cmp::Ordering::Equal));
+    order[..k].iter().map(|&i| i as i64).collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_rank_values() {
+        let x = vec![3.0, 1.0, 2.0];
+        let ranks = rank_values(&x);
+        assert!((ranks[0] - 3.0).abs() < 1e-10); // 3.0 is largest → rank 3
+        assert!((ranks[1] - 1.0).abs() < 1e-10); // 1.0 is smallest → rank 1
+        assert!((ranks[2] - 2.0).abs() < 1e-10); // 2.0 is middle → rank 2
+    }
+
+    #[test]
+    fn test_rank_values_ties() {
+        let x = vec![1.0, 2.0, 2.0, 4.0];
+        let ranks = rank_values(&x);
+        assert!((ranks[0] - 1.0).abs() < 1e-10);
+        assert!((ranks[1] - 2.5).abs() < 1e-10); // tied → average
+        assert!((ranks[2] - 2.5).abs() < 1e-10);
+        assert!((ranks[3] - 4.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_compose_rank() {
+        let col1 = vec![3.0, 1.0, 2.0];
+        let col2 = vec![1.0, 3.0, 2.0];
+        let signals: Vec<&[f64]> = vec![&col1, &col2];
+        let scores = compose_rank(&signals);
+        // Row 0: rank(3.0)=3 + rank(1.0)=1 = 4
+        // Row 1: rank(1.0)=1 + rank(3.0)=3 = 4
+        // Row 2: rank(2.0)=2 + rank(2.0)=2 = 4
+        assert!((scores[0] - 4.0).abs() < 1e-10);
+        assert!((scores[1] - 4.0).abs() < 1e-10);
+        assert!((scores[2] - 4.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_top_n_indices() {
+        let x = vec![10.0, 50.0, 30.0, 20.0, 40.0];
+        let result = top_n_indices(&x, 3);
+        assert_eq!(result.len(), 3);
+        assert_eq!(result[0], 1); // 50.0
+        assert_eq!(result[1], 4); // 40.0
+        assert_eq!(result[2], 2); // 30.0
+    }
+
+    #[test]
+    fn test_bottom_n_indices() {
+        let x = vec![10.0, 50.0, 30.0, 20.0, 40.0];
+        let result = bottom_n_indices(&x, 2);
+        assert_eq!(result.len(), 2);
+        assert_eq!(result[0], 0); // 10.0
+        assert_eq!(result[1], 3); // 20.0
+    }
+
+    #[test]
+    fn test_top_n_exceeds_len() {
+        let x = vec![1.0, 2.0];
+        let result = top_n_indices(&x, 5);
+        assert_eq!(result.len(), 2);
+    }
+}

crates/ferro_ta_core/src/statistic.rs

@@ -1,6 +1,14 @@
 //! Statistic functions.
 
-/// Standard deviation — population (`ddof = 0`).
+/// Compute the rolling population standard deviation, scaled by `nbdev`.
+///
+/// Uses population variance (`ddof = 0`). Returns `nbdev * stddev` for
+/// each window. The first `timeperiod - 1` values are `NaN`.
+///
+/// # Arguments
+/// * `real` - Input series.
+/// * `timeperiod` - Rolling window size (must be >= 1).
+/// * `nbdev` - Multiplier applied to the standard deviation (use 1.0 for raw stddev).
 pub fn stddev(real: &[f64], timeperiod: usize, nbdev: f64) -> Vec<f64> {
     let n = real.len();
     let mut result = vec![f64::NAN; n];

crates/ferro_ta_core/src/streaming.rs

@@ -0,0 +1,947 @@
+//! Streaming / Incremental Indicators — bar-by-bar stateful structs.
+//!
+//! Pure Rust implementations with no PyO3 dependency.  Each struct:
+//! - Accepts one value per call to `update()`.
+//! - Returns `NaN` (or a NaN tuple) during the warm-up window.
+//! - Exposes a `reset()` method to restart from scratch.
+//! - Has a `period()` accessor (where applicable).
+
+use std::collections::VecDeque;
+
+// ---------------------------------------------------------------------------
+// Error type
+// ---------------------------------------------------------------------------
+
+/// Validation error for streaming indicator parameters.
+#[derive(Debug, Clone)]
+pub struct StreamingError(pub String);
+
+impl std::fmt::Display for StreamingError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        write!(f, "{}", self.0)
+    }
+}
+
+impl std::error::Error for StreamingError {}
+
+fn validate_timeperiod(value: usize, name: &str, minimum: usize) -> Result<(), StreamingError> {
+    if value < minimum {
+        return Err(StreamingError(format!(
+            "{} must be >= {}, got {}",
+            name, minimum, value
+        )));
+    }
+    Ok(())
+}
+
+// ---------------------------------------------------------------------------
+// Internal helper: EMA state (used inside composite classes)
+// ---------------------------------------------------------------------------
+
+/// SMA-seeded EMA state machine.  Not exposed directly — used by
+/// `StreamingEMA`, `StreamingMACD`, etc.
+pub(crate) struct EmaState {
+    period: usize,
+    alpha: f64,
+    ema: f64,
+    seed_buf: Vec<f64>,
+    seeded: bool,
+}
+
+impl EmaState {
+    pub fn new(period: usize) -> Self {
+        Self {
+            period,
+            alpha: 2.0 / (period as f64 + 1.0),
+            ema: 0.0,
+            seed_buf: Vec::with_capacity(period),
+            seeded: false,
+        }
+    }
+
+    pub fn update(&mut self, value: f64) -> f64 {
+        if !self.seeded {
+            self.seed_buf.push(value);
+            if self.seed_buf.len() < self.period {
+                return f64::NAN;
+            }
+            let seed = self.seed_buf.iter().sum::<f64>() / self.period as f64;
+            self.ema = seed;
+            self.seeded = true;
+            return seed;
+        }
+        self.ema += self.alpha * (value - self.ema);
+        self.ema
+    }
+
+    pub fn reset(&mut self) {
+        self.ema = 0.0;
+        self.seed_buf.clear();
+        self.seeded = false;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Internal helper: ATR state (Wilder smoothing)
+// ---------------------------------------------------------------------------
+
+/// Wilder-smoothed ATR state machine.  Used by `StreamingATR` and
+/// `StreamingSupertrend`.
+pub(crate) struct AtrState {
+    period: usize,
+    prev_close: f64,
+    tr_buf: Vec<f64>,
+    atr: f64,
+    seeded: bool,
+    has_prev: bool,
+}
+
+impl AtrState {
+    pub fn new(period: usize) -> Self {
+        Self {
+            period,
+            prev_close: 0.0,
+            tr_buf: Vec::with_capacity(period),
+            atr: 0.0,
+            seeded: false,
+            has_prev: false,
+        }
+    }
+
+    pub fn update(&mut self, high: f64, low: f64, close: f64) -> f64 {
+        let tr = if self.has_prev {
+            let hl = high - low;
+            let hc = (high - self.prev_close).abs();
+            let lc = (low - self.prev_close).abs();
+            hl.max(hc).max(lc)
+        } else {
+            high - low
+        };
+        self.prev_close = close;
+        self.has_prev = true;
+
+        if !self.seeded {
+            self.tr_buf.push(tr);
+            if self.tr_buf.len() < self.period {
+                return f64::NAN;
+            }
+            let seed = self.tr_buf.iter().sum::<f64>() / self.period as f64;
+            self.atr = seed;
+            self.seeded = true;
+            return f64::NAN; // first `period` bars (including this one) return NaN
+        }
+        let pf = (self.period - 1) as f64;
+        self.atr = (self.atr * pf + tr) / self.period as f64;
+        self.atr
+    }
+
+    pub fn reset(&mut self) {
+        self.prev_close = 0.0;
+        self.has_prev = false;
+        self.tr_buf.clear();
+        self.atr = 0.0;
+        self.seeded = false;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingSMA
+// ---------------------------------------------------------------------------
+
+/// Simple Moving Average — O(1) per update via running sum.
+///
+/// Returns NaN during the first `period - 1` bars.
+pub struct StreamingSMA {
+    period: usize,
+    buf: VecDeque<f64>,
+    running_sum: f64,
+    count: usize,
+}
+
+impl StreamingSMA {
+    pub fn new(period: usize) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 1)?;
+        Ok(Self {
+            period,
+            buf: VecDeque::with_capacity(period + 1),
+            running_sum: 0.0,
+            count: 0,
+        })
+    }
+
+    /// Add a new bar and return the current SMA (NaN during warmup).
+    pub fn update(&mut self, value: f64) -> f64 {
+        if self.buf.len() == self.period {
+            if let Some(old) = self.buf.pop_front() {
+                self.running_sum -= old;
+            }
+        }
+        self.buf.push_back(value);
+        self.running_sum += value;
+        self.count += 1;
+        if self.count < self.period {
+            f64::NAN
+        } else {
+            self.running_sum / self.period as f64
+        }
+    }
+
+    /// Reset state to initial condition.
+    pub fn reset(&mut self) {
+        self.buf.clear();
+        self.running_sum = 0.0;
+        self.count = 0;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingEMA
+// ---------------------------------------------------------------------------
+
+/// Exponential Moving Average with SMA seeding.
+///
+/// Uses a simple SMA for the first `period` bars to seed the EMA, then
+/// switches to the standard EMA formula (alpha = 2 / (period + 1)).
+/// Returns NaN during the warmup window.
+pub struct StreamingEMA {
+    inner: EmaState,
+}
+
+impl StreamingEMA {
+    pub fn new(period: usize) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 1)?;
+        Ok(Self {
+            inner: EmaState::new(period),
+        })
+    }
+
+    /// Add a new bar and return the current EMA (NaN during warmup).
+    pub fn update(&mut self, value: f64) -> f64 {
+        self.inner.update(value)
+    }
+
+    pub fn reset(&mut self) {
+        self.inner.reset();
+    }
+
+    pub fn period(&self) -> usize {
+        self.inner.period()
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingRSI
+// ---------------------------------------------------------------------------
+
+/// Relative Strength Index with TA-Lib-compatible Wilder seeding.
+///
+/// Returns NaN during the first `period` bars.
+pub struct StreamingRSI {
+    period: usize,
+    prev: f64,
+    has_prev: bool,
+    gains: Vec<f64>,
+    losses: Vec<f64>,
+    avg_gain: f64,
+    avg_loss: f64,
+    seeded: bool,
+}
+
+impl StreamingRSI {
+    pub fn new(period: usize) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 1)?;
+        Ok(Self {
+            period,
+            prev: 0.0,
+            has_prev: false,
+            gains: Vec::with_capacity(period),
+            losses: Vec::with_capacity(period),
+            avg_gain: 0.0,
+            avg_loss: 0.0,
+            seeded: false,
+        })
+    }
+
+    /// Add a new close and return RSI in [0, 100] (NaN during warmup).
+    pub fn update(&mut self, value: f64) -> f64 {
+        if !self.has_prev {
+            self.prev = value;
+            self.has_prev = true;
+            return f64::NAN;
+        }
+        let delta = value - self.prev;
+        self.prev = value;
+        let gain = if delta > 0.0 { delta } else { 0.0 };
+        let loss = if delta < 0.0 { -delta } else { 0.0 };
+
+        if !self.seeded {
+            self.gains.push(gain);
+            self.losses.push(loss);
+            if self.gains.len() < self.period {
+                return f64::NAN;
+            }
+            self.avg_gain = self.gains.iter().sum::<f64>() / self.period as f64;
+            self.avg_loss = self.losses.iter().sum::<f64>() / self.period as f64;
+            self.seeded = true;
+        } else {
+            let pf = (self.period - 1) as f64;
+            self.avg_gain = (self.avg_gain * pf + gain) / self.period as f64;
+            self.avg_loss = (self.avg_loss * pf + loss) / self.period as f64;
+        }
+
+        if self.avg_loss == 0.0 {
+            return 100.0;
+        }
+        let rs = self.avg_gain / self.avg_loss;
+        100.0 - 100.0 / (1.0 + rs)
+    }
+
+    pub fn reset(&mut self) {
+        self.prev = 0.0;
+        self.has_prev = false;
+        self.gains.clear();
+        self.losses.clear();
+        self.avg_gain = 0.0;
+        self.avg_loss = 0.0;
+        self.seeded = false;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingATR
+// ---------------------------------------------------------------------------
+
+/// Average True Range with TA-Lib-compatible Wilder seeding.
+///
+/// Accepts (high, low, close) per bar.
+/// Returns NaN during the first `period` bars.
+pub struct StreamingATR {
+    inner: AtrState,
+}
+
+impl StreamingATR {
+    pub fn new(period: usize) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 1)?;
+        Ok(Self {
+            inner: AtrState::new(period),
+        })
+    }
+
+    /// Add a new bar (high, low, close) and return ATR (NaN during warmup).
+    pub fn update(&mut self, high: f64, low: f64, close: f64) -> f64 {
+        self.inner.update(high, low, close)
+    }
+
+    pub fn reset(&mut self) {
+        self.inner.reset();
+    }
+
+    pub fn period(&self) -> usize {
+        self.inner.period()
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingBBands
+// ---------------------------------------------------------------------------
+
+/// Bollinger Bands — streaming variant using Welford's online algorithm.
+///
+/// Returns (upper, middle, lower).
+/// NaN tuple during the warmup window.
+pub struct StreamingBBands {
+    period: usize,
+    nbdevup: f64,
+    nbdevdn: f64,
+    buf: VecDeque<f64>,
+    mean: f64,
+    m2: f64,
+}
+
+impl StreamingBBands {
+    pub fn new(period: usize, nbdevup: f64, nbdevdn: f64) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 2)?;
+        Ok(Self {
+            period,
+            nbdevup,
+            nbdevdn,
+            buf: VecDeque::with_capacity(period + 1),
+            mean: 0.0,
+            m2: 0.0,
+        })
+    }
+
+    /// Add a new bar; return (upper, middle, lower). NaN tuple during warmup.
+    pub fn update(&mut self, value: f64) -> (f64, f64, f64) {
+        let n = self.buf.len();
+
+        if n == self.period {
+            let x_old = self.buf.pop_front().unwrap();
+            let count = self.period as f64;
+            let delta_old = x_old - self.mean;
+            self.mean -= delta_old / (count - 1.0);
+            let delta2_old = x_old - self.mean;
+            self.m2 -= delta_old * delta2_old;
+        }
+
+        self.buf.push_back(value);
+        let count = self.buf.len() as f64;
+        let delta_new = value - self.mean;
+        self.mean += delta_new / count;
+        let delta2_new = value - self.mean;
+        self.m2 += delta_new * delta2_new;
+
+        if self.m2 < 0.0 {
+            self.m2 = 0.0;
+        }
+
+        if self.buf.len() < self.period {
+            return (f64::NAN, f64::NAN, f64::NAN);
+        }
+
+        let variance = self.m2 / (count - 1.0);
+        let std = variance.sqrt();
+        (
+            self.mean + self.nbdevup * std,
+            self.mean,
+            self.mean - self.nbdevdn * std,
+        )
+    }
+
+    pub fn reset(&mut self) {
+        self.buf.clear();
+        self.mean = 0.0;
+        self.m2 = 0.0;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingMACD
+// ---------------------------------------------------------------------------
+
+/// MACD — fast EMA, slow EMA, signal EMA.
+///
+/// Returns (macd_line, signal_line, histogram).
+/// NaN values during warmup.
+pub struct StreamingMACD {
+    fast: EmaState,
+    slow: EmaState,
+    signal: EmaState,
+}
+
+impl StreamingMACD {
+    pub fn new(
+        fastperiod: usize,
+        slowperiod: usize,
+        signalperiod: usize,
+    ) -> Result<Self, StreamingError> {
+        validate_timeperiod(fastperiod, "fastperiod", 1)?;
+        validate_timeperiod(slowperiod, "slowperiod", 1)?;
+        validate_timeperiod(signalperiod, "signalperiod", 1)?;
+        if fastperiod >= slowperiod {
+            return Err(StreamingError(
+                "fastperiod must be < slowperiod".to_string(),
+            ));
+        }
+        Ok(Self {
+            fast: EmaState::new(fastperiod),
+            slow: EmaState::new(slowperiod),
+            signal: EmaState::new(signalperiod),
+        })
+    }
+
+    /// Add a new close; return (macd_line, signal_line, histogram).
+    pub fn update(&mut self, value: f64) -> (f64, f64, f64) {
+        let fast_val = self.fast.update(value);
+        let slow_val = self.slow.update(value);
+
+        if slow_val.is_nan() {
+            return (f64::NAN, f64::NAN, f64::NAN);
+        }
+
+        let macd = fast_val - slow_val;
+        let signal = self.signal.update(macd);
+        if signal.is_nan() {
+            return (macd, f64::NAN, f64::NAN);
+        }
+        (macd, signal, macd - signal)
+    }
+
+    pub fn reset(&mut self) {
+        self.fast.reset();
+        self.slow.reset();
+        self.signal.reset();
+    }
+
+    pub fn fast_period(&self) -> usize {
+        self.fast.period()
+    }
+
+    pub fn slow_period(&self) -> usize {
+        self.slow.period()
+    }
+
+    pub fn signal_period(&self) -> usize {
+        self.signal.period()
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingStoch
+// ---------------------------------------------------------------------------
+
+/// Slow Stochastic (SMA-smoothed).
+///
+/// Returns (slowk, slowd).
+/// NaN tuple during warmup.
+pub struct StreamingStoch {
+    fastk_period: usize,
+    slowk_period: usize,
+    slowd_period: usize,
+    high_buf: VecDeque<f64>,
+    low_buf: VecDeque<f64>,
+    fastk_buf: VecDeque<f64>,
+    slowk_buf: VecDeque<f64>,
+}
+
+impl StreamingStoch {
+    pub fn new(
+        fastk_period: usize,
+        slowk_period: usize,
+        slowd_period: usize,
+    ) -> Result<Self, StreamingError> {
+        validate_timeperiod(fastk_period, "fastk_period", 1)?;
+        validate_timeperiod(slowk_period, "slowk_period", 1)?;
+        validate_timeperiod(slowd_period, "slowd_period", 1)?;
+        Ok(Self {
+            fastk_period,
+            slowk_period,
+            slowd_period,
+            high_buf: VecDeque::with_capacity(fastk_period + 1),
+            low_buf: VecDeque::with_capacity(fastk_period + 1),
+            fastk_buf: VecDeque::with_capacity(slowk_period + 1),
+            slowk_buf: VecDeque::with_capacity(slowd_period + 1),
+        })
+    }
+
+    /// Add a new bar (high, low, close); return (slowk, slowd).
+    pub fn update(&mut self, high: f64, low: f64, close: f64) -> (f64, f64) {
+        if self.high_buf.len() == self.fastk_period {
+            self.high_buf.pop_front();
+            self.low_buf.pop_front();
+        }
+        self.high_buf.push_back(high);
+        self.low_buf.push_back(low);
+
+        if self.high_buf.len() < self.fastk_period {
+            return (f64::NAN, f64::NAN);
+        }
+
+        let max_h = self
+            .high_buf
+            .iter()
+            .cloned()
+            .fold(f64::NEG_INFINITY, f64::max);
+        let min_l = self.low_buf.iter().cloned().fold(f64::INFINITY, f64::min);
+
+        let fastk = if max_h != min_l {
+            100.0 * (close - min_l) / (max_h - min_l)
+        } else {
+            0.0
+        };
+
+        if self.fastk_buf.len() == self.slowk_period {
+            self.fastk_buf.pop_front();
+        }
+        self.fastk_buf.push_back(fastk);
+        if self.fastk_buf.len() < self.slowk_period {
+            return (f64::NAN, f64::NAN);
+        }
+
+        let slowk = self.fastk_buf.iter().sum::<f64>() / self.slowk_period as f64;
+
+        if self.slowk_buf.len() == self.slowd_period {
+            self.slowk_buf.pop_front();
+        }
+        self.slowk_buf.push_back(slowk);
+        if self.slowk_buf.len() < self.slowd_period {
+            return (slowk, f64::NAN);
+        }
+
+        let slowd = self.slowk_buf.iter().sum::<f64>() / self.slowd_period as f64;
+        (slowk, slowd)
+    }
+
+    pub fn reset(&mut self) {
+        self.high_buf.clear();
+        self.low_buf.clear();
+        self.fastk_buf.clear();
+        self.slowk_buf.clear();
+    }
+
+    pub fn period(&self) -> usize {
+        self.fastk_period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingVWAP
+// ---------------------------------------------------------------------------
+
+/// Cumulative Volume Weighted Average Price.
+///
+/// Resets automatically when `reset()` is called (e.g. at session open).
+/// Accepts (high, low, close, volume) per bar.
+#[derive(Default)]
+pub struct StreamingVWAP {
+    cum_tpv: f64,
+    cum_vol: f64,
+}
+
+impl StreamingVWAP {
+    pub fn new() -> Self {
+        Self {
+            cum_tpv: 0.0,
+            cum_vol: 0.0,
+        }
+    }
+
+    /// Add a new bar (high, low, close, volume) and return cumulative VWAP.
+    pub fn update(&mut self, high: f64, low: f64, close: f64, volume: f64) -> f64 {
+        let tp = (high + low + close) / 3.0;
+        self.cum_tpv += tp * volume;
+        self.cum_vol += volume;
+        if self.cum_vol == 0.0 {
+            f64::NAN
+        } else {
+            self.cum_tpv / self.cum_vol
+        }
+    }
+
+    /// Reset for a new session.
+    pub fn reset(&mut self) {
+        self.cum_tpv = 0.0;
+        self.cum_vol = 0.0;
+    }
+}
+
+// ---------------------------------------------------------------------------
+// StreamingSupertrend
+// ---------------------------------------------------------------------------
+
+/// ATR-based Supertrend — streaming variant.
+///
+/// Accepts (high, low, close) per bar.
+/// Returns (supertrend_line, direction).
+/// direction: 1 = uptrend, -1 = downtrend, 0 = warmup.
+pub struct StreamingSupertrend {
+    period: usize,
+    multiplier: f64,
+    atr: AtrState,
+    upper_band: f64,
+    lower_band: f64,
+    has_bands: bool,
+    direction: i8,
+    prev_close: f64,
+    has_prev: bool,
+}
+
+impl StreamingSupertrend {
+    pub fn new(period: usize, multiplier: f64) -> Result<Self, StreamingError> {
+        validate_timeperiod(period, "period", 1)?;
+        Ok(Self {
+            period,
+            multiplier,
+            atr: AtrState::new(period),
+            upper_band: 0.0,
+            lower_band: 0.0,
+            has_bands: false,
+            direction: 0,
+            prev_close: 0.0,
+            has_prev: false,
+        })
+    }
+
+    /// Add a new bar (high, low, close); return (supertrend_line, direction).
+    pub fn update(&mut self, high: f64, low: f64, close: f64) -> (f64, i8) {
+        let atr = self.atr.update(high, low, close);
+        if atr.is_nan() {
+            self.prev_close = close;
+            self.has_prev = true;
+            return (f64::NAN, 0);
+        }
+
+        let hl2 = (high + low) / 2.0;
+        let upper_basic = hl2 + self.multiplier * atr;
+        let lower_basic = hl2 - self.multiplier * atr;
+
+        if !self.has_bands {
+            self.upper_band = upper_basic;
+            self.lower_band = lower_basic;
+            self.has_bands = true;
+            self.direction = -1;
+            self.prev_close = close;
+            self.has_prev = true;
+            return (self.upper_band, self.direction);
+        }
+
+        let prev_close = self.prev_close;
+
+        let new_lower = if lower_basic > self.lower_band || prev_close < self.lower_band {
+            lower_basic
+        } else {
+            self.lower_band
+        };
+        let new_upper = if upper_basic < self.upper_band || prev_close > self.upper_band {
+            upper_basic
+        } else {
+            self.upper_band
+        };
+
+        self.lower_band = new_lower;
+        self.upper_band = new_upper;
+
+        self.direction = if self.direction == -1 {
+            if close > new_upper {
+                1
+            } else {
+                -1
+            }
+        } else if close < new_lower {
+            -1
+        } else {
+            1
+        };
+
+        self.prev_close = close;
+        let line = if self.direction == 1 {
+            new_lower
+        } else {
+            new_upper
+        };
+        (line, self.direction)
+    }
+
+    pub fn reset(&mut self) {
+        self.atr.reset();
+        self.upper_band = 0.0;
+        self.lower_band = 0.0;
+        self.has_bands = false;
+        self.direction = 0;
+        self.prev_close = 0.0;
+        self.has_prev = false;
+    }
+
+    pub fn period(&self) -> usize {
+        self.period
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// Helper: compare two f64 values, treating NaN == NaN as true.
+    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
+        if a.is_nan() && b.is_nan() {
+            return true;
+        }
+        (a - b).abs() < tol
+    }
+
+    #[test]
+    fn test_sma_basic() {
+        let mut sma = StreamingSMA::new(3).unwrap();
+        assert!(sma.update(1.0).is_nan());
+        assert!(sma.update(2.0).is_nan());
+        let v = sma.update(3.0);
+        assert!(approx_eq(v, 2.0, 1e-10));
+        let v = sma.update(4.0);
+        assert!(approx_eq(v, 3.0, 1e-10));
+        let v = sma.update(5.0);
+        assert!(approx_eq(v, 4.0, 1e-10));
+        assert_eq!(sma.period(), 3);
+    }
+
+    #[test]
+    fn test_sma_reset() {
+        let mut sma = StreamingSMA::new(2).unwrap();
+        sma.update(10.0);
+        sma.update(20.0);
+        sma.reset();
+        assert!(sma.update(5.0).is_nan());
+        let v = sma.update(7.0);
+        assert!(approx_eq(v, 6.0, 1e-10));
+    }
+
+    #[test]
+    fn test_ema_warmup_and_decay() {
+        let mut ema = StreamingEMA::new(3).unwrap();
+        assert!(ema.update(2.0).is_nan());
+        assert!(ema.update(4.0).is_nan());
+        // Third bar: SMA seed = (2+4+6)/3 = 4.0
+        let v = ema.update(6.0);
+        assert!(approx_eq(v, 4.0, 1e-10));
+        // Fourth bar: alpha = 0.5, ema = 4.0 + 0.5*(8.0-4.0) = 6.0
+        let v = ema.update(8.0);
+        assert!(approx_eq(v, 6.0, 1e-10));
+    }
+
+    #[test]
+    fn test_rsi_warmup() {
+        let mut rsi = StreamingRSI::new(3).unwrap();
+        // First bar: no prev
+        assert!(rsi.update(44.0).is_nan());
+        // Bars 2-4: collecting gains/losses
+        assert!(rsi.update(44.5).is_nan());
+        assert!(rsi.update(43.5).is_nan());
+        // Bar 5: seeded
+        let v = rsi.update(44.5);
+        assert!(!v.is_nan());
+        assert!(v >= 0.0 && v <= 100.0);
+    }
+
+    #[test]
+    fn test_atr_warmup() {
+        let mut atr = StreamingATR::new(3).unwrap();
+        // First 3 bars return NaN (period = 3, seed happens on bar 3 but still NaN)
+        assert!(atr.update(10.0, 9.0, 9.5).is_nan());
+        assert!(atr.update(11.0, 9.5, 10.5).is_nan());
+        assert!(atr.update(10.5, 9.0, 9.5).is_nan());
+        // Bar 4: first real value
+        let v = atr.update(11.0, 10.0, 10.5);
+        assert!(!v.is_nan());
+        assert!(v > 0.0);
+    }
+
+    #[test]
+    fn test_bbands_warmup() {
+        let mut bb = StreamingBBands::new(3, 2.0, 2.0).unwrap();
+        let (u, m, l) = bb.update(10.0);
+        assert!(u.is_nan() && m.is_nan() && l.is_nan());
+        let (u, m, l) = bb.update(11.0);
+        assert!(u.is_nan() && m.is_nan() && l.is_nan());
+        let (u, m, l) = bb.update(12.0);
+        assert!(!u.is_nan() && !m.is_nan() && !l.is_nan());
+        assert!(approx_eq(m, 11.0, 1e-10));
+        assert!(u > m && l < m);
+    }
+
+    #[test]
+    fn test_macd_basic() {
+        let mut macd = StreamingMACD::new(3, 5, 2).unwrap();
+        // Feed enough bars for the slow (5) to seed
+        for i in 0..4 {
+            let (m, s, h) = macd.update(100.0 + i as f64);
+            assert!(m.is_nan());
+        }
+        // Bar 5: slow seeds
+        let (m, s, _h) = macd.update(104.0);
+        assert!(!m.is_nan());
+    }
+
+    #[test]
+    fn test_macd_fast_ge_slow_rejected() {
+        assert!(StreamingMACD::new(5, 3, 2).is_err());
+        assert!(StreamingMACD::new(5, 5, 2).is_err());
+    }
+
+    #[test]
+    fn test_stoch_basic() {
+        let mut stoch = StreamingStoch::new(3, 2, 2).unwrap();
+        // Need fastk_period bars, then slowk_period, then slowd_period
+        let (k, d) = stoch.update(10.0, 8.0, 9.0);
+        assert!(k.is_nan() && d.is_nan());
+        let (k, d) = stoch.update(11.0, 9.0, 10.0);
+        assert!(k.is_nan() && d.is_nan());
+        // Bar 3: fastk ready, collecting slowk
+        let (k, d) = stoch.update(12.0, 10.0, 11.0);
+        assert!(k.is_nan());
+        // Bar 4
+        let (k, d) = stoch.update(13.0, 11.0, 12.0);
+        assert!(!k.is_nan());
+    }
+
+    #[test]
+    fn test_vwap_basic() {
+        let mut vwap = StreamingVWAP::new();
+        let v = vwap.update(10.0, 8.0, 9.0, 100.0);
+        // tp = (10+8+9)/3 = 9.0, vwap = 9.0*100/100 = 9.0
+        assert!(approx_eq(v, 9.0, 1e-10));
+        let v = vwap.update(12.0, 10.0, 11.0, 200.0);
+        // tp2 = 11.0, cum_tpv = 900+2200=3100, cum_vol=300, vwap=10.333..
+        assert!(approx_eq(v, 3100.0 / 300.0, 1e-10));
+    }
+
+    #[test]
+    fn test_vwap_zero_volume() {
+        let mut vwap = StreamingVWAP::new();
+        let v = vwap.update(10.0, 8.0, 9.0, 0.0);
+        assert!(v.is_nan());
+    }
+
+    #[test]
+    fn test_supertrend_warmup() {
+        let mut st = StreamingSupertrend::new(3, 2.0).unwrap();
+        let (line, dir) = st.update(10.0, 9.0, 9.5);
+        assert!(line.is_nan() && dir == 0);
+        let (line, dir) = st.update(11.0, 9.5, 10.5);
+        assert!(line.is_nan() && dir == 0);
+        let (line, dir) = st.update(10.5, 9.0, 9.5);
+        assert!(line.is_nan() && dir == 0);
+        // Bar 4: first real value
+        let (line, dir) = st.update(11.0, 10.0, 10.5);
+        assert!(!line.is_nan());
+        assert!(dir == 1 || dir == -1);
+    }
+
+    #[test]
+    fn test_streaming_sma_matches_batch() {
+        // Compare streaming SMA against a simple batch computation
+        let data = vec![1.0, 3.0, 5.0, 7.0, 9.0, 11.0, 13.0];
+        let period = 3;
+        let mut sma = StreamingSMA::new(period).unwrap();
+        let streaming: Vec<f64> = data.iter().map(|&v| sma.update(v)).collect();
+
+        // Batch SMA
+        for i in 0..data.len() {
+            if i + 1 < period {
+                assert!(streaming[i].is_nan(), "bar {} should be NaN", i);
+            } else {
+                let batch: f64 =
+                    data[i + 1 - period..=i].iter().sum::<f64>() / period as f64;
+                assert!(
+                    approx_eq(streaming[i], batch, 1e-10),
+                    "bar {}: streaming={} batch={}",
+                    i,
+                    streaming[i],
+                    batch
+                );
+            }
+        }
+    }
+}

crates/ferro_ta_core/src/volatility.rs

@@ -1,10 +1,15 @@
 //! Volatility indicators.
 
-/// Average True Range — Wilder smoothed (TA-Lib compatible).
+/// Compute the Average True Range (ATR), Wilder smoothed (TA-Lib compatible).
 ///
-/// Seeds ATR with SMA of TR[1..=timeperiod] (bar 0 is skipped, matching TA-Lib).
-/// First valid output is at index `timeperiod`; indices 0..timeperiod are NaN.
-/// TR is computed on-the-fly (no separate tr Vec allocation).
+/// ATR measures market volatility by smoothing the True Range with Wilder's
+/// method. Seeded with the SMA of `TR[1..=timeperiod]` (bar 0 is skipped,
+/// matching TA-Lib). Returns non-negative values; the first `timeperiod`
+/// indices are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `timeperiod` - Smoothing period (typically 14).
 pub fn atr(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f64> {
     let n = high.len();
     let mut result = vec![f64::NAN; n];
@@ -33,7 +38,14 @@ pub fn atr(high: &[f64], low: &[f64], close: &[f64], timeperiod: usize) -> Vec<f
     result
 }
 
-/// True Range — max(H-L, |H-Cprev|, |L-Cprev|).
+/// Compute the True Range for each bar.
+///
+/// `TR = max(H - L, |H - C_prev|, |L - C_prev|)`. For bar 0, TR is
+/// simply `H - L` (no previous close available). Returns non-negative
+/// values for every bar (no `NaN` warmup).
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
 pub fn trange(high: &[f64], low: &[f64], close: &[f64]) -> Vec<f64> {
     let n = high.len();
     let mut result = vec![f64::NAN; n];

crates/ferro_ta_core/src/volume.rs

@@ -1,6 +1,14 @@
 //! Volume indicators.
 
-/// On-Balance Volume.
+/// Compute On-Balance Volume (OBV).
+///
+/// OBV is a cumulative indicator that adds volume on up-close bars and
+/// subtracts volume on down-close bars. Unchanged closes contribute zero.
+/// Returns a `Vec<f64>` of length `n` with no `NaN` values.
+///
+/// # Arguments
+/// * `close` - Price series.
+/// * `volume` - Volume series (same length as `close`).
 pub fn obv(close: &[f64], volume: &[f64]) -> Vec<f64> {
     let n = close.len();
     let mut result = vec![0.0_f64; n];
@@ -21,11 +29,17 @@ pub fn obv(close: &[f64], volume: &[f64]) -> Vec<f64> {
     result
 }
 
-/// Money Flow Index — O(n) sliding-window implementation without per-bar allocation.
+/// Compute the Money Flow Index (MFI).
 ///
-/// MFI = 100 - 100 / (1 + positive_flow / negative_flow) over `timeperiod` bars.
-/// typical_price = (high + low + close) / 3; raw_money_flow = typical_price * volume.
-/// Leading `timeperiod` values are NaN.
+/// MFI is a volume-weighted RSI, returning values in `[0, 100]`.
+/// `typical_price = (H + L + C) / 3`; money flow is positive when
+/// typical price rises, negative when it falls. The first `timeperiod`
+/// values are `NaN`.
+///
+/// # Arguments
+/// * `high` / `low` / `close` - OHLC price series (same length).
+/// * `volume` - Volume series (same length).
+/// * `timeperiod` - Lookback window (typically 14).
 pub fn mfi(
     high: &[f64],
     low: &[f64],
@@ -78,6 +92,51 @@ pub fn mfi(
     result
 }
 
+/// Chaikin Accumulation/Distribution Line.
+///
+/// Cumulates `(close - low - (high - close)) / (high - low) * volume`.
+pub fn ad(high: &[f64], low: &[f64], close: &[f64], volume: &[f64]) -> Vec<f64> {
+    let n = high.len();
+    let mut result = vec![0.0_f64; n];
+    let mut ad_val = 0.0_f64;
+    for i in 0..n {
+        let hl = high[i] - low[i];
+        let clv = if hl != 0.0 {
+            ((close[i] - low[i]) - (high[i] - close[i])) / hl
+        } else {
+            0.0
+        };
+        ad_val += clv * volume[i];
+        result[i] = ad_val;
+    }
+    result
+}
+
+/// Chaikin A/D Oscillator: fast EMA of AD minus slow EMA of AD.
+///
+/// Uses the core EMA implementation from `overlap::ema`.
+pub fn adosc(
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    volume: &[f64],
+    fastperiod: usize,
+    slowperiod: usize,
+) -> Vec<f64> {
+    let n = high.len();
+    let ad_vals = ad(high, low, close, volume);
+    let fast_ema = crate::overlap::ema(&ad_vals, fastperiod);
+    let slow_ema = crate::overlap::ema(&ad_vals, slowperiod);
+    let warmup = slowperiod - 1;
+    let mut result = vec![f64::NAN; n];
+    for i in warmup..n {
+        if !fast_ema[i].is_nan() && !slow_ema[i].is_nan() {
+            result[i] = fast_ema[i] - slow_ema[i];
+        }
+    }
+    result
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -92,6 +151,33 @@ mod tests {
         assert!((result[2] - 600.0).abs() < 1e-10);
     }
 
+    #[test]
+    fn ad_basic() {
+        let h = vec![10.0, 12.0, 11.0];
+        let l = vec![8.0, 9.0, 9.0];
+        let c = vec![9.0, 11.0, 10.0];
+        let v = vec![1000.0, 2000.0, 1500.0];
+        let result = ad(&h, &l, &c, &v);
+        assert_eq!(result.len(), 3);
+        // CLV[0] = ((9-8) - (10-9)) / (10-8) = (1 - 1) / 2 = 0
+        assert!((result[0] - 0.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn adosc_basic() {
+        let n = 30;
+        let h: Vec<f64> = (1..=n).map(|i| i as f64 + 1.0).collect();
+        let l: Vec<f64> = (1..=n).map(|i| i as f64 - 1.0).collect();
+        let c: Vec<f64> = (1..=n).map(|i| i as f64).collect();
+        let v: Vec<f64> = vec![1000.0; n];
+        let result = adosc(&h, &l, &c, &v, 3, 10);
+        assert_eq!(result.len(), n);
+        // Warmup period should be NaN
+        for i in 0..9 {
+            assert!(result[i].is_nan());
+        }
+    }
+
     #[test]
     fn mfi_range() {
         let n = 50;