Add 10 pairwise stat-arb indicators to Price Statistics (#154)

Adds ten pairwise `(f64, f64)` indicators to the **Price Statistics** family, completing the A1 stat-arb expansion block.

## Indicators

**Scalar output:**
- **RollingCorrelation** — rolling Pearson correlation of period-over-period *returns* (distinct from level-based `PearsonCorrelation`).
- **RollingCovariance** — rolling covariance of returns.
- **OuHalfLife** — Ornstein–Uhlenbeck half-life of mean reversion of the spread `a − b`.
- **SpreadHurst** — Hurst exponent of the spread (variance-of-lagged-differences fit) for regime detection.
- **DistanceSsd** — Gatev sum-of-squared-deviations between two start-normalised series.
- **BetaNeutralSpread** — rolling OLS regression residual `a − (α + β·b)`.
- **VarianceRatio** — Lo–MacKinlay variance-ratio test on the spread (two params: `period`, `q`).
- **GrangerCausality** — F-statistic for whether `b` predicts `a` (two params: `period`, `lag`).

**Struct output (custom bindings):**
- **KalmanHedgeRatio** — dynamic hedge ratio via a Kalman filter → `{ hedgeRatio, intercept, spread }`.
- **SpreadBollingerBands** — Bollinger bands on the spread → `{ middle, upper, lower, percentB }`.

## Notes
- No new traits or input families: all use the native `Indicator<Input = (f64, f64)>` (precedent `Beta`, `Cointegration`).
- Adds `Error::InvalidParameter` for floating-point constructor parameters (Kalman `delta`/`observation_var`, `num_std`).
- Full Python/Node/WASM bindings; the two struct-output indicators are hand-written, the rest use the pair macros.
- Indicator count 315 → 325; README, family rows, `__init__`, fuzz target, and CHANGELOG updated.

## Verification
- `cargo test --workspace --all-features` — green (2676 core lib + 308 doc).
- `cargo clippy --workspace --all-targets --all-features -- -D warnings` — clean.
- Node: `npm run build && npm test` — 410 passing (`index.d.ts`/`index.js` regenerated).
- Python: `pytest` — 684 passing.
This commit is contained in:
kingchenc
2026-06-03 15:39:55 +02:00
committed by GitHub
parent 53941b7b07
commit a3a1ae4dba
25 changed files with 4313 additions and 51 deletions
@@ -0,0 +1,247 @@
//! Beta-neutral spread: the rolling OLS regression residual of two series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// The beta-neutral spread between two assets — the residual of a rolling
/// ordinary-least-squares regression of `a` on `b`.
///
/// Each `update` takes one `(a, b)` price pair. Over the trailing window of
/// `period` pairs the indicator fits the hedge ratio `β` (and intercept `α`) by
/// OLS and reports the **current** residual:
///
/// ```text
/// β = cov(a, b) / var(b) α = ā β · b̄
/// spread = a_now (α + β · b_now)
/// ```
///
/// Subtracting `β · b` removes `a`'s exposure to `b`, so the spread is market-
/// (beta-)neutral: it is what is left after the common factor is hedged out.
/// Positive means `a` is rich relative to its hedge, negative means cheap — the
/// raw signal a pairs trade fades. Where [`crate::PairSpreadZScore`] standardises
/// this residual into a z-score and [`crate::Cointegration`] bundles it with an
/// ADF test, this indicator returns the residual itself, in price units.
///
/// If `b` is flat over the window (`var(b) = 0`) there is no defined slope; the
/// indicator falls back to `β = 0`, so the spread becomes `a_now ā`.
///
/// Each `update` is `O(1)`: four running sums (`Σa`, `Σb`, `Σb²`, `Σab`) are
/// maintained as the window slides.
///
/// # Example
///
/// ```
/// use wickra_core::{BetaNeutralSpread, Indicator};
///
/// let mut s = BetaNeutralSpread::new(20).unwrap();
/// let mut last = None;
/// for t in 0..40 {
/// let b = 100.0 + f64::from(t);
/// // a = 2·b + 5 exactly ⇒ the regression explains a fully ⇒ spread ≈ 0.
/// last = s.update((2.0 * b + 5.0, b));
/// }
/// assert!(last.unwrap().abs() < 1e-6);
/// ```
#[derive(Debug, Clone)]
pub struct BetaNeutralSpread {
period: usize,
window: VecDeque<(f64, f64)>,
sum_a: f64,
sum_b: f64,
sum_bb: f64,
sum_ab: f64,
}
impl BetaNeutralSpread {
/// Construct a new beta-neutral spread.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 2` — a regression slope
/// needs at least two points.
pub fn new(period: usize) -> Result<Self> {
if period < 2 {
return Err(Error::InvalidPeriod {
message: "beta-neutral spread needs period >= 2",
});
}
Ok(Self {
period,
window: VecDeque::with_capacity(period),
sum_a: 0.0,
sum_b: 0.0,
sum_bb: 0.0,
sum_ab: 0.0,
})
}
/// Configured look-back window.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for BetaNeutralSpread {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (a, b) = input;
if self.window.len() == self.period {
let (oa, ob) = self.window.pop_front().expect("non-empty");
self.sum_a -= oa;
self.sum_b -= ob;
self.sum_bb -= ob * ob;
self.sum_ab -= oa * ob;
}
self.window.push_back((a, b));
self.sum_a += a;
self.sum_b += b;
self.sum_bb += b * b;
self.sum_ab += a * b;
if self.window.len() < self.period {
return None;
}
let n = self.period as f64;
let mean_a = self.sum_a / n;
let mean_b = self.sum_b / n;
let var_b = (self.sum_bb / n - mean_b * mean_b).max(0.0);
let (beta, intercept) = if var_b == 0.0 {
(0.0, mean_a)
} else {
let cov = self.sum_ab / n - mean_a * mean_b;
let slope = cov / var_b;
(slope, mean_a - slope * mean_b)
};
Some(a - (intercept + beta * b))
}
fn reset(&mut self) {
self.window.clear();
self.sum_a = 0.0;
self.sum_b = 0.0;
self.sum_bb = 0.0;
self.sum_ab = 0.0;
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"BetaNeutralSpread"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_period_below_two() {
assert!(BetaNeutralSpread::new(1).is_err());
assert!(BetaNeutralSpread::new(2).is_ok());
}
#[test]
fn accessors_and_metadata() {
let s = BetaNeutralSpread::new(20).unwrap();
assert_eq!(s.period(), 20);
assert_eq!(s.warmup_period(), 20);
assert_eq!(s.name(), "BetaNeutralSpread");
assert!(!s.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut s = BetaNeutralSpread::new(3).unwrap();
assert_eq!(s.update((1.0, 1.0)), None);
assert_eq!(s.update((2.0, 2.0)), None);
assert!(s.update((3.0, 3.0)).is_some());
assert!(s.is_ready());
}
#[test]
fn perfect_linear_relationship_has_zero_spread() {
let pairs: Vec<(f64, f64)> = (0..40)
.map(|t| {
let b = 100.0 + f64::from(t);
(2.0 * b + 5.0, b)
})
.collect();
let last = BetaNeutralSpread::new(20)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.0, epsilon = 1e-6);
}
#[test]
fn dislocation_produces_nonzero_spread() {
// a tracks 2·b, then the last bar jumps up ⇒ positive residual.
let mut pairs: Vec<(f64, f64)> = (0..19)
.map(|t| {
let b = 100.0 + f64::from(t);
(2.0 * b + 5.0, b)
})
.collect();
pairs.push((2.0 * 119.0 + 5.0 + 10.0, 119.0));
let last = BetaNeutralSpread::new(20)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last > 1.0, "spread {last}");
}
#[test]
fn flat_b_falls_back_to_demeaned_a() {
// b constant ⇒ β = 0 ⇒ spread = a mean(a). Last window of a = 0..9,
// mean = 4.5, last a = 9 ⇒ spread = 4.5.
let pairs: Vec<(f64, f64)> = (0..10).map(|t| (f64::from(t), 7.0)).collect();
let last = BetaNeutralSpread::new(10)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 4.5, epsilon = 1e-12);
}
#[test]
fn reset_clears_state() {
let mut s = BetaNeutralSpread::new(4).unwrap();
s.batch(&[(1.0, 2.0), (2.0, 4.0), (3.0, 5.0), (4.0, 9.0), (5.0, 2.0)]);
assert!(s.is_ready());
s.reset();
assert!(!s.is_ready());
assert_eq!(s.update((1.0, 1.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..60)
.map(|t| {
let b = 30.0 + 0.7 * f64::from(t);
(1.8 * b + 2.0 + (f64::from(t) * 0.4).sin(), b)
})
.collect();
let batch = BetaNeutralSpread::new(20).unwrap().batch(&pairs);
let mut s = BetaNeutralSpread::new(20).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| s.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,235 @@
//! Gatev distance (sum of squared deviations) between two normalised series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Sum of squared deviations between two price series, normalised to a common
/// start — the classic Gatev et al. pairs-selection distance.
///
/// Each `update` takes one `(a, b)` price pair. Over the trailing window of
/// `period` pairs each series is rebased to `1` at the window's first bar and
/// the squared gap between the two normalised paths is summed:
///
/// ```text
/// ãᵢ = aᵢ / a_first b̃ᵢ = bᵢ / b_first
/// SSD = Σ (ãᵢ b̃ᵢ)²
/// ```
///
/// Rebasing puts the two series on the same scale (both start at `1`), so the
/// distance measures how far their *relative* paths drift apart. A **small**
/// SSD means the two assets track each other tightly — the screen Gatev,
/// Goetzmann and Rouwenhorst use to pick tradeable pairs; a large SSD means
/// they have decoupled. The output is always `≥ 0`. If either series is `0` at
/// the start of the window the normalisation is undefined and the indicator
/// returns `0`.
///
/// Each `update` is `O(period)`, bounded by the fixed window.
///
/// # Example
///
/// ```
/// use wickra_core::{DistanceSsd, Indicator};
///
/// let mut d = DistanceSsd::new(20).unwrap();
/// let mut last = None;
/// for t in 0..40 {
/// let base = 100.0 + f64::from(t);
/// // Two near-identical paths ⇒ tiny distance.
/// last = d.update((base, base * 1.0001));
/// }
/// assert!(last.unwrap() < 1e-3);
/// ```
#[derive(Debug, Clone)]
pub struct DistanceSsd {
period: usize,
window: VecDeque<(f64, f64)>,
}
impl DistanceSsd {
/// Construct a new Gatev distance estimator.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 2` — a distance needs at
/// least two points.
pub fn new(period: usize) -> Result<Self> {
if period < 2 {
return Err(Error::InvalidPeriod {
message: "distance SSD needs period >= 2",
});
}
Ok(Self {
period,
window: VecDeque::with_capacity(period),
})
}
/// Configured look-back window.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for DistanceSsd {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
if self.window.len() == self.period {
self.window.pop_front();
}
self.window.push_back(input);
if self.window.len() < self.period {
return None;
}
let &(a_first, b_first) = self.window.front().expect("window is full");
if a_first == 0.0 || b_first == 0.0 {
// Cannot rebase a series that starts at zero.
return Some(0.0);
}
let ssd = self
.window
.iter()
.map(|&(a, b)| {
let gap = a / a_first - b / b_first;
gap * gap
})
.sum();
Some(ssd)
}
fn reset(&mut self) {
self.window.clear();
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"DistanceSsd"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_period_below_two() {
assert!(DistanceSsd::new(1).is_err());
assert!(DistanceSsd::new(2).is_ok());
}
#[test]
fn accessors_and_metadata() {
let d = DistanceSsd::new(20).unwrap();
assert_eq!(d.period(), 20);
assert_eq!(d.warmup_period(), 20);
assert_eq!(d.name(), "DistanceSsd");
assert!(!d.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut d = DistanceSsd::new(3).unwrap();
assert_eq!(d.update((1.0, 1.0)), None);
assert_eq!(d.update((2.0, 2.0)), None);
assert!(d.update((3.0, 3.0)).is_some());
assert!(d.is_ready());
}
#[test]
fn identical_normalised_paths_have_zero_distance() {
// b = 2·a ⇒ both rebase to the same path ⇒ SSD = 0.
let pairs: Vec<(f64, f64)> = (0..20)
.map(|t| {
let a = 100.0 + f64::from(t);
(a, 2.0 * a)
})
.collect();
let last = DistanceSsd::new(10)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.0, epsilon = 1e-12);
}
#[test]
fn diverging_paths_have_positive_distance() {
let pairs: Vec<(f64, f64)> = (0..20)
.map(|t| (100.0 + f64::from(t), 100.0 + 3.0 * f64::from(t)))
.collect();
let last = DistanceSsd::new(10)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last > 0.0, "ssd {last}");
}
#[test]
fn hand_computed_value() {
// Window of three pairs, a_first = b_first = 1:
// (1,1) → 0; (2,4) → (24)² = 4; (3,9) → (39)² = 36 ⇒ SSD = 40.
let pairs = [(1.0, 1.0), (2.0, 4.0), (3.0, 9.0)];
let last = DistanceSsd::new(3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 40.0, epsilon = 1e-12);
}
#[test]
fn zero_start_returns_zero() {
// First bar of the window has a = 0 ⇒ rebasing undefined ⇒ 0.
let pairs = [(0.0, 1.0), (2.0, 2.0), (3.0, 3.0)];
let last = DistanceSsd::new(3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 0.0);
}
#[test]
fn reset_clears_state() {
let mut d = DistanceSsd::new(4).unwrap();
d.batch(&[(1.0, 1.0), (2.0, 2.0), (3.0, 4.0), (4.0, 5.0), (5.0, 6.0)]);
assert!(d.is_ready());
d.reset();
assert!(!d.is_ready());
assert_eq!(d.update((1.0, 1.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..60)
.map(|t| {
let a = 100.0 + f64::from(t);
(a, 100.0 + 1.2 * f64::from(t) + (f64::from(t) * 0.5).sin())
})
.collect();
let batch = DistanceSsd::new(15).unwrap().batch(&pairs);
let mut d = DistanceSsd::new(15).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| d.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,335 @@
//! Granger causality F-statistic: does series `b` help predict series `a`?
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Granger causality of `b` on `a` over a rolling window, as an F-statistic.
///
/// Each `update` takes one `(a, b)` pair. Over the trailing window of `period`
/// observations the indicator fits two autoregressions of `a` and compares them
/// with an F-test:
///
/// ```text
/// restricted: aₜ = c + Σ φᵢ·aₜ₋ᵢ (a's own lags only)
/// unrestricted: aₜ = c + Σ φᵢ·aₜ₋ᵢ + Σ ψᵢ·bₜ₋ᵢ (+ b's lags)
/// F = ((RSSᵣ RSSᵤ) / lag) / (RSSᵤ / (n 2·lag 1))
/// ```
///
/// If adding `b`'s lags significantly reduces the residual sum of squares, `b`
/// **Granger-causes** `a`: past values of `b` carry information about the future
/// of `a` beyond what `a`'s own past holds. A **larger** F means stronger
/// predictive causality (leadlag structure a stat-arb model can trade); a
/// value near `0` means `b` adds nothing. Note Granger causality is purely
/// predictive — it is not structural cause and effect.
///
/// The statistic is `0` when a regression is degenerate — a collinear or flat
/// window makes the normal equations singular. The output is always `≥ 0`.
///
/// Each `update` is `O(period · lag² + lag³)`, bounded by the fixed parameters.
///
/// # Example
///
/// ```
/// use wickra_core::{GrangerCausality, Indicator};
///
/// let mut g = GrangerCausality::new(60, 1).unwrap();
/// let mut last = None;
/// for t in 0..120 {
/// let drive = (f64::from(t) * 0.3).sin();
/// // a echoes b's previous value plus noise ⇒ b Granger-causes a.
/// let b = drive;
/// let a = 0.5 * (f64::from(t.max(1) - 1) * 0.3).sin() + 0.1 * (f64::from(t) * 0.9).cos();
/// last = g.update((a, b));
/// }
/// assert!(last.unwrap() >= 0.0);
/// ```
#[derive(Debug, Clone)]
pub struct GrangerCausality {
period: usize,
lag: usize,
window: VecDeque<(f64, f64)>,
}
impl GrangerCausality {
/// Construct a new Granger causality test.
///
/// `period` is the look-back window; `lag` is the autoregressive order
/// (number of own/cross lags in each model).
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `lag < 1` or if `period < 3·lag + 2`
/// (the smallest window that leaves the unrestricted regression at least one
/// residual degree of freedom).
pub fn new(period: usize, lag: usize) -> Result<Self> {
if lag < 1 {
return Err(Error::InvalidPeriod {
message: "granger causality needs lag >= 1",
});
}
if period < 3 * lag + 2 {
return Err(Error::InvalidPeriod {
message: "granger causality needs period >= 3*lag + 2",
});
}
Ok(Self {
period,
lag,
window: VecDeque::with_capacity(period),
})
}
/// Configured look-back window.
pub const fn period(&self) -> usize {
self.period
}
/// Configured autoregressive order.
pub const fn lag(&self) -> usize {
self.lag
}
}
impl Indicator for GrangerCausality {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
if self.window.len() == self.period {
self.window.pop_front();
}
self.window.push_back(input);
if self.window.len() < self.period {
return None;
}
let lag = self.lag;
let a: Vec<f64> = self.window.iter().map(|&(av, _)| av).collect();
let b: Vec<f64> = self.window.iter().map(|&(_, bv)| bv).collect();
let num_obs = self.period - lag;
let mut target = Vec::with_capacity(num_obs);
let mut restricted = Vec::with_capacity(num_obs);
let mut unrestricted = Vec::with_capacity(num_obs);
for k in 0..num_obs {
let now = lag + k;
target.push(a[now]);
let mut row_r = Vec::with_capacity(lag + 1);
row_r.push(1.0);
for back in 1..=lag {
row_r.push(a[now - back]);
}
let mut row_u = row_r.clone();
for back in 1..=lag {
row_u.push(b[now - back]);
}
restricted.push(row_r);
unrestricted.push(row_u);
}
let Some(rss_r) = ols_rss(&restricted, &target, lag + 1) else {
return Some(0.0);
};
let Some(rss_u) = ols_rss(&unrestricted, &target, 2 * lag + 1) else {
return Some(0.0);
};
let dof = (num_obs - (2 * lag + 1)) as f64;
let numerator = (rss_r - rss_u) / lag as f64;
let denominator = rss_u / dof;
Some((numerator / denominator).max(0.0))
}
fn reset(&mut self) {
self.window.clear();
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"GrangerCausality"
}
}
/// Residual sum of squares of the OLS fit of `target` on the design `rows`
/// (each a length-`num_reg` regressor vector). Returns `None` if the normal
/// equations are singular.
fn ols_rss(rows: &[Vec<f64>], target: &[f64], num_reg: usize) -> Option<f64> {
let mut xtx = vec![vec![0.0; num_reg]; num_reg];
let mut xty = vec![0.0; num_reg];
for (row, &observed) in rows.iter().zip(target) {
for (ri, &left) in row.iter().enumerate() {
xty[ri] += left * observed;
for (ci, &right) in row.iter().enumerate() {
xtx[ri][ci] += left * right;
}
}
}
let theta = solve(xtx, xty)?;
let mut rss = 0.0;
for (row, &observed) in rows.iter().zip(target) {
let pred: f64 = row
.iter()
.zip(&theta)
.map(|(coeff, value)| coeff * value)
.sum();
let resid = observed - pred;
rss += resid * resid;
}
Some(rss)
}
/// Solve the linear system `mat·x = rhs` by Gaussian elimination, returning
/// `None` if the matrix is (numerically) singular. `mat` is row-major.
fn solve(mut mat: Vec<Vec<f64>>, mut rhs: Vec<f64>) -> Option<Vec<f64>> {
let dim = rhs.len();
for col in 0..dim {
let pivot = mat[col][col];
if pivot.abs() < 1e-12 {
return None;
}
let pivot_row = mat[col].clone();
for row in (col + 1)..dim {
let factor = mat[row][col] / pivot;
for (cell, &above) in mat[row].iter_mut().zip(&pivot_row).skip(col) {
*cell -= factor * above;
}
rhs[row] -= factor * rhs[col];
}
}
let mut sol = vec![0.0; dim];
for row in (0..dim).rev() {
let known: f64 = mat[row]
.iter()
.zip(&sol)
.skip(row + 1)
.map(|(coeff, value)| coeff * value)
.sum();
sol[row] = (rhs[row] - known) / mat[row][row];
}
Some(sol)
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
#[test]
fn rejects_bad_parameters() {
assert!(GrangerCausality::new(10, 0).is_err()); // lag must be >= 1
assert!(GrangerCausality::new(4, 1).is_err()); // period must be >= 3*lag + 2
assert!(GrangerCausality::new(5, 1).is_ok());
}
#[test]
fn accessors_and_metadata() {
let g = GrangerCausality::new(60, 2).unwrap();
assert_eq!(g.period(), 60);
assert_eq!(g.lag(), 2);
assert_eq!(g.warmup_period(), 60);
assert_eq!(g.name(), "GrangerCausality");
assert!(!g.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut g = GrangerCausality::new(5, 1).unwrap();
for t in 0..4 {
assert_eq!(g.update((f64::from(t), f64::from(t) * 0.5)), None);
}
assert!(g.update((4.0, 2.0)).is_some());
assert!(g.is_ready());
}
#[test]
fn b_leading_a_has_positive_statistic() {
// a[t] is driven by b[t-1] plus a little of its own past ⇒ b helps.
let mut prev_drive = 0.0;
let pairs: Vec<(f64, f64)> = (0..120)
.map(|t| {
let drive = (f64::from(t) * 0.3).sin() + 0.4 * (f64::from(t) * 0.11).cos();
let a = 0.8 * prev_drive + 0.05 * (f64::from(t) * 0.7).sin();
prev_drive = drive;
(a, drive)
})
.collect();
let last = GrangerCausality::new(60, 1)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last > 1.0, "F {last}");
}
#[test]
fn constant_b_is_singular_and_returns_zero() {
// b is constant ⇒ its lag columns are collinear with the intercept ⇒
// the unrestricted normal equations are singular ⇒ 0.
let pairs: Vec<(f64, f64)> = (0..40)
.map(|t| (f64::from(t) + (f64::from(t) * 0.6).sin(), 3.0))
.collect();
let last = GrangerCausality::new(20, 1)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 0.0);
}
#[test]
fn constant_a_restricted_singular_returns_zero() {
// a is constant ⇒ its own lag columns are collinear with the intercept
// ⇒ the restricted normal equations are singular ⇒ 0.
let pairs: Vec<(f64, f64)> = (0..40).map(|t| (5.0, (f64::from(t) * 0.4).sin())).collect();
let last = GrangerCausality::new(20, 1)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 0.0);
}
#[test]
fn reset_clears_state() {
let mut g = GrangerCausality::new(8, 1).unwrap();
for t in 0..12 {
g.update((
f64::from(t) + (f64::from(t) * 0.7).sin(),
(f64::from(t) * 0.3).cos(),
));
}
assert!(g.is_ready());
g.reset();
assert!(!g.is_ready());
assert_eq!(g.update((1.0, 1.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..80)
.map(|t| {
let b = (f64::from(t) * 0.4).sin();
(
0.6 * (f64::from(t.max(1) - 1) * 0.4).sin() + 0.1 * f64::from(t % 3),
b,
)
})
.collect();
let batch = GrangerCausality::new(30, 2).unwrap().batch(&pairs);
let mut g = GrangerCausality::new(30, 2).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| g.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,290 @@
//! Kalman-filter dynamic hedge ratio between two series.
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Output of [`KalmanHedgeRatio`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct KalmanHedgeRatioOutput {
/// Current hedge ratio `β` — the filtered slope of `a` on `b`.
pub hedge_ratio: f64,
/// Current intercept `α` — the filtered level offset.
pub intercept: f64,
/// Forecast error `a (α + β·b)`: how far the latest `a` sits from the
/// Kalman-predicted relationship. This is the tradeable spread signal.
pub spread: f64,
}
/// Dynamic hedge ratio between two series, estimated online with a Kalman filter.
///
/// Each `update` takes one `(a, b)` price pair and treats the linear relation
/// `aₜ = αₜ + βₜ·bₜ + noise` as a state-space model whose hidden state
/// `[βₜ, αₜ]` follows a random walk. The filter updates the state from every
/// observation, so the hedge ratio **adapts continuously** instead of being a
/// flat OLS slope over a fixed window:
///
/// ```text
/// state xₜ = [βₜ, αₜ], drifts as a random walk with covariance Vw·I
/// observe aₜ = [bₜ, 1]·xₜ + εₜ, Var(εₜ) = observation_var
/// Vw = delta / (1 delta)
/// ```
///
/// `delta` controls how fast the hedge ratio is allowed to move: a larger
/// `delta` tracks regime changes faster but is noisier; a smaller `delta` is
/// smoother but slower. `observation_var` is the measurement-noise variance.
/// The reported `spread` (the filter's forecast error) is the mean-reverting
/// signal a pairs trade fades — the Kalman analogue of the
/// [`crate::Cointegration`] residual, but with a hedge ratio that breathes.
///
/// The filter emits an estimate from the **first** update (warmup of one bar);
/// early estimates are diffuse and settle as observations accumulate. Each
/// `update` is `O(1)` over the fixed 2×2 covariance.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, KalmanHedgeRatio};
///
/// let mut k = KalmanHedgeRatio::new(1e-2, 1e-3).unwrap();
/// let mut last = None;
/// for t in 0..400 {
/// // `b` sweeps a wide range so the slope and intercept are identifiable.
/// let b = 100.0 + (f64::from(t) * 0.5).sin() * 95.0;
/// last = k.update((2.0 * b + 5.0, b)); // a = 2·b + 5
/// }
/// let out = last.unwrap();
/// assert!((out.hedge_ratio - 2.0).abs() < 0.05);
/// assert!(out.spread.abs() < 0.05);
/// ```
#[derive(Debug, Clone)]
pub struct KalmanHedgeRatio {
delta: f64,
transition_var: f64,
observation_var: f64,
beta: f64,
alpha: f64,
// State covariance, row-major 2×2: [[p00, p01], [p10, p11]].
cov: [[f64; 2]; 2],
count: usize,
}
impl KalmanHedgeRatio {
/// Construct a new Kalman hedge-ratio filter.
///
/// `delta` is the state-drift ratio in `(0, 1)`; `observation_var` is the
/// measurement-noise variance (`> 0`).
///
/// # Errors
/// Returns [`Error::InvalidParameter`] if `delta` is not in `(0, 1)` or if
/// `observation_var` is not strictly positive (both must also be finite).
pub fn new(delta: f64, observation_var: f64) -> Result<Self> {
if !delta.is_finite() || delta <= 0.0 || delta >= 1.0 {
return Err(Error::InvalidParameter {
message: "kalman hedge ratio needs delta in (0, 1)",
});
}
if !observation_var.is_finite() || observation_var <= 0.0 {
return Err(Error::InvalidParameter {
message: "kalman hedge ratio needs observation_var > 0",
});
}
Ok(Self {
delta,
transition_var: delta / (1.0 - delta),
observation_var,
beta: 0.0,
alpha: 0.0,
cov: [[0.0; 2]; 2],
count: 0,
})
}
/// Configured state-drift ratio `delta`.
pub const fn delta(&self) -> f64 {
self.delta
}
/// Configured measurement-noise variance.
pub const fn observation_var(&self) -> f64 {
self.observation_var
}
}
impl Indicator for KalmanHedgeRatio {
type Input = (f64, f64);
type Output = KalmanHedgeRatioOutput;
fn update(&mut self, input: (f64, f64)) -> Option<KalmanHedgeRatioOutput> {
let (a, b) = input;
// Predicted state covariance: add the transition noise to the diagonal
// (the very first observation starts from a zero prior).
let mut cov_pred = self.cov;
if self.count > 0 {
cov_pred[0][0] += self.transition_var;
cov_pred[1][1] += self.transition_var;
}
// Observation row is F = [b, 1].
let predicted = self.beta * b + self.alpha;
let innovation = a - predicted;
// F·cov_pred (a 1×2 row).
let fr0 = b * cov_pred[0][0] + cov_pred[1][0];
let fr1 = b * cov_pred[0][1] + cov_pred[1][1];
// Innovation variance S = F·cov_pred·Fᵀ + observation_var ≥ observation_var > 0.
let innovation_var = fr0 * b + fr1 + self.observation_var;
// Kalman gain = cov_pred·Fᵀ / S.
let rft0 = cov_pred[0][0] * b + cov_pred[0][1];
let rft1 = cov_pred[1][0] * b + cov_pred[1][1];
let gain0 = rft0 / innovation_var;
let gain1 = rft1 / innovation_var;
// State update.
self.beta += gain0 * innovation;
self.alpha += gain1 * innovation;
// Covariance update P = cov_pred gain·(F·cov_pred).
self.cov[0][0] = cov_pred[0][0] - gain0 * fr0;
self.cov[0][1] = cov_pred[0][1] - gain0 * fr1;
self.cov[1][0] = cov_pred[1][0] - gain1 * fr0;
self.cov[1][1] = cov_pred[1][1] - gain1 * fr1;
self.count += 1;
Some(KalmanHedgeRatioOutput {
hedge_ratio: self.beta,
intercept: self.alpha,
spread: innovation,
})
}
fn reset(&mut self) {
self.beta = 0.0;
self.alpha = 0.0;
self.cov = [[0.0; 2]; 2];
self.count = 0;
}
fn warmup_period(&self) -> usize {
1
}
fn is_ready(&self) -> bool {
self.count >= 1
}
fn name(&self) -> &'static str {
"KalmanHedgeRatio"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
#[test]
fn rejects_bad_parameters() {
assert!(KalmanHedgeRatio::new(0.0, 1.0).is_err());
assert!(KalmanHedgeRatio::new(1.0, 1.0).is_err());
assert!(KalmanHedgeRatio::new(-0.1, 1.0).is_err());
assert!(KalmanHedgeRatio::new(f64::NAN, 1.0).is_err());
assert!(KalmanHedgeRatio::new(0.001, 0.0).is_err());
assert!(KalmanHedgeRatio::new(0.001, -1.0).is_err());
assert!(KalmanHedgeRatio::new(0.001, f64::INFINITY).is_err());
assert!(KalmanHedgeRatio::new(0.001, 0.001).is_ok());
}
#[test]
fn accessors_and_metadata() {
let k = KalmanHedgeRatio::new(0.001, 0.01).unwrap();
assert_eq!(k.delta(), 0.001);
assert_eq!(k.observation_var(), 0.01);
assert_eq!(k.warmup_period(), 1);
assert_eq!(k.name(), "KalmanHedgeRatio");
assert!(!k.is_ready());
}
#[test]
fn emits_from_first_update() {
let mut k = KalmanHedgeRatio::new(0.001, 0.001).unwrap();
let first = k.update((10.0, 5.0)).unwrap();
// The diffuse prior leaves the first state at the origin.
assert_eq!(first.hedge_ratio, 0.0);
assert_eq!(first.intercept, 0.0);
assert_eq!(first.spread, 10.0);
assert!(k.is_ready());
}
#[test]
fn converges_to_static_relationship() {
// a = 2·b + 5 ⇒ the filter should recover β ≈ 2, α ≈ 5, spread ≈ 0.
// `b` sweeps a wide range so β and α are jointly identifiable.
let pairs: Vec<(f64, f64)> = (0..500)
.map(|t| {
let b = 100.0 + (f64::from(t) * 0.5).sin() * 95.0;
(2.0 * b + 5.0, b)
})
.collect();
let out = KalmanHedgeRatio::new(1e-2, 1e-3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(
(out.hedge_ratio - 2.0).abs() < 0.05,
"beta {}",
out.hedge_ratio
);
assert!((out.intercept - 5.0).abs() < 1.0, "alpha {}", out.intercept);
assert!(out.spread.abs() < 0.05, "spread {}", out.spread);
}
#[test]
fn tracks_a_changing_hedge_ratio() {
// Hedge ratio steps from 2 to 3 partway through; the filter should move
// toward the new ratio.
let mut pairs: Vec<(f64, f64)> = (0..300)
.map(|t| {
let b = 100.0 + (f64::from(t) * 0.5).sin() * 95.0;
(2.0 * b + 5.0, b)
})
.collect();
pairs.extend((0..300).map(|t| {
let b = 100.0 + (f64::from(t) * 0.5).cos() * 95.0;
(3.0 * b + 5.0, b)
}));
let out = KalmanHedgeRatio::new(1e-2, 1e-3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(out.hedge_ratio > 2.5, "beta {}", out.hedge_ratio);
}
#[test]
fn reset_clears_state() {
let mut k = KalmanHedgeRatio::new(0.001, 0.001).unwrap();
for t in 0..50 {
let b = 100.0 + f64::from(t);
k.update((2.0 * b, b));
}
assert!(k.is_ready());
k.reset();
assert!(!k.is_ready());
let first = k.update((10.0, 5.0)).unwrap();
assert_eq!(first.hedge_ratio, 0.0);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..120)
.map(|t| {
let b = 30.0 + 0.7 * f64::from(t);
(1.8 * b + 2.0 + (f64::from(t) * 0.4).sin(), b)
})
.collect();
let batch = KalmanHedgeRatio::new(1e-3, 1e-2).unwrap().batch(&pairs);
let mut k = KalmanHedgeRatio::new(1e-3, 1e-2).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| k.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
+31 -1
View File
@@ -33,6 +33,7 @@ mod awesome_oscillator_histogram;
mod balance_of_power;
mod belt_hold;
mod beta;
mod beta_neutral_spread;
mod bollinger;
mod bollinger_bandwidth;
mod breakaway;
@@ -67,6 +68,7 @@ mod demand_index;
mod demark_pivots;
mod depth_slope;
mod detrended_std_dev;
mod distance_ssd;
mod doji;
mod doji_star;
mod donchian;
@@ -101,6 +103,7 @@ mod funding_rate_zscore;
mod gain_loss_ratio;
mod gap_side_by_side_white;
mod garman_klass;
mod granger_causality;
mod gravestone_doji;
mod hammer;
mod hanging_man;
@@ -130,6 +133,7 @@ mod inverse_fisher_transform;
mod inverted_hammer;
mod jma;
mod kagi_bars;
mod kalman_hedge_ratio;
mod kama;
mod kelly_criterion;
mod keltner;
@@ -187,6 +191,7 @@ mod omega_ratio;
mod on_neck;
mod opening_marubozu;
mod opening_range;
mod ou_half_life;
mod pain_index;
mod pair_spread_zscore;
mod pairwise_beta;
@@ -218,6 +223,8 @@ mod rocp;
mod rocr;
mod rocr100;
mod rogers_satchell;
mod rolling_correlation;
mod rolling_covariance;
mod roofing_filter;
mod rsi;
mod rvi;
@@ -237,6 +244,8 @@ mod smma;
mod sortino_ratio;
mod spearman_correlation;
mod spinning_top;
mod spread_bollinger_bands;
mod spread_hurst;
mod stalled_pattern;
mod standard_error;
mod standard_error_bands;
@@ -295,6 +304,7 @@ mod upside_gap_two_crows;
mod value_area;
mod value_at_risk;
mod variance;
mod variance_ratio;
mod vertical_horizontal_filter;
mod vidya;
mod volty_stop;
@@ -348,6 +358,7 @@ pub use awesome_oscillator_histogram::AwesomeOscillatorHistogram;
pub use balance_of_power::BalanceOfPower;
pub use belt_hold::BeltHold;
pub use beta::Beta;
pub use beta_neutral_spread::BetaNeutralSpread;
pub use bollinger::{BollingerBands, BollingerOutput};
pub use bollinger_bandwidth::BollingerBandwidth;
pub use breakaway::Breakaway;
@@ -382,6 +393,7 @@ pub use demand_index::DemandIndex;
pub use demark_pivots::{DemarkPivots, DemarkPivotsOutput};
pub use depth_slope::DepthSlope;
pub use detrended_std_dev::DetrendedStdDev;
pub use distance_ssd::DistanceSsd;
pub use doji::Doji;
pub use doji_star::DojiStar;
pub use donchian::{Donchian, DonchianOutput};
@@ -416,6 +428,7 @@ pub use funding_rate_zscore::FundingRateZScore;
pub use gain_loss_ratio::GainLossRatio;
pub use gap_side_by_side_white::GapSideBySideWhite;
pub use garman_klass::GarmanKlassVolatility;
pub use granger_causality::GrangerCausality;
pub use gravestone_doji::GravestoneDoji;
pub use hammer::Hammer;
pub use hanging_man::HangingMan;
@@ -445,6 +458,7 @@ pub use inverse_fisher_transform::InverseFisherTransform;
pub use inverted_hammer::InvertedHammer;
pub use jma::Jma;
pub use kagi_bars::{KagiBar, KagiBars};
pub use kalman_hedge_ratio::{KalmanHedgeRatio, KalmanHedgeRatioOutput};
pub use kama::Kama;
pub use kelly_criterion::KellyCriterion;
pub use keltner::{Keltner, KeltnerOutput};
@@ -502,6 +516,7 @@ pub use omega_ratio::OmegaRatio;
pub use on_neck::OnNeck;
pub use opening_marubozu::OpeningMarubozu;
pub use opening_range::{OpeningRange, OpeningRangeOutput};
pub use ou_half_life::OuHalfLife;
pub use pain_index::PainIndex;
pub use pair_spread_zscore::PairSpreadZScore;
pub use pairwise_beta::PairwiseBeta;
@@ -533,6 +548,8 @@ pub use rocp::Rocp;
pub use rocr::Rocr;
pub use rocr100::Rocr100;
pub use rogers_satchell::RogersSatchellVolatility;
pub use rolling_correlation::RollingCorrelation;
pub use rolling_covariance::RollingCovariance;
pub use roofing_filter::RoofingFilter;
pub use rsi::Rsi;
pub use rvi::Rvi;
@@ -552,6 +569,8 @@ pub use smma::Smma;
pub use sortino_ratio::SortinoRatio;
pub use spearman_correlation::SpearmanCorrelation;
pub use spinning_top::SpinningTop;
pub use spread_bollinger_bands::{SpreadBollingerBands, SpreadBollingerBandsOutput};
pub use spread_hurst::SpreadHurst;
pub use stalled_pattern::StalledPattern;
pub use standard_error::StandardError;
pub use standard_error_bands::{StandardErrorBands, StandardErrorBandsOutput};
@@ -610,6 +629,7 @@ pub use upside_gap_two_crows::UpsideGapTwoCrows;
pub use value_area::{ValueArea, ValueAreaOutput};
pub use value_at_risk::ValueAtRisk;
pub use variance::Variance;
pub use variance_ratio::VarianceRatio;
pub use vertical_horizontal_filter::VerticalHorizontalFilter;
pub use vidya::Vidya;
pub use volty_stop::VoltyStop;
@@ -848,6 +868,16 @@ pub const FAMILIES: &[(&str, &[&str])] = &[
"AvgPrice",
"LinRegIntercept",
"Tsf",
"RollingCorrelation",
"RollingCovariance",
"OuHalfLife",
"SpreadHurst",
"DistanceSsd",
"BetaNeutralSpread",
"VarianceRatio",
"GrangerCausality",
"KalmanHedgeRatio",
"SpreadBollingerBands",
],
),
(
@@ -1069,6 +1099,6 @@ mod family_tests {
// the actual indicator count is the early-warning signal that an
// indicator was added without being assigned a family.
let total: usize = FAMILIES.iter().map(|(_, ns)| ns.len()).sum();
assert_eq!(total, 315, "FAMILIES total drifted from indicator count");
assert_eq!(total, 325, "FAMILIES total drifted from indicator count");
}
}
@@ -0,0 +1,245 @@
//! OrnsteinUhlenbeck half-life of mean reversion for the spread of two series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Half-life of mean reversion of the spread `a b`, from an OrnsteinUhlenbeck
/// fit.
///
/// Each `update` takes one `(a, b)` price pair and forms the spread
/// `sₜ = aₜ bₜ`. Over the trailing window of `period` spreads the indicator
/// fits the discrete OrnsteinUhlenbeck (mean-reverting AR(1)) model by
/// ordinary least squares of the change on the level:
///
/// ```text
/// Δsₜ = λ · sₜ₋₁ + c + εₜ
/// half_life = ln(2) / λ (only when λ < 0)
/// ```
///
/// `λ` is the speed of mean reversion: a more negative `λ` pulls the spread back
/// to its mean faster. The **half-life** is the number of bars for a deviation
/// to decay by half — the single most useful number for sizing a pairs trade's
/// holding period and look-back. When the spread is not mean-reverting
/// (`λ ≥ 0`, a random walk or a trend) or the regression is degenerate (a flat
/// spread), the indicator returns `0`, meaning "no finite half-life".
///
/// Each `update` is `O(period)`: the OLS slope is recomputed from the window's
/// running geometry. Output is in bars and is always `≥ 0`.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, OuHalfLife};
///
/// let mut hl = OuHalfLife::new(40).unwrap();
/// let mut last = None;
/// for t in 0..120 {
/// let b = 100.0 + f64::from(t);
/// // `a` hugs `b` with a fast mean-reverting wobble ⇒ short half-life.
/// let a = b + 2.0 * (f64::from(t) * 0.9).sin();
/// last = hl.update((a, b));
/// }
/// let half_life = last.unwrap();
/// assert!(half_life > 0.0 && half_life < 40.0);
/// ```
#[derive(Debug, Clone)]
pub struct OuHalfLife {
period: usize,
window: VecDeque<f64>,
}
impl OuHalfLife {
/// Construct a new OrnsteinUhlenbeck half-life estimator.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 3` — the AR(1) regression
/// needs at least two observations (a slope and an intercept).
pub fn new(period: usize) -> Result<Self> {
if period < 3 {
return Err(Error::InvalidPeriod {
message: "OU half-life needs period >= 3",
});
}
Ok(Self {
period,
window: VecDeque::with_capacity(period),
})
}
/// Configured look-back window of spreads.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for OuHalfLife {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (a, b) = input;
if self.window.len() == self.period {
self.window.pop_front();
}
self.window.push_back(a - b);
if self.window.len() < self.period {
return None;
}
// OLS slope λ of Δsₜ on sₜ₋₁ over the window.
let spreads: Vec<f64> = self.window.iter().copied().collect();
let count = (spreads.len() - 1) as f64;
let mut sum_level = 0.0;
let mut sum_delta = 0.0;
let mut sum_ll = 0.0;
let mut sum_ld = 0.0;
for pair in spreads.windows(2) {
let level = pair[0];
let delta = pair[1] - pair[0];
sum_level += level;
sum_delta += delta;
sum_ll += level * level;
sum_ld += level * delta;
}
let mean_level = sum_level / count;
let mean_delta = sum_delta / count;
let var_level = sum_ll / count - mean_level * mean_level;
if var_level <= 0.0 {
// Flat spread: the regression has no defined slope.
return Some(0.0);
}
let cov = sum_ld / count - mean_level * mean_delta;
let lambda = cov / var_level;
if lambda >= 0.0 {
// Not mean-reverting (random walk or diverging): no finite half-life.
return Some(0.0);
}
Some(-std::f64::consts::LN_2 / lambda)
}
fn reset(&mut self) {
self.window.clear();
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"OuHalfLife"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
#[test]
fn rejects_period_below_three() {
assert!(OuHalfLife::new(2).is_err());
assert!(OuHalfLife::new(3).is_ok());
}
#[test]
fn accessors_and_metadata() {
let hl = OuHalfLife::new(30).unwrap();
assert_eq!(hl.period(), 30);
assert_eq!(hl.warmup_period(), 30);
assert_eq!(hl.name(), "OuHalfLife");
assert!(!hl.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut hl = OuHalfLife::new(4).unwrap();
assert_eq!(hl.update((1.0, 0.0)), None);
assert_eq!(hl.update((2.0, 0.0)), None);
assert_eq!(hl.update((3.0, 0.0)), None);
assert!(hl.update((4.0, 0.0)).is_some());
assert!(hl.is_ready());
}
#[test]
fn mean_reverting_spread_has_positive_half_life() {
// Fast sinusoidal spread around zero ⇒ strong mean reversion.
let pairs: Vec<(f64, f64)> = (0..120)
.map(|t| {
let b = 100.0 + f64::from(t);
let a = b + 2.0 * (f64::from(t) * 0.9).sin();
(a, b)
})
.collect();
let last = OuHalfLife::new(40)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last > 0.0 && last < 40.0, "half-life {last}");
}
#[test]
fn trending_spread_has_zero_half_life() {
// Spread = a b grows monotonically (λ ≥ 0) ⇒ no finite half-life.
let pairs: Vec<(f64, f64)> = (0..40)
.map(|t| (2.0 * f64::from(t), f64::from(t)))
.collect();
let last = OuHalfLife::new(20)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 0.0);
}
#[test]
fn flat_spread_returns_zero() {
// a b is constant ⇒ var(level) = 0 ⇒ undefined ⇒ 0.
let pairs: Vec<(f64, f64)> = (0..30)
.map(|t| (5.0 + f64::from(t), f64::from(t)))
.collect();
let last = OuHalfLife::new(10)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 0.0);
}
#[test]
fn reset_clears_state() {
let mut hl = OuHalfLife::new(5).unwrap();
for t in 0..10 {
hl.update((f64::from(t) + (f64::from(t) * 0.7).sin(), f64::from(t)));
}
assert!(hl.is_ready());
hl.reset();
assert!(!hl.is_ready());
assert_eq!(hl.update((1.0, 0.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..80)
.map(|t| {
let b = 50.0 + 0.5 * f64::from(t);
(b + (f64::from(t) * 0.6).sin(), b)
})
.collect();
let batch = OuHalfLife::new(25).unwrap().batch(&pairs);
let mut hl = OuHalfLife::new(25).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| hl.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,275 @@
//! Rolling Pearson correlation of the period-over-period *returns* of two series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Rolling correlation of the **returns** of two synchronised series.
///
/// Where [`crate::PearsonCorrelation`] correlates the raw *levels* `(x, y)`,
/// this indicator first differences each channel into a one-step return and
/// correlates those returns over the trailing window:
///
/// ```text
/// rxₜ = xₜ xₜ₋₁ ryₜ = yₜ − yₜ₋₁
/// corr = cov(rx, ry) / √(var(rx) · var(ry))
/// ```
///
/// Return correlation is the quantity that matters for hedging and portfolio
/// risk: two assets can trend together (high level correlation) while their
/// day-to-day moves are nearly independent (low return correlation). The output
/// is in `[1, +1]`; a flat return channel makes the ratio undefined and the
/// indicator reports `0` rather than `NaN`. The value is clamped to `[1, +1]`
/// to absorb tiny floating-point overshoots near the boundaries.
///
/// Each `update` is O(1): the five running sums (`Σrx`, `Σry`, `Σrx²`, `Σry²`,
/// `Σrxry`) are maintained as the window of returns slides. The first level in
/// each channel produces no return, so a `period`-pair correlation needs
/// `period + 1` updates of warmup.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, RollingCorrelation};
///
/// let mut rc = RollingCorrelation::new(10).unwrap();
/// let mut last = None;
/// for i in 0..40 {
/// // A varying path where y always moves with x ⇒ return correlation +1.
/// let x = (f64::from(i) * 0.5).sin() * 10.0;
/// last = rc.update((x, 2.0 * x));
/// }
/// assert!((last.unwrap() - 1.0).abs() < 1e-9);
/// ```
#[derive(Debug, Clone)]
pub struct RollingCorrelation {
period: usize,
prev: Option<(f64, f64)>,
window: VecDeque<(f64, f64)>,
sum_x: f64,
sum_y: f64,
sum_xx: f64,
sum_yy: f64,
sum_xy: f64,
}
impl RollingCorrelation {
/// Construct a new rolling return-correlation.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 2` — correlation is
/// undefined for fewer than two return pairs.
pub fn new(period: usize) -> Result<Self> {
if period < 2 {
return Err(Error::InvalidPeriod {
message: "rolling correlation needs period >= 2",
});
}
Ok(Self {
period,
prev: None,
window: VecDeque::with_capacity(period),
sum_x: 0.0,
sum_y: 0.0,
sum_xx: 0.0,
sum_yy: 0.0,
sum_xy: 0.0,
})
}
/// Configured window of returns.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for RollingCorrelation {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (x, y) = input;
let Some((px, py)) = self.prev else {
// First level in each channel: store it, no return yet.
self.prev = Some((x, y));
return None;
};
self.prev = Some((x, y));
let (rx, ry) = (x - px, y - py);
if self.window.len() == self.period {
let (ox, oy) = self.window.pop_front().expect("non-empty");
self.sum_x -= ox;
self.sum_y -= oy;
self.sum_xx -= ox * ox;
self.sum_yy -= oy * oy;
self.sum_xy -= ox * oy;
}
self.window.push_back((rx, ry));
self.sum_x += rx;
self.sum_y += ry;
self.sum_xx += rx * rx;
self.sum_yy += ry * ry;
self.sum_xy += rx * ry;
if self.window.len() < self.period {
return None;
}
let n = self.period as f64;
let mean_x = self.sum_x / n;
let mean_y = self.sum_y / n;
let var_x = (self.sum_xx / n - mean_x * mean_x).max(0.0);
let var_y = (self.sum_yy / n - mean_y * mean_y).max(0.0);
let cov = self.sum_xy / n - mean_x * mean_y;
let denom = (var_x * var_y).sqrt();
if denom == 0.0 {
// At least one return channel is flat: correlation is undefined.
return Some(0.0);
}
Some((cov / denom).clamp(-1.0, 1.0))
}
fn reset(&mut self) {
self.prev = None;
self.window.clear();
self.sum_x = 0.0;
self.sum_y = 0.0;
self.sum_xx = 0.0;
self.sum_yy = 0.0;
self.sum_xy = 0.0;
}
fn warmup_period(&self) -> usize {
self.period + 1
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"RollingCorrelation"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_period_below_two() {
assert!(RollingCorrelation::new(0).is_err());
assert!(RollingCorrelation::new(1).is_err());
assert!(RollingCorrelation::new(2).is_ok());
}
#[test]
fn accessors_and_metadata() {
let rc = RollingCorrelation::new(14).unwrap();
assert_eq!(rc.period(), 14);
assert_eq!(rc.warmup_period(), 15);
assert_eq!(rc.name(), "RollingCorrelation");
assert!(!rc.is_ready());
}
#[test]
fn warmup_needs_period_plus_one() {
let mut rc = RollingCorrelation::new(3).unwrap();
// First update only seeds the previous level ⇒ None.
assert_eq!(rc.update((1.0, 1.0)), None);
assert_eq!(rc.update((2.0, 3.0)), None); // 1 return
assert_eq!(rc.update((3.0, 5.0)), None); // 2 returns
assert!(rc.update((4.0, 7.0)).is_some()); // 3 returns ⇒ ready
assert!(rc.is_ready());
}
#[test]
fn comoving_returns_are_plus_one() {
// y always moves by 2x x's move ⇒ perfectly correlated returns.
let pairs: Vec<(f64, f64)> = (0..20)
.map(|i| {
let x = (f64::from(i) * 0.5).sin() * 10.0;
(x, 2.0 * x + 100.0)
})
.collect();
let last = RollingCorrelation::new(8)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 1.0, epsilon = 1e-9);
}
#[test]
fn opposing_returns_are_minus_one() {
let pairs: Vec<(f64, f64)> = (0..20)
.map(|i| {
let x = (f64::from(i) * 0.5).sin() * 10.0;
(x, -1.5 * x + 50.0)
})
.collect();
let last = RollingCorrelation::new(8)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, -1.0, epsilon = 1e-9);
}
#[test]
fn flat_return_channel_yields_zero() {
// y is constant ⇒ its returns are all zero ⇒ undefined ⇒ 0.
let pairs: Vec<(f64, f64)> = (0..20).map(|i| (f64::from(i), 7.0)).collect();
let last = RollingCorrelation::new(6)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.0, epsilon = 1e-12);
}
#[test]
fn output_in_range() {
let pairs: Vec<(f64, f64)> = (0..80)
.map(|i| {
let t = f64::from(i);
(100.0 + t.sin() * 5.0, 50.0 + (t * 0.3).cos() * 3.0)
})
.collect();
let mut rc = RollingCorrelation::new(20).unwrap();
for v in rc.batch(&pairs).into_iter().flatten() {
assert!((-1.0..=1.0).contains(&v));
}
}
#[test]
fn reset_clears_state() {
let mut rc = RollingCorrelation::new(4).unwrap();
rc.batch(&[(1.0, 2.0), (2.0, 4.0), (3.0, 6.0), (4.0, 8.0), (5.0, 10.0)]);
assert!(rc.is_ready());
rc.reset();
assert!(!rc.is_ready());
assert_eq!(rc.update((1.0, 1.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..60)
.map(|i| {
let t = f64::from(i);
(t.sin(), (t * 0.5).cos())
})
.collect();
let batch = RollingCorrelation::new(14).unwrap().batch(&pairs);
let mut rc = RollingCorrelation::new(14).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| rc.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,243 @@
//! Rolling covariance of the period-over-period *returns* of two series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Rolling covariance of the **returns** of two synchronised series.
///
/// Each `update` takes one `(x, y)` level pair, differences each channel into a
/// one-step return, and reports the population covariance of those returns over
/// the trailing window of `period` return pairs:
///
/// ```text
/// rxₜ = xₜ xₜ₋₁ ryₜ = yₜ − yₜ₋₁
/// cov = (1/n) · Σ rx·ry r̄x · r̄y
/// ```
///
/// Unlike [`crate::RollingCorrelation`] the result is **not** normalised to
/// `[1, 1]`: it carries the units of the two return streams multiplied
/// together, so it scales with volatility. It is the raw building block behind
/// correlation, beta and portfolio variance — positive when the two return
/// streams tend to move the same way, negative when they offset.
///
/// Each `update` is O(1): three running sums (`Σrx`, `Σry`, `Σrxry`) are
/// maintained as the window slides. The first level in each channel produces no
/// return, so a `period`-pair covariance needs `period + 1` updates of warmup.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, RollingCovariance};
///
/// let mut rc = RollingCovariance::new(5).unwrap();
/// let mut last = None;
/// for i in 0..20 {
/// let x = f64::from(i);
/// last = rc.update((x, 3.0 * x)); // y's return is 3× x's return
/// }
/// // cov(rx, ry) = cov(1, 3) over constant unit returns = 3 · var(rx) = 0
/// // for a constant return; use a varying path in practice. Here returns are
/// // constant (1 and 3) ⇒ covariance 0.
/// assert!(last.unwrap().abs() < 1e-9);
/// ```
#[derive(Debug, Clone)]
pub struct RollingCovariance {
period: usize,
prev: Option<(f64, f64)>,
window: VecDeque<(f64, f64)>,
sum_x: f64,
sum_y: f64,
sum_xy: f64,
}
impl RollingCovariance {
/// Construct a new rolling return-covariance.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 2` — covariance is
/// undefined for fewer than two return pairs.
pub fn new(period: usize) -> Result<Self> {
if period < 2 {
return Err(Error::InvalidPeriod {
message: "rolling covariance needs period >= 2",
});
}
Ok(Self {
period,
prev: None,
window: VecDeque::with_capacity(period),
sum_x: 0.0,
sum_y: 0.0,
sum_xy: 0.0,
})
}
/// Configured window of returns.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for RollingCovariance {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (x, y) = input;
let Some((px, py)) = self.prev else {
self.prev = Some((x, y));
return None;
};
self.prev = Some((x, y));
let (rx, ry) = (x - px, y - py);
if self.window.len() == self.period {
let (ox, oy) = self.window.pop_front().expect("non-empty");
self.sum_x -= ox;
self.sum_y -= oy;
self.sum_xy -= ox * oy;
}
self.window.push_back((rx, ry));
self.sum_x += rx;
self.sum_y += ry;
self.sum_xy += rx * ry;
if self.window.len() < self.period {
return None;
}
let n = self.period as f64;
let mean_x = self.sum_x / n;
let mean_y = self.sum_y / n;
Some(self.sum_xy / n - mean_x * mean_y)
}
fn reset(&mut self) {
self.prev = None;
self.window.clear();
self.sum_x = 0.0;
self.sum_y = 0.0;
self.sum_xy = 0.0;
}
fn warmup_period(&self) -> usize {
self.period + 1
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"RollingCovariance"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_period_below_two() {
assert!(RollingCovariance::new(0).is_err());
assert!(RollingCovariance::new(1).is_err());
assert!(RollingCovariance::new(2).is_ok());
}
#[test]
fn accessors_and_metadata() {
let rc = RollingCovariance::new(14).unwrap();
assert_eq!(rc.period(), 14);
assert_eq!(rc.warmup_period(), 15);
assert_eq!(rc.name(), "RollingCovariance");
assert!(!rc.is_ready());
}
#[test]
fn warmup_needs_period_plus_one() {
let mut rc = RollingCovariance::new(3).unwrap();
assert_eq!(rc.update((1.0, 1.0)), None);
assert_eq!(rc.update((2.0, 3.0)), None);
assert_eq!(rc.update((3.0, 5.0)), None);
assert!(rc.update((4.0, 7.0)).is_some());
assert!(rc.is_ready());
}
#[test]
fn hand_computed_value() {
// Levels x = 0,1,3,6,10 ⇒ returns 1,2,3,4; y = 2x ⇒ returns 2,4,6,8.
// With period = 3 the final window is rx = [2,3,4], ry = [4,6,8]:
// Σrx·ry/3 = 58/3, r̄x·r̄y = 3·6 = 18 ⇒ cov = 58/3 18 = 4/3.
let pairs = [
(0.0, 0.0),
(1.0, 2.0),
(3.0, 6.0),
(6.0, 12.0),
(10.0, 20.0),
];
let last = RollingCovariance::new(3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 4.0 / 3.0, epsilon = 1e-9);
}
#[test]
fn opposing_returns_give_negative_covariance() {
let pairs: Vec<(f64, f64)> = (0..30)
.map(|i| {
let x = (f64::from(i) * 0.4).sin() * 10.0;
(x, -x)
})
.collect();
let last = RollingCovariance::new(10)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last < 0.0, "cov {last}");
}
#[test]
fn flat_channel_gives_zero() {
let pairs: Vec<(f64, f64)> = (0..20).map(|i| (f64::from(i), 7.0)).collect();
let last = RollingCovariance::new(6)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.0, epsilon = 1e-12);
}
#[test]
fn reset_clears_state() {
let mut rc = RollingCovariance::new(4).unwrap();
rc.batch(&[(1.0, 2.0), (2.0, 4.0), (3.0, 1.0), (4.0, 9.0), (5.0, 2.0)]);
assert!(rc.is_ready());
rc.reset();
assert!(!rc.is_ready());
assert_eq!(rc.update((1.0, 1.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..60)
.map(|i| {
let t = f64::from(i);
(t.sin() * 4.0, (t * 0.5).cos() * 2.0)
})
.collect();
let batch = RollingCovariance::new(12).unwrap().batch(&pairs);
let mut rc = RollingCovariance::new(12).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| rc.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,278 @@
//! Bollinger bands on the spread of two series, for pairs mean-reversion trading.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Output of [`SpreadBollingerBands`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SpreadBollingerBandsOutput {
/// Middle band: the rolling mean of the spread.
pub middle: f64,
/// Upper band: `middle + num_std · σ`.
pub upper: f64,
/// Lower band: `middle num_std · σ`.
pub lower: f64,
/// `%b`: where the current spread sits across the band, `(s lower) /
/// (upper lower)`. `0` is the lower band, `1` the upper, `0.5` the middle.
/// Reported as `0.5` when the band has zero width (a flat spread).
pub percent_b: f64,
}
/// Bollinger bands on the spread `a b` of two series.
///
/// Each `update` takes one `(a, b)` price pair and forms the spread
/// `sₜ = aₜ bₜ`. Over the trailing window of `period` spreads it builds a
/// classic Bollinger envelope:
///
/// ```text
/// middle = mean(s) σ = stddev(s)
/// upper = middle + num_std · σ
/// lower = middle num_std · σ
/// %b = (s_now lower) / (upper lower)
/// ```
///
/// Applied to a spread rather than a price, the bands are a ready-made pairs
/// mean-reversion signal: the spread riding the **upper** band is stretched
/// rich (a short-the-spread setup), the **lower** band stretched cheap, and a
/// return to the **middle** is the exit. `%b` compresses the location into one
/// number for thresholding. The spread is the raw difference `a b`, so feed
/// already-comparable legs (e.g. a hedged pair, two yields, or log prices); pair
/// this with [`crate::BetaNeutralSpread`] when the legs need a hedge ratio first.
///
/// A flat spread yields a zero-width band; `%b` is then reported as the neutral
/// `0.5`. Each `update` is `O(1)`: the mean and variance come from two running
/// sums maintained as the window slides.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, SpreadBollingerBands};
///
/// let mut bb = SpreadBollingerBands::new(20, 2.0).unwrap();
/// let mut last = None;
/// for t in 0..60 {
/// let b = 100.0 + f64::from(t);
/// let a = b + 2.0 * (f64::from(t) * 0.5).sin();
/// last = bb.update((a, b));
/// }
/// let out = last.unwrap();
/// assert!(out.lower <= out.middle && out.middle <= out.upper);
/// ```
#[derive(Debug, Clone)]
pub struct SpreadBollingerBands {
period: usize,
num_std: f64,
window: VecDeque<f64>,
sum: f64,
sum_sq: f64,
}
impl SpreadBollingerBands {
/// Construct new spread Bollinger bands.
///
/// `period` is the look-back window; `num_std` is the band width in standard
/// deviations.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 2`, or
/// [`Error::InvalidParameter`] if `num_std` is not strictly positive (and
/// finite).
pub fn new(period: usize, num_std: f64) -> Result<Self> {
if period < 2 {
return Err(Error::InvalidPeriod {
message: "spread bollinger bands needs period >= 2",
});
}
if !num_std.is_finite() || num_std <= 0.0 {
return Err(Error::InvalidParameter {
message: "spread bollinger bands needs num_std > 0",
});
}
Ok(Self {
period,
num_std,
window: VecDeque::with_capacity(period),
sum: 0.0,
sum_sq: 0.0,
})
}
/// Configured look-back window.
pub const fn period(&self) -> usize {
self.period
}
/// Configured band width in standard deviations.
pub const fn num_std(&self) -> f64 {
self.num_std
}
}
impl Indicator for SpreadBollingerBands {
type Input = (f64, f64);
type Output = SpreadBollingerBandsOutput;
fn update(&mut self, input: (f64, f64)) -> Option<SpreadBollingerBandsOutput> {
let (a, b) = input;
let spread = a - b;
if self.window.len() == self.period {
let old = self.window.pop_front().expect("non-empty");
self.sum -= old;
self.sum_sq -= old * old;
}
self.window.push_back(spread);
self.sum += spread;
self.sum_sq += spread * spread;
if self.window.len() < self.period {
return None;
}
let n = self.period as f64;
let middle = self.sum / n;
let variance = (self.sum_sq / n - middle * middle).max(0.0);
let sigma = variance.sqrt();
let half_width = self.num_std * sigma;
let upper = middle + half_width;
let lower = middle - half_width;
let percent_b = if half_width == 0.0 {
0.5
} else {
(spread - lower) / (upper - lower)
};
Some(SpreadBollingerBandsOutput {
middle,
upper,
lower,
percent_b,
})
}
fn reset(&mut self) {
self.window.clear();
self.sum = 0.0;
self.sum_sq = 0.0;
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"SpreadBollingerBands"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_bad_parameters() {
assert!(SpreadBollingerBands::new(1, 2.0).is_err());
assert!(SpreadBollingerBands::new(20, 0.0).is_err());
assert!(SpreadBollingerBands::new(20, -1.0).is_err());
assert!(SpreadBollingerBands::new(20, f64::NAN).is_err());
assert!(SpreadBollingerBands::new(2, 2.0).is_ok());
}
#[test]
fn accessors_and_metadata() {
let bb = SpreadBollingerBands::new(20, 2.5).unwrap();
assert_eq!(bb.period(), 20);
assert_eq!(bb.num_std(), 2.5);
assert_eq!(bb.warmup_period(), 20);
assert_eq!(bb.name(), "SpreadBollingerBands");
assert!(!bb.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut bb = SpreadBollingerBands::new(3, 2.0).unwrap();
assert_eq!(bb.update((1.0, 0.0)), None);
assert_eq!(bb.update((2.0, 0.0)), None);
assert!(bb.update((3.0, 0.0)).is_some());
assert!(bb.is_ready());
}
#[test]
fn hand_computed_value() {
// Spreads 1,2,3,4 (b = 0), period 4, num_std 2:
// mean = 2.5, σ = √1.25, upper = 2.5 + 2√1.25, lower = 2.5 2√1.25,
// %b at s = 4 ⇒ 0.8354102.
let pairs = [(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)];
let out = SpreadBollingerBands::new(4, 2.0)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(out.middle, 2.5, epsilon = 1e-9);
assert_relative_eq!(out.upper, 4.736_067_977_499_79, epsilon = 1e-9);
assert_relative_eq!(out.lower, 0.263_932_022_500_21, epsilon = 1e-9);
assert_relative_eq!(out.percent_b, 0.835_410_196_624_97, epsilon = 1e-9);
}
#[test]
fn flat_spread_collapses_band() {
// a b constant ⇒ σ = 0 ⇒ upper = middle = lower, %b = 0.5.
let pairs: Vec<(f64, f64)> = (0..10)
.map(|t| (5.0 + f64::from(t), f64::from(t)))
.collect();
let out = SpreadBollingerBands::new(5, 2.0)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(out.upper, out.middle, epsilon = 1e-12);
assert_relative_eq!(out.lower, out.middle, epsilon = 1e-12);
assert_relative_eq!(out.percent_b, 0.5, epsilon = 1e-12);
}
#[test]
fn bands_are_ordered() {
let pairs: Vec<(f64, f64)> = (0..80)
.map(|t| {
let b = 100.0 + f64::from(t);
(b + 3.0 * (f64::from(t) * 0.4).sin(), b)
})
.collect();
let mut bb = SpreadBollingerBands::new(20, 2.0).unwrap();
for out in bb.batch(&pairs).into_iter().flatten() {
assert!(out.lower <= out.middle && out.middle <= out.upper);
}
}
#[test]
fn reset_clears_state() {
let mut bb = SpreadBollingerBands::new(4, 2.0).unwrap();
bb.batch(&[(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (5.0, 0.0)]);
assert!(bb.is_ready());
bb.reset();
assert!(!bb.is_ready());
assert_eq!(bb.update((1.0, 0.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..60)
.map(|t| {
let b = 30.0 + 0.7 * f64::from(t);
(b + (f64::from(t) * 0.4).sin() * 1.5, b)
})
.collect();
let batch = SpreadBollingerBands::new(15, 2.0).unwrap().batch(&pairs);
let mut bb = SpreadBollingerBands::new(15, 2.0).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| bb.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,271 @@
//! Hurst exponent of the spread of two series, for pairs-trading regime detection.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// Hurst exponent of the spread `a b` over a rolling window.
///
/// Each `update` takes one `(a, b)` price pair and forms the spread
/// `sₜ = aₜ bₜ`. Over the trailing window of `period` spreads the indicator
/// estimates the Hurst exponent `H` from how the variance of `τ`-lagged
/// differences grows with the lag `τ`:
///
/// ```text
/// V(τ) = mean_t (s_{t+τ} s_t)² ∝ τ^(2H)
/// H = slope of log V(τ) on log τ, divided by two
/// ```
///
/// `H` classifies the spread's regime:
///
/// * `H < 0.5` — **mean-reverting** (anti-persistent): the spread snaps back,
/// the regime pairs traders want.
/// * `H ≈ 0.5` — a **random walk**: no exploitable structure.
/// * `H > 0.5` — **trending** (persistent): the spread keeps diverging.
///
/// The fit uses lags `1..=period/4` (at least two). When the spread is flat —
/// every lagged difference is zero, so the log-regression has fewer than two
/// usable points — the indicator returns the neutral `0.5`. The output is
/// clamped to `[0, 1]`.
///
/// Each `update` is `O(period · period/4)`, bounded by the fixed window.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, SpreadHurst};
///
/// let mut h = SpreadHurst::new(60).unwrap();
/// let mut last = None;
/// for t in 0..200 {
/// let b = 100.0 + f64::from(t);
/// // A tight oscillating spread is anti-persistent ⇒ H < 0.5.
/// let a = b + 3.0 * (f64::from(t) * 0.8).sin();
/// last = h.update((a, b));
/// }
/// assert!(last.unwrap() < 0.5);
/// ```
#[derive(Debug, Clone)]
pub struct SpreadHurst {
period: usize,
max_lag: usize,
window: VecDeque<f64>,
}
impl SpreadHurst {
/// Construct a new spread Hurst estimator.
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `period < 8` — fewer than eight
/// spreads cannot support a two-lag loglog regression.
pub fn new(period: usize) -> Result<Self> {
if period < 8 {
return Err(Error::InvalidPeriod {
message: "spread Hurst needs period >= 8",
});
}
Ok(Self {
period,
max_lag: (period / 4).max(2),
window: VecDeque::with_capacity(period),
})
}
/// Configured look-back window of spreads.
pub const fn period(&self) -> usize {
self.period
}
}
impl Indicator for SpreadHurst {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (a, b) = input;
if self.window.len() == self.period {
self.window.pop_front();
}
self.window.push_back(a - b);
if self.window.len() < self.period {
return None;
}
let spreads: Vec<f64> = self.window.iter().copied().collect();
// Collect (log τ, log V(τ)) for every lag whose variance is positive.
let mut log_lag = Vec::with_capacity(self.max_lag);
let mut log_var = Vec::with_capacity(self.max_lag);
for lag in 1..=self.max_lag {
let mut sum_sq = 0.0;
let mut count = 0.0;
for pair in spreads.windows(lag + 1) {
let diff = pair[lag] - pair[0];
sum_sq += diff * diff;
count += 1.0;
}
let var = sum_sq / count;
if var > 0.0 {
log_lag.push((lag as f64).ln());
log_var.push(var.ln());
}
}
if log_lag.len() < 2 {
// Degenerate (flat) spread: report the random-walk midpoint.
return Some(0.5);
}
let n = log_lag.len() as f64;
let mean_lag = log_lag.iter().sum::<f64>() / n;
let mean_var = log_var.iter().sum::<f64>() / n;
let mut cov = 0.0;
let mut var_lag = 0.0;
for (lx, lv) in log_lag.iter().zip(&log_var) {
cov += (lx - mean_lag) * (lv - mean_var);
var_lag += (lx - mean_lag) * (lx - mean_lag);
}
// `log_lag` holds at least two *distinct* lag logarithms, so the lag
// variance is strictly positive — no degenerate-slope guard is needed.
let slope = cov / var_lag;
Some((slope / 2.0).clamp(0.0, 1.0))
}
fn reset(&mut self) {
self.window.clear();
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"SpreadHurst"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_period_below_eight() {
assert!(SpreadHurst::new(7).is_err());
assert!(SpreadHurst::new(8).is_ok());
}
#[test]
fn accessors_and_metadata() {
let h = SpreadHurst::new(40).unwrap();
assert_eq!(h.period(), 40);
assert_eq!(h.warmup_period(), 40);
assert_eq!(h.name(), "SpreadHurst");
assert!(!h.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut h = SpreadHurst::new(8).unwrap();
for t in 0..7 {
assert_eq!(h.update((f64::from(t), 0.0)), None);
}
assert!(h.update((7.0, 0.0)).is_some());
assert!(h.is_ready());
}
#[test]
fn oscillating_spread_is_anti_persistent() {
let pairs: Vec<(f64, f64)> = (0..200)
.map(|t| {
let b = 100.0 + f64::from(t);
(b + 3.0 * (f64::from(t) * 0.8).sin(), b)
})
.collect();
let last = SpreadHurst::new(60)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last < 0.5, "H {last}");
}
#[test]
fn linear_trend_spread_is_persistent() {
// Spread = a b = t ⇒ τ-lagged differences are all τ ⇒ V(τ) = τ² ⇒ H = 1.
let pairs: Vec<(f64, f64)> = (0..40)
.map(|t| (2.0 * f64::from(t), f64::from(t)))
.collect();
let last = SpreadHurst::new(20)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 1.0, epsilon = 1e-9);
}
#[test]
fn flat_spread_returns_midpoint() {
// a b is constant ⇒ all lagged differences zero ⇒ neutral 0.5.
let pairs: Vec<(f64, f64)> = (0..30)
.map(|t| (5.0 + f64::from(t), f64::from(t)))
.collect();
let last = SpreadHurst::new(16)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.5, epsilon = 1e-12);
}
#[test]
fn output_in_unit_range() {
let pairs: Vec<(f64, f64)> = (0..150)
.map(|t| {
let b = 50.0 + 0.3 * f64::from(t);
(
b + (f64::from(t) * 0.5).sin() * 2.0 + (f64::from(t) * 0.13).cos(),
b,
)
})
.collect();
let mut h = SpreadHurst::new(48).unwrap();
for v in h.batch(&pairs).into_iter().flatten() {
assert!((0.0..=1.0).contains(&v));
}
}
#[test]
fn reset_clears_state() {
let mut h = SpreadHurst::new(8).unwrap();
for t in 0..12 {
h.update((f64::from(t) + (f64::from(t) * 0.7).sin(), f64::from(t)));
}
assert!(h.is_ready());
h.reset();
assert!(!h.is_ready());
assert_eq!(h.update((1.0, 0.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..100)
.map(|t| {
let b = 30.0 + 0.7 * f64::from(t);
(b + (f64::from(t) * 0.4).sin() * 1.5, b)
})
.collect();
let batch = SpreadHurst::new(32).unwrap().batch(&pairs);
let mut h = SpreadHurst::new(32).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| h.update(*p)).collect();
assert_eq!(batch, streamed);
}
}
@@ -0,0 +1,267 @@
//! LoMacKinlay variance-ratio test on the spread of two series.
use std::collections::VecDeque;
use crate::error::{Error, Result};
use crate::traits::Indicator;
/// LoMacKinlay variance ratio of the spread `a b` at horizon `q`.
///
/// Each `update` takes one `(a, b)` price pair and forms the spread
/// `sₜ = aₜ bₜ`. Over the trailing window of `period` spreads the indicator
/// compares the variance of `q`-step changes against `q` times the variance of
/// one-step changes:
///
/// ```text
/// rₜ = sₜ sₜ₋₁ (one-step changes)
/// VR(q) = Var(Σ of q consecutive r) / (q · Var(r))
/// ```
///
/// Under a random walk the variance of returns grows linearly with the horizon,
/// so `VR(q) = 1`. Departures reveal autocorrelation structure:
///
/// * `VR(q) < 1` — **mean reversion** (negatively autocorrelated changes): the
/// spread's moves partly cancel, the regime pairs traders exploit.
/// * `VR(q) ≈ 1` — a **random walk**: no exploitable structure.
/// * `VR(q) > 1` — **momentum / trending** (positively autocorrelated changes).
///
/// The estimator uses overlapping `q`-step windows. When the one-step changes
/// have zero variance (a flat spread) the ratio is undefined and the indicator
/// returns the null value `1`. The output is always `≥ 0`.
///
/// Each `update` is `O(period)`, bounded by the fixed window.
///
/// # Example
///
/// ```
/// use wickra_core::{Indicator, VarianceRatio};
///
/// let mut vr = VarianceRatio::new(60, 2).unwrap();
/// let mut last = None;
/// for t in 0..200 {
/// let b = 100.0 + f64::from(t);
/// // A fast, choppy spread mean-reverts (negatively autocorrelated
/// // changes) ⇒ VR(2) < 1.
/// let a = b + 2.0 * (f64::from(t) * 2.5).sin();
/// last = vr.update((a, b));
/// }
/// assert!(last.unwrap() < 1.0);
/// ```
#[derive(Debug, Clone)]
pub struct VarianceRatio {
period: usize,
q: usize,
window: VecDeque<f64>,
}
impl VarianceRatio {
/// Construct a new variance-ratio test.
///
/// `period` is the look-back window of spreads; `q` is the aggregation
/// horizon (number of one-step changes summed per long-horizon change).
///
/// # Errors
/// Returns [`Error::InvalidPeriod`] if `q < 2` or if `period < q + 2`
/// (which would leave fewer than two long-horizon observations).
pub fn new(period: usize, q: usize) -> Result<Self> {
if q < 2 {
return Err(Error::InvalidPeriod {
message: "variance ratio needs q >= 2",
});
}
if period < q + 2 {
return Err(Error::InvalidPeriod {
message: "variance ratio needs period >= q + 2",
});
}
Ok(Self {
period,
q,
window: VecDeque::with_capacity(period),
})
}
/// Configured look-back window of spreads.
pub const fn period(&self) -> usize {
self.period
}
/// Configured aggregation horizon `q`.
pub const fn q(&self) -> usize {
self.q
}
}
impl Indicator for VarianceRatio {
type Input = (f64, f64);
type Output = f64;
fn update(&mut self, input: (f64, f64)) -> Option<f64> {
let (a, b) = input;
if self.window.len() == self.period {
self.window.pop_front();
}
self.window.push_back(a - b);
if self.window.len() < self.period {
return None;
}
let spreads: Vec<f64> = self.window.iter().copied().collect();
// One-step changes.
let returns: Vec<f64> = spreads.windows(2).map(|w| w[1] - w[0]).collect();
let m = returns.len() as f64;
let mean = returns.iter().sum::<f64>() / m;
let var_one = returns.iter().map(|r| (r - mean) * (r - mean)).sum::<f64>() / m;
if var_one <= 0.0 {
// Flat spread: the random-walk null value.
return Some(1.0);
}
// Overlapping q-step changes; their mean is q·mean by construction.
let q_mean = self.q as f64 * mean;
let long: Vec<f64> = returns.windows(self.q).map(|w| w.iter().sum()).collect();
let count = long.len() as f64;
let var_q = long
.iter()
.map(|y| (y - q_mean) * (y - q_mean))
.sum::<f64>()
/ count;
Some(var_q / (self.q as f64 * var_one))
}
fn reset(&mut self) {
self.window.clear();
}
fn warmup_period(&self) -> usize {
self.period
}
fn is_ready(&self) -> bool {
self.window.len() == self.period
}
fn name(&self) -> &'static str {
"VarianceRatio"
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::traits::BatchExt;
use approx::assert_relative_eq;
#[test]
fn rejects_bad_parameters() {
assert!(VarianceRatio::new(10, 1).is_err()); // q must be >= 2
assert!(VarianceRatio::new(3, 2).is_err()); // period must be >= q + 2
assert!(VarianceRatio::new(4, 2).is_ok());
}
#[test]
fn accessors_and_metadata() {
let vr = VarianceRatio::new(60, 4).unwrap();
assert_eq!(vr.period(), 60);
assert_eq!(vr.q(), 4);
assert_eq!(vr.warmup_period(), 60);
assert_eq!(vr.name(), "VarianceRatio");
assert!(!vr.is_ready());
}
#[test]
fn warmup_returns_none() {
let mut vr = VarianceRatio::new(4, 2).unwrap();
assert_eq!(vr.update((1.0, 0.0)), None);
assert_eq!(vr.update((2.0, 0.0)), None);
assert_eq!(vr.update((3.0, 0.0)), None);
assert!(vr.update((4.0, 0.0)).is_some());
assert!(vr.is_ready());
}
#[test]
fn alternating_changes_give_zero_ratio() {
// Spreads 0,2,1,3,2 ⇒ changes 2,-1,2,-1; q = 2 overlapping sums are all
// 1 (constant) ⇒ Var(q) = 0 ⇒ VR = 0 (perfect mean reversion).
let pairs = [(0.0, 0.0), (2.0, 0.0), (1.0, 0.0), (3.0, 0.0), (2.0, 0.0)];
let last = VarianceRatio::new(5, 2)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_relative_eq!(last, 0.0, epsilon = 1e-12);
}
#[test]
fn oscillating_spread_is_below_one() {
let pairs: Vec<(f64, f64)> = (0..200)
.map(|t| {
let b = 100.0 + f64::from(t);
(b + 2.0 * (f64::from(t) * 2.5).sin(), b)
})
.collect();
let last = VarianceRatio::new(60, 2)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert!(last < 1.0, "VR {last}");
}
#[test]
fn flat_spread_returns_one() {
let pairs: Vec<(f64, f64)> = (0..30)
.map(|t| (5.0 + f64::from(t), f64::from(t)))
.collect();
let last = VarianceRatio::new(10, 3)
.unwrap()
.batch(&pairs)
.into_iter()
.flatten()
.last()
.unwrap();
assert_eq!(last, 1.0);
}
#[test]
fn output_non_negative() {
let pairs: Vec<(f64, f64)> = (0..150)
.map(|t| {
let b = 50.0 + 0.3 * f64::from(t);
(b + (f64::from(t) * 0.5).sin() * 2.0, b)
})
.collect();
let mut vr = VarianceRatio::new(40, 4).unwrap();
for v in vr.batch(&pairs).into_iter().flatten() {
assert!(v >= 0.0, "VR {v}");
}
}
#[test]
fn reset_clears_state() {
let mut vr = VarianceRatio::new(6, 2).unwrap();
for t in 0..12 {
vr.update((f64::from(t) + (f64::from(t) * 0.7).sin(), f64::from(t)));
}
assert!(vr.is_ready());
vr.reset();
assert!(!vr.is_ready());
assert_eq!(vr.update((1.0, 0.0)), None);
}
#[test]
fn batch_equals_streaming() {
let pairs: Vec<(f64, f64)> = (0..100)
.map(|t| {
let b = 30.0 + 0.7 * f64::from(t);
(b + (f64::from(t) * 0.4).sin() * 1.5, b)
})
.collect();
let batch = VarianceRatio::new(32, 3).unwrap().batch(&pairs);
let mut vr = VarianceRatio::new(32, 3).unwrap();
let streamed: Vec<_> = pairs.iter().map(|p| vr.update(*p)).collect();
assert_eq!(batch, streamed);
}
}