Struct Gpd 

pub struct Gpd { /* private fields */ }

The generalized Pareto distribution with shape xi, scale beta and location u (0 unless set), on x >= u (and x <= u - beta / xi when xi < 0), as SciPy’s genpareto(c=xi, loc=u, scale=beta).

P(X > x) = (1 + xi (x - u) / beta)^(-1/xi), and exp(-(x - u) / beta) at xi = 0.

As a Severity it has closed-form layer means and second moments for every xi, so it also serves as Riegel’s generalized Pareto for treaty pricing (Gpd::riegel).

use prospicio_prob::{Distribution, evt::Gpd};

let g = Gpd::new(0.5, 2.0).unwrap();
assert!((g.mean() - 4.0).abs() < 1e-12); // beta / (1 - xi)
assert!((g.cdf(g.quantile(0.99).unwrap()) - 0.99).abs() < 1e-12);

Implementations

impl Gpd

pub fn new(xi: f64, beta: f64) -> Result<Self>

pub fn shifted(self, location: f64) -> Result<Self>

The same distribution shifted to start at location (finite).

pub fn riegel(t: f64, alpha_ini: f64, alpha_tail: f64) -> Result<Self>

Riegel’s generalized Pareto with threshold t, initial alpha alpha_ini and tail alpha alpha_tail (all positive):

P(X > x) = (1 + (alpha_ini / alpha_tail) (x / t - 1))^(-alpha_tail),  x >= t,

whose local Pareto alpha moves from alpha_ini at t to alpha_tail as x grows. It is this GPD with xi = 1/alpha_tail, beta = t/alpha_ini and location t, and matches pGenPareto(x, t, alpha_ini, alpha_tail) in the R package Pareto.

use prospicio_prob::{Distribution, Severity, evt::Gpd};

let g = Gpd::riegel(1000.0, 2.0, 1.5).unwrap();
// P(X > 2000) = (1 + 2/1.5)^-1.5.
assert!((g.survival(2000.0) - (7.0f64 / 3.0).powf(-1.5)).abs() < 1e-15);
assert!(g.layer(4000.0, 1000.0) > 0.0);

pub fn xi(&self) -> f64

Shape ξ: heavier tails as it grows; moments of order 1/ξ and above are infinite.

pub fn beta(&self) -> f64

Scale β.

pub fn location(&self) -> f64

Location u, where the support starts.

pub fn fit(exceedances: &[f64]) -> Result<Self>

Maximum likelihood fit to exceedances (values over a threshold, minus the threshold), all non-negative.

Maximizes the profile likelihood in θ = ξ / β (Grimshaw 1993): for fixed θ the likelihood is maximized by ξ(θ) = mean(ln(1 + θ y)), leaving a one-dimensional search: a scan over θ finds the maximum, then bisection on the score pins it. ξ is restricted to ξ > -1, where the maximum likelihood estimate exists. Needs at least 3 exceedances, not all equal.

use prospicio_core::StreamRng;
use prospicio_prob::{Distribution, evt::Gpd};

let truth = Gpd::new(0.3, 10.0).unwrap();
let y = truth.sample(&mut StreamRng::new(1, 0), 20_000);
let fit = Gpd::fit(&y).unwrap();
assert!((fit.xi() - 0.3).abs() < 0.05);
assert!((fit.beta() / 10.0 - 1.0).abs() < 0.05);

impl Gpd

pub fn fit_riegel(t: f64, data: &LargeLosses) -> Result<Self>

Maximum likelihood fit of Riegel’s generalized Pareto (Gpd::riegel) with threshold t to losses at or above t, with the reporting thresholds, censoring and weights of data.

With z = x/t − 1 and k = α_ini / α_tail, S(x) = (1 + k z)^(−α_tail). For fixed k the likelihood is maximized by α_tail(k) = Σ_uncensored w / Σ w ln((1 + k z_y) / (1 + k z_r)), so the fit is a one-dimensional profile likelihood in k: a scan on a log grid over [1e-6, 1e6] brackets the maximum, and bisection on the analytic score pins it. k = 1 is the Pareto.

use prospicio_core::StreamRng;
use prospicio_prob::{Distribution, LargeLosses, evt::Gpd};

let truth = Gpd::riegel(1000.0, 3.0, 1.5).unwrap();
let draws = truth.sample(&mut StreamRng::new(4, 0), 100_000);
let fit = Gpd::fit_riegel(1000.0, &LargeLosses::new(draws).unwrap()).unwrap();
// ξ = 1/α_tail, β = t/α_ini.
assert!((1.0 / fit.xi() / 1.5 - 1.0).abs() < 0.05);
assert!((1000.0 / fit.beta() / 3.0 - 1.0).abs() < 0.05);

Trait Implementations

impl Clone for Gpd

fn clone(&self) -> Gpd

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Copy for Gpd

impl Debug for Gpd

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Distribution for Gpd

fn mean(&self) -> f64

u + beta / (1 - xi); infinite for xi >= 1.

fn variance(&self) -> f64

beta^2 / ((1 - xi)^2 (1 - 2 xi)); infinite for xi >= 1/2.

fn cdf(&self, x: f64) -> f64

P(X <= x).

fn survival(&self, x: f64) -> f64

P(X > x). Representations with a direct form override the default 1 - cdf(x), which loses all precision far in the tail.

fn quantile(&self, p: f64) -> Result<f64>

Smallest x with cdf(x) >= p.

fn std_dev(&self) -> f64

Standard deviation.

fn sample(&self, rng: &mut StreamRng, n: usize) -> Vec<f64>

n draws from stream rng, by inverse transform unless the family overrides it (crate::Gamma draws by Marsaglia and Tsang).

fn is_parallel_safe(&self) -> bool

Whether the distribution may be evaluated on several threads at once. True for every native family; false for a crate::Custom whose callbacks must stay on the calling thread (an R function), so the parallel simulations run single-threaded when they meet one.

impl From<Gpd> for Dist

fn from(d: Gpd) -> Self

Converts to this type from the input type.

impl PartialEq for Gpd

fn eq(&self, other: &Gpd) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl Severity for Gpd

fn layer(&self, limit: f64, attachment: f64) -> f64

∫_a^b S(x) dx: the part below the location pays in full.

fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64

2 ∫_a^b (x − a) S(x) dx.

fn lev(&self, limit: f64) -> f64

Limited expected value E[min(X, limit)].

fn stop_loss(&self, retention: f64) -> f64

Expected excess over a retention, E[max(X - retention, 0)].

fn layer_variance(&self, limit: f64, attachment: f64) -> f64

Variance of the loss to the layer limit xs attachment.

impl StructuralPartialEq for Gpd