prospicio_prob/
pareto.rs

1//! The single-parameter (European) Pareto distribution, optionally
2//! truncated: the standard large-loss severity in reinsurance pricing (see
3//! `docs/design/pareto.md`).
4
5use prospicio_core::{Error, Result};
6
7use crate::distribution::{Distribution, check_probability};
8use crate::severity::Severity;
9
10/// Pareto distribution `Pareto(t, α)`: `P(X > x) = (t/x)^α` for `x ≥ t`,
11/// and 1 below `t`. Also known as Pareto type I or single-parameter
12/// Pareto; matches the R package Pareto (`pPareto(x, t, alpha)`).
13///
14/// With a `truncation` `T > t` it is the Pareto conditioned on `X < T`:
15/// `P(X > x) = ((t/x)^α − (t/T)^α) / (1 − (t/T)^α)` for `t ≤ x < T`.
16///
17/// Layer integrals are in closed form for every `α > 0`, including the
18/// logarithmic cases `α = 1` (limited mean) and `α = 2` (limited second
19/// moment), which need no special handling: integrals of `x^{k−α}` are
20/// evaluated as `a^{k+1−α} · L · exprel((k + 1 − α) L)` with
21/// `L = ln(b/a)` and `exprel(z) = (e^z − 1)/z`.
22///
23/// ```
24/// use prospicio_prob::{Distribution, Pareto, Severity};
25///
26/// let p = Pareto::new(500.0, 2.0).unwrap();
27/// // The R package's vignette example: layer 4000 xs 1000.
28/// assert!((p.layer(4000.0, 1000.0) - 200.0).abs() < 1e-12);
29/// assert!((p.mean() - 1000.0).abs() < 1e-12);
30/// ```
31#[derive(Debug, Clone, Copy, PartialEq)]
32pub struct Pareto {
33    t: f64,
34    alpha: f64,
35    truncation: Option<f64>,
36}
37
38impl Pareto {
39    /// Pareto with threshold `t > 0` and alpha `α > 0`.
40    pub fn new(t: f64, alpha: f64) -> Result<Self> {
41        if !t.is_finite() || t <= 0.0 {
42            return Err(invalid("t", t, "must be finite and positive"));
43        }
44        if !alpha.is_finite() || alpha <= 0.0 {
45            return Err(invalid("alpha", alpha, "must be finite and positive"));
46        }
47        Ok(Self {
48            t,
49            alpha,
50            truncation: None,
51        })
52    }
53
54    /// The same Pareto conditioned on `X < truncation`, with
55    /// `t < truncation < ∞`.
56    pub fn truncated(self, truncation: f64) -> Result<Self> {
57        if !truncation.is_finite() || truncation <= self.t {
58            return Err(invalid(
59                "truncation",
60                truncation,
61                "must be finite and above t",
62            ));
63        }
64        Ok(Self {
65            truncation: Some(truncation),
66            ..self
67        })
68    }
69
70    /// Threshold `t`.
71    pub fn t(&self) -> f64 {
72        self.t
73    }
74
75    /// Pareto alpha.
76    pub fn alpha(&self) -> f64 {
77        self.alpha
78    }
79
80    /// Truncation point, if any.
81    pub fn truncation(&self) -> Option<f64> {
82        self.truncation
83    }
84
85    /// `P(X > x)`.
86    fn survival_at(&self, x: f64) -> f64 {
87        if x < self.t {
88            return 1.0;
89        }
90        let raw = (self.t / x).powf(self.alpha);
91        match self.truncation {
92            None => raw,
93            Some(tr) if x >= tr => 0.0,
94            Some(tr) => {
95                // (t/x)^α − (t/T)^α = (t/x)^α (1 − (x/T)^α), without cancellation.
96                let (_, one_minus_q) = self.truncated_mass(tr);
97                raw * -(self.alpha * (x / tr).ln()).exp_m1() / one_minus_q
98            }
99        }
100    }
101
102    /// `(q, 1 − q)` with `q = (t/T)^α`, `1 − q` without cancellation.
103    fn truncated_mass(&self, tr: f64) -> (f64, f64) {
104        let log_q = self.alpha * (self.t / tr).ln();
105        (log_q.exp(), -log_q.exp_m1())
106    }
107
108    /// `∫_a^b x^k S(x) dx` for `k ∈ {0, 1}` and `0 ≤ a ≤ b ≤ ∞`, with `S`
109    /// the survival function (truncated if set).
110    fn integral(&self, k: i32, a: f64, b: f64) -> f64 {
111        match self.truncation {
112            None => self.raw_integral(k, a, b),
113            Some(tr) => {
114                let (a, b) = (a.min(tr), b.min(tr));
115                if a >= b {
116                    return 0.0;
117                }
118                let (q, one_minus_q) = self.truncated_mass(tr);
119                (self.raw_integral(k, a, b) - q * power_integral(k, a, b)) / one_minus_q
120            }
121        }
122    }
123
124    /// `∫_a^b x^k S(x) dx` for the untruncated Pareto.
125    fn raw_integral(&self, k: i32, a: f64, b: f64) -> f64 {
126        raw_integral(k, self.t, self.alpha, a, b)
127    }
128}
129
130/// `∫_a^b x^k S(x) dx` for `S` the untruncated `Pareto(t, α)` survival
131/// function, `k ∈ {0, 1}`, `0 ≤ a ≤ b ≤ ∞` and `α ≥ 0`.
132pub(crate) fn raw_integral(k: i32, t: f64, alpha: f64, a: f64, b: f64) -> f64 {
133    if a >= b {
134        return 0.0;
135    }
136    let below = if a < t {
137        power_integral(k, a, b.min(t))
138    } else {
139        0.0
140    };
141    let lo = a.max(t);
142    if lo >= b {
143        return below;
144    }
145    // t^α ∫_lo^b x^(k−α) dx = (t/lo)^α lo^(k+1) ∫_1^(b/lo) u^(k−α) du.
146    let e = f64::from(k) + 1.0 - alpha;
147    let scale = (t / lo).powf(alpha) * lo.powi(k + 1);
148    let above = if b == f64::INFINITY {
149        if e < 0.0 { scale / -e } else { f64::INFINITY }
150    } else {
151        let l = (b / lo).ln();
152        scale * l * exprel(e * l)
153    };
154    below + above
155}
156
157impl Distribution for Pareto {
158    /// `α t / (α − 1)` untruncated; infinite for `α ≤ 1`.
159    fn mean(&self) -> f64 {
160        match self.truncation {
161            None if self.alpha > 1.0 => self.alpha * self.t / (self.alpha - 1.0),
162            None => f64::INFINITY,
163            Some(_) => self.integral(0, 0.0, f64::INFINITY),
164        }
165    }
166
167    /// `α t² / ((α − 1)² (α − 2))` untruncated; infinite for `α ≤ 2`.
168    fn variance(&self) -> f64 {
169        match self.truncation {
170            None if self.alpha > 2.0 => {
171                let a = self.alpha;
172                a * self.t * self.t / ((a - 1.0) * (a - 1.0) * (a - 2.0))
173            }
174            None => f64::INFINITY,
175            Some(_) => {
176                let m = self.mean();
177                2.0 * self.integral(1, 0.0, f64::INFINITY) - m * m
178            }
179        }
180    }
181
182    fn cdf(&self, x: f64) -> f64 {
183        1.0 - self.survival_at(x)
184    }
185
186    fn survival(&self, x: f64) -> f64 {
187        self.survival_at(x)
188    }
189
190    /// `t (1 − p)^(−1/α)` untruncated; `p = 1` gives `+∞` (or the
191    /// truncation point).
192    fn quantile(&self, p: f64) -> Result<f64> {
193        check_probability(p)?;
194        let s = match self.truncation {
195            None => 1.0 - p,
196            Some(tr) => {
197                let (q, one_minus_q) = self.truncated_mass(tr);
198                q + (1.0 - p) * one_minus_q
199            }
200        };
201        if s <= 0.0 {
202            return Ok(f64::INFINITY);
203        }
204        Ok(self.t * s.powf(-1.0 / self.alpha))
205    }
206}
207
208impl Severity for Pareto {
209    fn lev(&self, limit: f64) -> f64 {
210        if limit <= 0.0 {
211            return limit;
212        }
213        self.integral(0, 0.0, limit)
214    }
215
216    fn stop_loss(&self, retention: f64) -> f64 {
217        if retention <= 0.0 {
218            return self.mean() - retention;
219        }
220        self.integral(0, retention, f64::INFINITY)
221    }
222
223    fn layer(&self, limit: f64, attachment: f64) -> f64 {
224        let a = attachment.max(0.0);
225        self.integral(0, a, a + limit)
226    }
227
228    /// `2 ∫_a^b (x − a) S(x) dx` with `b = a + limit`.
229    fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64 {
230        let a = attachment.max(0.0);
231        let b = a + limit;
232        2.0 * (self.integral(1, a, b) - a * self.integral(0, a, b))
233    }
234}
235
236/// `∫_a^b x^k dx` for `k ∈ {0, 1}`.
237pub(crate) fn power_integral(k: i32, a: f64, b: f64) -> f64 {
238    match k {
239        0 => b - a,
240        _ => 0.5 * (b - a) * (b + a),
241    }
242}
243
244/// `(e^z − 1) / z`, and 1 at `z = 0`.
245fn exprel(z: f64) -> f64 {
246    if z == 0.0 { 1.0 } else { z.exp_m1() / z }
247}
248
249pub(crate) fn invalid(name: &'static str, value: f64, reason: &'static str) -> Error {
250    Error::InvalidParameter {
251        name,
252        value,
253        reason,
254    }
255}
256
257#[cfg(test)]
258mod tests {
259    use super::*;
260
261    fn close(a: f64, b: f64, rel: f64) -> bool {
262        (a - b).abs() <= rel * b.abs().max(1e-300)
263    }
264
265    #[test]
266    fn layer_means_match_the_closed_forms() {
267        // I(a, b) = t^α (b^(1−α) − a^(1−α)) / (1 − α), and t ln(b/a) at α = 1.
268        for alpha in [0.5, 1.0, 1.5, 2.0, 3.0] {
269            let p = Pareto::new(500.0, alpha).unwrap();
270            let (a, b) = (1000.0f64, 5000.0f64);
271            let want = if alpha == 1.0 {
272                500.0 * (b / a).ln()
273            } else {
274                500f64.powf(alpha) * (b.powf(1.0 - alpha) - a.powf(1.0 - alpha)) / (1.0 - alpha)
275            };
276            assert!(close(p.layer(4000.0, 1000.0), want, 1e-13), "alpha {alpha}");
277        }
278        // Below the threshold every loss pays in full.
279        let p = Pareto::new(500.0, 2.0).unwrap();
280        assert_eq!(p.layer(300.0, 100.0), 300.0);
281        assert!(close(
282            p.lev(1000.0),
283            500.0 + 500.0 * 500.0 * (1.0 / 500.0 - 1.0 / 1000.0),
284            1e-14
285        ));
286    }
287
288    #[test]
289    fn severity_identities() {
290        for p in [
291            Pareto::new(1000.0, 2.5).unwrap(),
292            Pareto::new(1000.0, 1.5)
293                .unwrap()
294                .truncated(20_000.0)
295                .unwrap(),
296            Pareto::new(1000.0, 0.8).unwrap().truncated(1e6).unwrap(),
297        ] {
298            for d in [500.0, 1000.0, 3000.0, 15_000.0] {
299                assert!(
300                    close(p.lev(d) + p.stop_loss(d), p.mean(), 1e-12),
301                    "{p:?} at {d}"
302                );
303                let layer = p.layer(4000.0, d);
304                assert!(close(
305                    layer,
306                    p.stop_loss(d) - p.stop_loss(d + 4000.0),
307                    1e-10
308                ));
309            }
310            // Second moment of the whole loss gives the variance.
311            let m2 = p.layer_second_moment(f64::INFINITY, 0.0);
312            assert!(
313                close(m2 - p.mean() * p.mean(), p.variance(), 1e-10),
314                "{p:?}"
315            );
316            assert!(close(
317                p.layer_variance(f64::INFINITY, 0.0),
318                p.variance(),
319                1e-10
320            ));
321        }
322    }
323
324    #[test]
325    fn untruncated_moments_and_infinite_cases() {
326        let p = Pareto::new(100.0, 3.0).unwrap();
327        assert!(close(p.mean(), 150.0, 1e-15));
328        // E[X^2] = α t^2 / (α − 2) = 30,000; variance 30,000 − 150^2.
329        assert!(close(p.variance(), 7500.0, 1e-15));
330        assert!(close(
331            p.layer_second_moment(f64::INFINITY, 0.0),
332            3.0e4,
333            1e-13
334        ));
335        let heavy = Pareto::new(100.0, 1.0).unwrap();
336        assert_eq!(heavy.mean(), f64::INFINITY);
337        assert_eq!(heavy.stop_loss(1000.0), f64::INFINITY);
338        assert!(heavy.layer(1000.0, 1000.0).is_finite());
339        assert_eq!(Pareto::new(100.0, 2.0).unwrap().variance(), f64::INFINITY);
340        // α = 2: the limited second moment has a logarithm.
341        let p2 = Pareto::new(1.0, 2.0).unwrap();
342        // 2 ∫_1^e x · x^-2 dx + 2 ∫_0^1 x dx = 2 + 1.
343        assert!(close(
344            p2.layer_second_moment(std::f64::consts::E, 0.0),
345            3.0,
346            1e-14
347        ));
348    }
349
350    #[test]
351    fn cdf_quantile_round_trip() {
352        for p in [
353            Pareto::new(1000.0, 2.0).unwrap(),
354            Pareto::new(1000.0, 0.7)
355                .unwrap()
356                .truncated(50_000.0)
357                .unwrap(),
358        ] {
359            for q in [0.0, 0.1, 0.5, 0.9, 0.999] {
360                let x = p.quantile(q).unwrap();
361                assert!(
362                    close(p.cdf(x), q, 1e-12) || (q == 0.0 && p.cdf(x) == 0.0),
363                    "{p:?} {q}"
364                );
365            }
366            assert_eq!(p.cdf(999.0), 0.0);
367        }
368        assert_eq!(
369            Pareto::new(1.0, 2.0).unwrap().quantile(1.0).unwrap(),
370            f64::INFINITY
371        );
372        let t = Pareto::new(1.0, 2.0).unwrap().truncated(10.0).unwrap();
373        assert!(close(t.quantile(1.0).unwrap(), 10.0, 1e-15));
374        assert_eq!(t.survival(10.0), 0.0);
375    }
376
377    #[test]
378    fn large_alpha_does_not_overflow() {
379        let p = Pareto::new(1e6, 80.0).unwrap();
380        let l = p.layer(1e6, 2e6);
381        assert!(l.is_finite() && l > 0.0 && l < 1e-15);
382        assert!(close(p.mean(), 80.0 * 1e6 / 79.0, 1e-15));
383    }
384
385    #[test]
386    fn rejects_bad_parameters() {
387        assert!(Pareto::new(0.0, 1.0).is_err());
388        assert!(Pareto::new(1.0, 0.0).is_err());
389        assert!(Pareto::new(1.0, f64::INFINITY).is_err());
390        assert!(Pareto::new(10.0, 1.0).unwrap().truncated(5.0).is_err());
391    }
392}
prospicio_prob/pareto.rs

prospicio_prob/
pareto.rs
