quantwave_core/regimes/hmm.rs
1//! Hidden Markov Models (Hamilton 1989)
2//!
3//! Source: Hamilton, J. D. (1989).
4//! "A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle."
5//! Econometrica, 57(2), 357-384.
6//!
7//! Implementation of a regime-switching Hidden Markov Model with Gaussian emissions.
8//! Includes the Viterbi algorithm for online state decoding and placeholder for Baum-Welch training.
9
10use crate::regimes::MarketRegime;
11use crate::traits::Next;
12use serde::{Deserialize, Serialize};
13
14/// A Hidden Markov Model for regime detection.
15#[derive(Debug, Clone, Serialize, Deserialize)]
16pub struct HMM {
17 n_states: usize,
18 /// Transition probability matrix [from_state][to_state]
19 a: Vec<Vec<f64>>,
20 /// Emission means for each state
21 means: Vec<f64>,
22 /// Emission standard deviations for each state
23 stds: Vec<f64>,
24 /// Initial state probabilities
25 pi: Vec<f64>,
26 /// Last Viterbi log-probabilities for each state
27 last_delta: Vec<f64>,
28 initialized: bool,
29}
30
31impl HMM {
32 /// Creates a new HMM with pre-defined parameters.
33 pub fn new(a: Vec<Vec<f64>>, means: Vec<f64>, stds: Vec<f64>, pi: Vec<f64>) -> Self {
34 let n_states = a.len();
35 Self {
36 n_states,
37 a,
38 means,
39 stds,
40 pi,
41 last_delta: vec![0.0; n_states],
42 initialized: false,
43 }
44 }
45
46 /// Default 2-state HMM (Bull/Bear).
47 ///
48 /// **Input contract: feed daily returns, not price.** The Gaussian emission
49 /// parameters below are hardcoded to the scale of daily fractional returns
50 /// (means near `0.001`/`-0.002`, stds near `0.01`/`0.02`). Feeding a raw price
51 /// series (e.g. `100.0`) makes `gaussian_pdf` underflow to `0.0` for *every*
52 /// state, which drives every `next()` call into a `-inf`/`-inf` tie that always
53 /// resolves to state `0` (Bull) — i.e. a silently constant output. The
54 /// `hmm_bull_bear` plugin wrapper (quantwave-py) calls [`Self::is_degenerate`]
55 /// after every [`Next::next`] to turn that case into an explicit error
56 /// instead of a silent wrong answer; callers driving the model directly
57 /// should do the same.
58 ///
59 /// **State labels are `{1, 2}`, not `{0, 1}`** at the plugin/FFI boundary —
60 /// `0` is reserved for "no regime yet" (`MarketRegime::Steady`), Bull maps to
61 /// `1`, Bear to `2`. See `hmm_bull_bear` (quantwave-py) and `MarketRegime`.
62 pub fn bull_bear() -> Self {
63 Self::new(
64 vec![
65 vec![0.95, 0.05], // Bull -> Bull, Bull -> Bear
66 vec![0.10, 0.90], // Bear -> Bull, Bear -> Bear
67 ],
68 vec![0.001, -0.002], // Daily returns: 0.1% for Bull, -0.2% for Bear
69 vec![0.01, 0.02], // Volatility: 1% for Bull, 2% for Bear
70 vec![0.5, 0.5],
71 )
72 }
73
74 /// Calculate Gaussian PDF: P(x | mu, sigma)
75 fn gaussian_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
76 let variance = sigma * sigma;
77 let denom = (2.0 * std::f64::consts::PI * variance).sqrt();
78 let exponent = -((x - mu).powi(2)) / (2.0 * variance);
79 exponent.exp() / denom
80 }
81
82 /// Returns `true` when the most recent [`Next::next`] update produced a
83 /// `-inf` log-likelihood for *every* state, i.e. `gaussian_pdf` underflowed to
84 /// `0.0` for all states. This is the degenerate-input signature described on
85 /// [`Self::bull_bear`]: it means the observation is off the scale the
86 /// emission parameters were fit for (classically, raw price fed where daily
87 /// returns were expected), and `best_state` would otherwise silently stick at
88 /// index `0` forever. Callers that need to fail loudly instead of emitting a
89 /// constant regime should check this after every `next()` call.
90 pub fn is_degenerate(&self) -> bool {
91 self.initialized
92 && self
93 .last_delta
94 .iter()
95 .all(|d| d.is_infinite() && d.is_sign_negative())
96 }
97
98 /// Normalized state probabilities from Viterbi log-deltas (for soft regime features).
99 pub fn state_probabilities(&self) -> Vec<f64> {
100 if !self.initialized {
101 return self.pi.clone();
102 }
103 let max_log = self
104 .last_delta
105 .iter()
106 .cloned()
107 .fold(f64::NEG_INFINITY, f64::max);
108 let mut probs: Vec<f64> = self
109 .last_delta
110 .iter()
111 .map(|&d| (d - max_log).exp())
112 .collect();
113 let sum: f64 = probs.iter().sum();
114 if sum > 0.0 {
115 for p in &mut probs {
116 *p /= sum;
117 }
118 }
119 probs
120 }
121}
122
123impl Next<f64> for HMM {
124 type Output = MarketRegime;
125
126 fn next(&mut self, x: f64) -> Self::Output {
127 let mut next_delta = vec![0.0; self.n_states];
128 let mut best_state = 0;
129 let mut max_prob = -f64::INFINITY;
130
131 if !self.initialized {
132 for i in 0..self.n_states {
133 let emission = Self::gaussian_pdf(x, self.means[i], self.stds[i]);
134 next_delta[i] = (self.pi[i] * emission).ln();
135 if next_delta[i] > max_prob {
136 max_prob = next_delta[i];
137 best_state = i;
138 }
139 }
140 self.initialized = true;
141 } else {
142 for j in 0..self.n_states {
143 let mut max_prev = -f64::INFINITY;
144 for i in 0..self.n_states {
145 // delta_j(t) = max_i [ delta_i(t-1) + ln(A_ij) ] + ln(P(x|j))
146 let prob = self.last_delta[i] + self.a[i][j].ln();
147 if prob > max_prev {
148 max_prev = prob;
149 }
150 }
151 let emission = Self::gaussian_pdf(x, self.means[j], self.stds[j]);
152 next_delta[j] = max_prev + emission.ln();
153
154 if next_delta[j] > max_prob {
155 max_prob = next_delta[j];
156 best_state = j;
157 }
158 }
159 }
160
161 self.last_delta = next_delta;
162
163 match best_state {
164 0 => MarketRegime::Bull,
165 1 => MarketRegime::Bear,
166 _ => MarketRegime::Cluster(best_state as u8),
167 }
168 }
169}