prospicio_prob/
distortion.rs

1//! Distortion risk measures: `ρ(X) = ∫ g(S(x)) dx` for a concave
2//! distortion `g` of the survival function, with `g(0) = 0`, `g(1) = 1`.
3//!
4//! On a discrete distribution the integral is a weighted sum of its
5//! values: the value `x_k` gets weight `g(P(X >= x_k)) - g(P(X > x_k))`.
6//! The weights depend only on ranks, which is what lets capital
7//! allocation reuse them on a joint distribution (see
8//! `docs/design/risk.md`).
9
10use prospicio_core::{Error, Result};
11use prospicio_math::special::{norm_cdf, norm_quantile};
12
13use crate::distribution::check_probability;
14
15/// A distortion of the survival function. Every variant is concave, so
16/// every measure here is coherent.
17///
18/// | Variant | `g(s)` | Parameter |
19/// |---|---|---|
20/// | `Tvar(p)` | `min(s / (1 - p), 1)` | `p` in `[0, 1]` |
21/// | `Wang(λ)` | `Φ(Φ⁻¹(s) + λ)` | `λ >= 0` |
22/// | `ProportionalHazard(ρ)` | `s^ρ` | `ρ` in `(0, 1]` |
23/// | `DualPower(β)` | `1 - (1 - s)^β` | `β >= 1` |
24/// | `Exponential(k)` | `(1 - e^(-k s)) / (1 - e^(-k))` | `k > 0` |
25///
26/// Each has a parameter value that gives the mean (`Tvar(0)`, `Wang(0)`,
27/// `ProportionalHazard(1)`, `DualPower(1)`, and `Exponential(k)` as
28/// `k → 0`).
29///
30/// A concave distortion is a spectral risk measure,
31/// `ρ = ∫₀¹ φ(u) VaR_u du` with the non-decreasing risk-aversion spectrum
32/// `φ(u) = g'(1 - u)`. `Exponential(k)` is the spectral measure with
33/// exponential risk aversion, `φ(u) = k e^(-k(1-u)) / (1 - e^(-k))`
34/// (Acerbi, 2002; Dowd, Cotter and Sorwar, 2008); `Tvar(p)` is the one
35/// whose spectrum is flat above `p`.
36///
37/// ```
38/// use prospicio_prob::Distortion;
39///
40/// let x = [1.0, 2.0, 3.0, 4.0];
41/// // TVaR at 50%: the mean of the top half.
42/// assert_eq!(Distortion::tvar(0.5).unwrap().apply_sorted(&x), 3.5);
43/// // Wang with λ = 0 is the mean.
44/// assert!((Distortion::wang(0.0).unwrap().apply_sorted(&x) - 2.5).abs() < 1e-15);
45/// ```
46#[derive(Debug, Clone, Copy, PartialEq)]
47pub enum Distortion {
48    Tvar(f64),
49    Wang(f64),
50    ProportionalHazard(f64),
51    DualPower(f64),
52    Exponential(f64),
53}
54
55impl Distortion {
56    /// Tail value at risk at level `p`.
57    pub fn tvar(p: f64) -> Result<Self> {
58        check_probability(p)?;
59        Ok(Self::Tvar(p))
60    }
61
62    /// Wang transform with market price of risk `lambda >= 0`.
63    pub fn wang(lambda: f64) -> Result<Self> {
64        if !lambda.is_finite() || lambda < 0.0 {
65            return Err(invalid("lambda", lambda, "must be finite and non-negative"));
66        }
67        Ok(Self::Wang(lambda))
68    }
69
70    /// Proportional hazard transform with `rho` in `(0, 1]`.
71    pub fn proportional_hazard(rho: f64) -> Result<Self> {
72        if !(rho > 0.0 && rho <= 1.0) {
73            return Err(invalid("rho", rho, "must be in (0, 1]"));
74        }
75        Ok(Self::ProportionalHazard(rho))
76    }
77
78    /// Dual power transform with `beta >= 1`.
79    pub fn dual_power(beta: f64) -> Result<Self> {
80        if !beta.is_finite() || beta < 1.0 {
81            return Err(invalid("beta", beta, "must be finite and at least 1"));
82        }
83        Ok(Self::DualPower(beta))
84    }
85
86    /// The spectral measure with exponential risk aversion `k > 0`: the
87    /// larger `k`, the more weight on the worst outcomes.
88    pub fn exponential(k: f64) -> Result<Self> {
89        if !(k.is_finite() && k > 0.0) {
90            return Err(invalid("k", k, "must be finite and positive"));
91        }
92        Ok(Self::Exponential(k))
93    }
94
95    /// The distortion `g(s)` of a survival probability `s` in `[0, 1]`.
96    pub fn g(&self, s: f64) -> f64 {
97        if s <= 0.0 {
98            return 0.0;
99        }
100        if s >= 1.0 {
101            return 1.0;
102        }
103        match *self {
104            Self::Tvar(p) if p >= 1.0 => 1.0,
105            Self::Tvar(p) => (s / (1.0 - p)).min(1.0),
106            Self::Wang(lambda) => norm_cdf(norm_quantile(s) + lambda),
107            Self::ProportionalHazard(rho) => s.powf(rho),
108            Self::DualPower(beta) => -((-s).ln_1p() * beta).exp_m1(),
109            Self::Exponential(k) => (-k * s).exp_m1() / (-k).exp_m1(),
110        }
111    }
112
113    /// Weights for `n` equally likely values sorted ascending: value `i`
114    /// (0-based) gets `g((n - i) / n) - g((n - i - 1) / n)`. They are
115    /// non-negative and sum to 1.
116    ///
117    /// ```
118    /// use prospicio_prob::Distortion;
119    ///
120    /// // TVaR at 50% on four values: the top two, equally.
121    /// assert_eq!(Distortion::tvar(0.5).unwrap().weights(4), [0.0, 0.0, 0.5, 0.5]);
122    /// ```
123    pub fn weights(&self, n: usize) -> Vec<f64> {
124        let nf = n as f64;
125        (0..n)
126            .map(|i| self.g((n - i) as f64 / nf) - self.g((n - i - 1) as f64 / nf))
127            .collect()
128    }
129
130    /// The risk measure of equally likely draws sorted ascending.
131    ///
132    /// For `Tvar(p)` this agrees with [`tvar_sorted`](crate::risk::tvar_sorted)
133    /// up to rounding, including the fractional weight at the VaR. `sorted`
134    /// must be non-empty and sorted ascending; this is checked only in
135    /// debug builds.
136    pub fn apply_sorted(&self, sorted: &[f64]) -> f64 {
137        debug_assert!(!sorted.is_empty(), "draws must not be empty");
138        debug_assert!(
139            sorted.windows(2).all(|w| w[0] <= w[1]),
140            "draws must be sorted ascending"
141        );
142        self.weights(sorted.len())
143            .iter()
144            .zip(sorted)
145            .map(|(w, x)| w * x)
146            .sum()
147    }
148
149    /// The risk measure of a discrete distribution with ascending `values`
150    /// and their `probs` (which should sum to 1).
151    ///
152    /// Survival probabilities are summed from the top, so small tail
153    /// probabilities keep their precision.
154    ///
155    /// ```
156    /// use prospicio_prob::Distortion;
157    ///
158    /// let ph = Distortion::proportional_hazard(0.5).unwrap();
159    /// // P(X = 0) = 0.75, P(X = 1) = 0.25: ρ = g(0.25) = 0.5.
160    /// assert_eq!(ph.apply_discrete(&[0.0, 1.0], &[0.75, 0.25]), 0.5);
161    /// ```
162    pub fn apply_discrete(&self, values: &[f64], probs: &[f64]) -> f64 {
163        debug_assert_eq!(values.len(), probs.len());
164        debug_assert!(
165            values.windows(2).all(|w| w[0] <= w[1]),
166            "values must be sorted ascending"
167        );
168        let mut above = 0.0; // P(X > x_k)
169        let mut g_above = 0.0;
170        let mut total = 0.0;
171        for (x, p) in values.iter().zip(probs).rev() {
172            let at_or_above = above + p;
173            let g_at_or_above = self.g(at_or_above);
174            total += (g_at_or_above - g_above) * x;
175            above = at_or_above;
176            g_above = g_at_or_above;
177        }
178        total
179    }
180}
181
182fn invalid(name: &'static str, value: f64, reason: &'static str) -> Error {
183    Error::InvalidParameter {
184        name,
185        value,
186        reason,
187    }
188}
189
190#[cfg(test)]
191mod tests {
192    use super::*;
193
194    /// On a uniform, `∫₀¹ φ(u) u du = (e^k (k - 1) + 1) / (k (e^k - 1))`.
195    #[test]
196    fn exponential_spectral_on_a_uniform() {
197        let n = 100_000;
198        let u: Vec<f64> = (0..n).map(|i| (i as f64 + 0.5) / n as f64).collect();
199        for k in [0.5f64, 3.0, 20.0] {
200            let want = (k.exp() * (k - 1.0) + 1.0) / (k * k.exp_m1());
201            let got = Distortion::exponential(k).unwrap().apply_sorted(&u);
202            assert!((got - want).abs() < 1e-6, "{k}: {got} vs {want}");
203        }
204        assert!(Distortion::exponential(0.0).is_err());
205    }
206    use crate::risk::tvar_sorted;
207
208    const X: [f64; 5] = [10.0, 20.0, 30.0, 40.0, 50.0];
209
210    fn all() -> Vec<Distortion> {
211        vec![
212            Distortion::tvar(0.7).unwrap(),
213            Distortion::wang(0.5).unwrap(),
214            Distortion::proportional_hazard(0.6).unwrap(),
215            Distortion::dual_power(2.5).unwrap(),
216        ]
217    }
218
219    #[test]
220    fn tvar_matches_the_risk_module() {
221        for i in 0..=100 {
222            let p = i as f64 / 100.0;
223            let d = Distortion::tvar(p).unwrap().apply_sorted(&X);
224            let t = tvar_sorted(&X, p).unwrap();
225            assert!((d - t).abs() <= 1e-12 * t, "p = {p}: {d} vs {t}");
226        }
227    }
228
229    #[test]
230    fn identity_parameters_give_the_mean() {
231        for d in [
232            Distortion::tvar(0.0).unwrap(),
233            Distortion::wang(0.0).unwrap(),
234            Distortion::proportional_hazard(1.0).unwrap(),
235            Distortion::dual_power(1.0).unwrap(),
236        ] {
237            assert!((d.apply_sorted(&X) - 30.0).abs() < 1e-12, "{d:?}");
238        }
239    }
240
241    #[test]
242    fn weights_are_a_probability_vector_rising_to_the_tail() {
243        for d in all() {
244            let w = d.weights(1000);
245            assert!((w.iter().sum::<f64>() - 1.0).abs() < 1e-12, "{d:?}");
246            assert!(w.iter().all(|&w| w >= 0.0));
247            // Concave g: weights never fall towards the larger values.
248            assert!(w.windows(2).all(|p| p[1] >= p[0] - 1e-15), "{d:?}");
249        }
250    }
251
252    #[test]
253    fn coherence_properties() {
254        let mean = 30.0;
255        for d in all() {
256            let r = d.apply_sorted(&X);
257            assert!(r >= mean && r <= 50.0, "{d:?}: {r}");
258            // Translation and positive scaling.
259            let shifted: Vec<f64> = X.iter().map(|x| 2.0 * x + 7.0).collect();
260            assert!((d.apply_sorted(&shifted) - (2.0 * r + 7.0)).abs() < 1e-12);
261        }
262        // A larger parameter is more conservative.
263        let w = |l| Distortion::wang(l).unwrap().apply_sorted(&X);
264        assert!(w(0.2) < w(0.5) && w(0.5) < w(1.0));
265        let ph = |r| Distortion::proportional_hazard(r).unwrap().apply_sorted(&X);
266        assert!(ph(0.9) < ph(0.5));
267        let dp = |b| Distortion::dual_power(b).unwrap().apply_sorted(&X);
268        assert!(dp(1.5) < dp(3.0));
269    }
270
271    #[test]
272    fn discrete_agrees_with_equal_weights_and_merges_ties() {
273        let p = [0.2; 5];
274        for d in all() {
275            let a = d.apply_discrete(&X, &p);
276            let b = d.apply_sorted(&X);
277            assert!((a - b).abs() < 1e-12, "{d:?}");
278        }
279        // Tied draws and one atom with their total mass are the same.
280        let d = Distortion::dual_power(2.0).unwrap();
281        let ties = d.apply_sorted(&[1.0, 5.0, 5.0, 5.0]);
282        let atom = d.apply_discrete(&[1.0, 5.0], &[0.25, 0.75]);
283        assert!((ties - atom).abs() < 1e-15);
284        // Closed form: g(0.75) × 5 + (1 - g(0.75)) × 1.
285        assert!((atom - (0.9375 * 5.0 + 0.0625)).abs() < 1e-15);
286    }
287
288    #[test]
289    fn g_closed_forms() {
290        let s = 0.3;
291        assert_eq!(Distortion::tvar(0.8).unwrap().g(s), 1.0);
292        assert!((Distortion::tvar(0.4).unwrap().g(s) - 0.5).abs() < 1e-15);
293        assert!((Distortion::proportional_hazard(0.5).unwrap().g(s) - s.sqrt()).abs() < 1e-15);
294        assert!((Distortion::dual_power(2.0).unwrap().g(s) - 0.51).abs() < 1e-15);
295        // Φ(Φ⁻¹(0.5) + 1) = Φ(1).
296        assert!((Distortion::wang(1.0).unwrap().g(0.5) - norm_cdf(1.0)).abs() < 1e-15);
297        assert_eq!(Distortion::tvar(1.0).unwrap().g(1e-300), 1.0);
298        for d in all() {
299            assert_eq!(d.g(0.0), 0.0);
300            assert_eq!(d.g(1.0), 1.0);
301        }
302    }
303
304    #[test]
305    fn rejects_bad_parameters() {
306        assert!(Distortion::tvar(1.5).is_err());
307        assert!(Distortion::wang(-0.1).is_err());
308        assert!(Distortion::wang(f64::INFINITY).is_err());
309        assert!(Distortion::proportional_hazard(0.0).is_err());
310        assert!(Distortion::proportional_hazard(1.5).is_err());
311        assert!(Distortion::dual_power(0.5).is_err());
312    }
313}
prospicio_prob/distortion.rs

prospicio_prob/
distortion.rs
