qualia_core_db/specialized_libs/computational_economics/
risk.rs

1//! Risk metrics over supplied return and scenario data.
2//!
3//! These routines do not fetch market data or infer missing histories. Callers
4//! provide the return series or shocks and the output buffers.
5
6#[derive(Debug, Clone, Copy, PartialEq, Eq)]
7pub enum RiskError {
8    InvalidInput,
9    OutputBufferTooSmall,
10}
11
12fn valid_return(x: f64) -> bool {
13    x.is_finite() && x > -1.0
14}
15
16#[cfg(not(test))]
17fn normal_quantile(probability: f64) -> Result<f64, RiskError> {
18    if !(0.0..1.0).contains(&probability) {
19        return Err(RiskError::InvalidInput);
20    }
21    Ok(crate::solvers::statistics::distributions::normal::standard_quantile(probability))
22}
23
24#[cfg(test)]
25fn normal_quantile(probability: f64) -> Result<f64, RiskError> {
26    if !(0.0..1.0).contains(&probability) {
27        return Err(RiskError::InvalidInput);
28    }
29
30    // Peter John Acklam's inverse-normal approximation coefficients.
31    const A: [f64; 6] = [
32        -3.969_683_028_665_376e1,
33        2.209_460_984_245_205e2,
34        -2.759_285_104_469_687e2,
35        1.383_577_518_672_69e2,
36        -3.066_479_806_614_716e1,
37        2.506_628_277_459_239,
38    ];
39    const B: [f64; 5] = [
40        -5.447_609_879_822_406e1,
41        1.615_858_368_580_409e2,
42        -1.556_989_798_598_866e2,
43        6.680_131_188_771_972e1,
44        -1.328_068_155_288_572e1,
45    ];
46    const C: [f64; 6] = [
47        -7.784_894_002_430_293e-3,
48        -3.223_964_580_411_365e-1,
49        -2.400_758_277_161_838,
50        -2.549_732_539_343_734,
51        4.374_664_141_464_968,
52        2.938_163_982_698_783,
53    ];
54    const D: [f64; 4] = [
55        7.784_695_709_041_462e-3,
56        3.224_671_290_700_398e-1,
57        2.445_134_137_142_996,
58        3.754_408_661_907_416,
59    ];
60
61    let plow = 0.02425;
62    let phigh = 1.0 - plow;
63    if probability < plow {
64        let q = (-2.0 * probability.ln()).sqrt();
65        let numerator = ((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5];
66        let denominator = (((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0;
67        Ok(numerator / denominator)
68    } else if probability <= phigh {
69        let q = probability - 0.5;
70        let r = q * q;
71        let numerator = (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q;
72        let denominator = ((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0;
73        Ok(numerator / denominator)
74    } else {
75        let q = (-2.0 * (1.0 - probability).ln()).sqrt();
76        let numerator = ((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5];
77        let denominator = (((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0;
78        Ok(-(numerator / denominator))
79    }
80}
81
82/// Copy and sort returns ascending into caller-owned scratch.
83pub fn sorted_returns_into(returns: &[f64], scratch: &mut [f64]) -> Result<usize, RiskError> {
84    if returns.is_empty() || scratch.len() < returns.len() {
85        return Err(RiskError::OutputBufferTooSmall);
86    }
87    for (idx, value) in returns.iter().enumerate() {
88        if !valid_return(*value) {
89            return Err(RiskError::InvalidInput);
90        }
91        scratch[idx] = *value;
92    }
93    scratch[..returns.len()].sort_by(|a, b| a.total_cmp(b));
94    Ok(returns.len())
95}
96
97/// Historical VaR as a positive loss fraction at `confidence`.
98pub fn historical_var(
99    returns: &[f64],
100    confidence: f64,
101    scratch: &mut [f64],
102) -> Result<f64, RiskError> {
103    if !(0.0..1.0).contains(&confidence) {
104        return Err(RiskError::InvalidInput);
105    }
106    let n = sorted_returns_into(returns, scratch)?;
107    let tail_probability = 1.0 - confidence;
108    let idx = ((tail_probability * n as f64).ceil() as usize).saturating_sub(1);
109    Ok((-scratch[idx]).max(0.0))
110}
111
112/// Historical expected shortfall/CVaR as a positive average tail loss.
113pub fn historical_cvar(
114    returns: &[f64],
115    confidence: f64,
116    scratch: &mut [f64],
117) -> Result<f64, RiskError> {
118    if !(0.0..1.0).contains(&confidence) {
119        return Err(RiskError::InvalidInput);
120    }
121    let n = sorted_returns_into(returns, scratch)?;
122    let tail_probability = 1.0 - confidence;
123    let tail_count = ((tail_probability * n as f64).ceil() as usize).max(1);
124    let mut loss = 0.0;
125    for value in scratch.iter().take(tail_count) {
126        loss += (-*value).max(0.0);
127    }
128    Ok(loss / tail_count as f64)
129}
130
131/// Parametric Gaussian VaR as a positive loss fraction.
132pub fn gaussian_var(mean: f64, std_dev: f64, confidence: f64) -> Result<f64, RiskError> {
133    if !mean.is_finite()
134        || !std_dev.is_finite()
135        || std_dev < 0.0
136        || !(0.0..1.0).contains(&confidence)
137    {
138        return Err(RiskError::InvalidInput);
139    }
140    let z_left_tail = normal_quantile(1.0 - confidence)?;
141    Ok((-(mean + z_left_tail * std_dev)).max(0.0))
142}
143
144/// Apply asset shocks to portfolio weights and return scenario portfolio loss.
145pub fn scenario_loss(weights: &[f64], shocks: &[f64]) -> Result<f64, RiskError> {
146    if weights.is_empty() || weights.len() != shocks.len() {
147        return Err(RiskError::InvalidInput);
148    }
149    let mut scenario_return = 0.0;
150    for idx in 0..weights.len() {
151        if !weights[idx].is_finite() || !valid_return(shocks[idx]) {
152            return Err(RiskError::InvalidInput);
153        }
154        scenario_return += weights[idx] * shocks[idx];
155    }
156    Ok((-scenario_return).max(0.0))
157}
158
159/// Scenario losses for row-major scenario shock matrix.
160pub fn scenario_losses_into(
161    weights: &[f64],
162    scenario_shocks: &[f64],
163    scenario_count: usize,
164    asset_count: usize,
165    out: &mut [f64],
166) -> Result<usize, RiskError> {
167    if scenario_count == 0
168        || asset_count == 0
169        || weights.len() != asset_count
170        || scenario_shocks.len() != scenario_count * asset_count
171    {
172        return Err(RiskError::InvalidInput);
173    }
174    if out.len() < scenario_count {
175        return Err(RiskError::OutputBufferTooSmall);
176    }
177
178    for scenario in 0..scenario_count {
179        let start = scenario * asset_count;
180        out[scenario] = scenario_loss(weights, &scenario_shocks[start..start + asset_count])?;
181    }
182    Ok(scenario_count)
183}
184
185#[cfg(test)]
186mod tests {
187    use super::*;
188
189    #[test]
190    fn historical_var_uses_left_tail_loss() {
191        let returns = [-0.10, -0.05, 0.0, 0.02, 0.03];
192        let mut scratch = [0.0; 5];
193        let var = historical_var(&returns, 0.80, &mut scratch).unwrap();
194        assert!((var - 0.10).abs() < 1e-12);
195    }
196
197    #[test]
198    fn historical_cvar_averages_tail_losses() {
199        let returns = [-0.10, -0.05, 0.0, 0.02, 0.03];
200        let mut scratch = [0.0; 5];
201        let cvar = historical_cvar(&returns, 0.60, &mut scratch).unwrap();
202        assert!((cvar - 0.075).abs() < 1e-12);
203    }
204
205    #[test]
206    fn gaussian_var_is_positive_for_left_tail() {
207        let var = gaussian_var(0.0, 0.02, 0.95).unwrap();
208        assert!(var > 0.032 && var < 0.034);
209    }
210
211    #[test]
212    fn scenario_loss_is_weighted_negative_return() {
213        let weights = [0.25, 0.75];
214        let shocks = [-0.10, -0.20];
215        let loss = scenario_loss(&weights, &shocks).unwrap();
216        assert!((loss - 0.175).abs() < 1e-12);
217    }
218
219    #[test]
220    fn scenario_losses_use_row_major_shocks() {
221        let weights = [0.5, 0.5];
222        let shocks = [-0.10, 0.0, 0.0, -0.20];
223        let mut out = [0.0; 2];
224        scenario_losses_into(&weights, &shocks, 2, 2, &mut out).unwrap();
225        assert!((out[0] - 0.05).abs() < 1e-12);
226        assert!((out[1] - 0.10).abs() < 1e-12);
227    }
228
229    #[test]
230    fn invalid_confidence_is_rejected() {
231        let returns = [0.0, 0.1];
232        let mut scratch = [0.0; 2];
233        assert_eq!(
234            historical_var(&returns, 1.0, &mut scratch),
235            Err(RiskError::InvalidInput)
236        );
237    }
238
239    #[test]
240    fn scratch_too_small_is_rejected() {
241        let returns = [0.0, 0.1];
242        let mut scratch = [0.0; 1];
243        assert_eq!(
244            sorted_returns_into(&returns, &mut scratch),
245            Err(RiskError::OutputBufferTooSmall)
246        );
247    }
248}
qualia_core_db/specialized_libs/computational_economics/risk.rs

qualia_core_db/specialized_libs/computational_economics/
risk.rs
