---------------------------------------------------------------------- This is the API documentation for the prospicio library. ---------------------------------------------------------------------- ## Distributions Parametric, discretized and sampled representations. Lognormal(meanlog, sdlog) Lognormal distribution: ``ln X ~ Normal(meanlog, sdlog**2)``. Parameters ---------- meanlog : float Mean of ``ln X``. sdlog : float Standard deviation of ``ln X``; must be positive. Raises ------ ValueError If ``meanlog`` is not finite or ``sdlog`` is not positive and finite. Examples -------- >>> from prospicio.distributions import Lognormal >>> d = Lognormal.from_mean_cv(1000.0, 0.5) >>> round(d.mean(), 6) 1000.0 Gamma(shape, scale) Gamma distribution with shape ``alpha`` and scale ``theta``: mean ``alpha * theta``, variance ``alpha * theta**2``. Parameters ---------- shape : float scale : float Raises ------ ValueError If a parameter is not finite and positive. Examples -------- >>> from prospicio.distributions import Gamma >>> g = Gamma.from_mean_cv(1000.0, 0.5) >>> g.shape, round(g.std(), 9) (4.0, 500.0) Tweedie(mean, dispersion, power) Tweedie distribution with mean ``mu``, dispersion ``phi`` and power ``1 < p < 2``: variance ``phi * mu**p``, a point mass at 0 and a continuous density above it. It is a Poisson number of gamma losses, the GLM family for pure premium. ``P(Y = 0) = exp(-lambda_)``. Parameters ---------- mean : float dispersion : float power : float In ``(1, 2)``. Raises ------ ValueError If a parameter is out of range. Examples -------- >>> import math >>> from prospicio.distributions import Tweedie >>> y = Tweedie(500.0, 40.0, 1.6) >>> abs(y.cdf(0.0) - math.exp(-y.lambda_)) < 1e-15 True Weibull(shape, scale) Weibull distribution with shape ``k`` and scale ``lam``: ``P(X > x) = exp(-(x / lam) ** k)``, as SciPy's ``weibull_min``. Parameters ---------- shape : float scale : float Examples -------- >>> from prospicio.distributions import Weibull >>> Weibull(1.0, 2.0).mean() 2.0 Loglogistic(shape, scale) Loglogistic (Fisk) distribution with shape ``alpha`` and scale ``theta`` (the median): ``F(x) = (x/theta)**alpha / (1 + (x/theta)**alpha)``, as SciPy's ``fisk``. Its ``cdf`` is Clark's loglogistic growth curve. The mean is infinite for ``alpha <= 1`` and the variance for ``alpha <= 2``; limited and layer moments always exist. Parameters ---------- shape : float scale : float Examples -------- >>> from prospicio.distributions import Loglogistic >>> Loglogistic(1.0, 2.0).cdf(3.0) 0.6 Mixture(components) A finite mixture of severities: component ``i`` with probability ``w_i``, such as attritional plus large losses. Parameters ---------- components : list of (float, severity) Weights (positive, summing to 1) and severities. Examples -------- >>> from prospicio.distributions import Lognormal, Mixture, Pareto >>> m = Mixture([(0.9, Lognormal.from_mean_cv(1e4, 1.0)), (0.1, Pareto(1e5, 2.0))]) >>> round(m.mean(), 6) 29000.0 Custom(cdf, quantile=None, name='custom') A loss severity defined by your own distribution function: the slow path for a distribution the library does not have. Give the cdf, and the quantile function if you have one (sampling inverts the cdf by bisection otherwise, about a hundred cdf calls per draw). The mean, variance, limited expected values and layer moments are computed by Gauss–Legendre quadrature of the survival function between the distribution's own quantiles, ignoring the probability above the ``1 - 1e-12`` quantile. A ``Custom`` goes anywhere a severity does (layers, compound distributions, simulated events, copula marginals, mixtures); calculations that meet one run single-threaded, since every value calls back into Python. Parameters ---------- cdf : callable ``cdf(x) -> float``: ``P(X <= x)`` for ``x >= 0``, in ``[0, 1]`` and non-decreasing. Losses are non-negative. quantile : callable, optional ``quantile(p) -> float``: the smallest ``x`` with ``cdf(x) >= p``. name : str, default "custom" Shown in errors and ``repr``. Raises ------ ValueError If a callable raises or returns a value out of range, or the cdf never reaches ``1 - 1e-12`` at a finite loss. An error in a later call makes that value ``nan``; ``last_error`` says why. Examples -------- >>> import math >>> from prospicio.distributions import Custom >>> d = Custom(lambda x: 1 - math.exp(-x / 100), name="exponential") >>> round(d.mean(), 6) 100.0 >>> round(d.lev(50), 6) == round(100 * (1 - math.exp(-0.5)), 6) True to_json(dist) A distribution as a JSON document: the family and the parameters its constructor takes, versioned, numbers bit for bit. ``from_json`` reads it back to an equal distribution of the same class. A ``Custom`` cannot be saved: it is a Python function. Parameters ---------- dist : a distribution Any distribution class, ``Sampled`` and ``Mixture`` included. Returns ------- str Raises ------ ValueError For a ``Custom``. Examples -------- >>> from prospicio.distributions import Lognormal, from_json, to_json >>> text = to_json(Lognormal(7.0, 0.5)) >>> from_json(text).mean() == Lognormal(7.0, 0.5).mean() True from_json(text) A distribution from a document written by ``to_json``, as the class of its family. Parameters ---------- text : str Returns ------- a distribution Raises ------ ValueError If the document is malformed, of another format or a newer version, or its parameters are out of range. Grid(step, probs) A distribution on the points ``0, step, 2*step, ...``: the discretized representation that FFT and Panjer aggregation work on. Parameters ---------- step : float Grid step; must be positive. probs : list of float Probabilities at ``0, step, ...``; non-negative, summing to 1. Raises ------ ValueError If the step is not positive or the probabilities are invalid. Examples -------- >>> from prospicio.distributions import Grid, Lognormal >>> grid, report = Grid.local_moment(Lognormal(7.0, 0.5), 100.0, 200) >>> report.tail_mass < 1e-8 True DiscretizationReport How a distribution was discretized, and the error that introduced. Returned with the grid by ``Grid.local_moment``, ``Grid.rounding`` and ``Grid.lower``. Sampled(draws) A distribution known only through equally weighted draws. Parameters ---------- draws : list of float Non-empty, all finite. Raises ------ ValueError If ``draws`` is empty or holds a value that is not finite. Examples -------- >>> from prospicio.distributions import Sampled >>> s = Sampled([1.0, 2.0, 3.0, 4.0]) >>> s.tvar(0.5) 3.5 ## Large-loss severities The Pareto family for layer and treaty pricing, with maximum likelihood fits. Pareto(t, alpha, truncation=None) Single-parameter Pareto: ``P(X > x) = (t / x) ** alpha`` for ``x >= t``, optionally truncated (conditioned on ``X < truncation``). Parameters ---------- t : float Threshold; finite and positive. alpha : float Pareto alpha; finite and positive. truncation : float, optional Truncation point above ``t``. Raises ------ ValueError If a parameter is out of range. Examples -------- >>> from prospicio.distributions import Pareto >>> p = Pareto(500.0, 2.0) >>> round(p.layer(4000.0, 1000.0), 9) 200.0 PiecewisePareto(t, alpha, truncation=None, truncation_type='lp') Piecewise Pareto: alpha ``alpha[k]`` above threshold ``t[k]``, the general large-loss model and the result of tower matching. Parameters ---------- t : list of float Strictly increasing positive thresholds. alpha : list of float One alpha per threshold; interior ones may be 0, the last must be positive. truncation : float, optional Truncation point above the last threshold. truncation_type : {"lp", "wd"}, default "lp" Truncate the last piece only, or the whole distribution. Raises ------ ValueError If a parameter is out of range. Examples -------- >>> from prospicio.distributions import PiecewisePareto >>> pp = PiecewisePareto([1000.0, 2000.0], [1.0, 2.0]) >>> round(pp.survival(4000.0), 12) 0.125 LogAffinePareto(t, alpha0, gamma) Log-affine local Pareto: the local alpha ``alpha0 * (1 + gamma * ln(x / t))`` rises linearly in the log of the amount, so ``P(X > x) = exp(-alpha0 L - alpha0 gamma L**2 / 2)`` with ``L = ln(x / t)``. Parameters ---------- t : float Threshold; finite and positive. alpha0 : float Local alpha at ``t``; finite and positive. gamma : float Non-negative; 0 gives the Pareto. Examples -------- >>> from prospicio.distributions import LogAffinePareto >>> d = LogAffinePareto.from_delta(1e6, 1.5, 0.5) >>> round(d.local_alpha(2e6), 12) 2.0 GeneralizedPareto(xi, beta, location=0.0) Generalized Pareto severity with a location (Riegel's parameterization via :meth:`GeneralizedPareto.riegel`): ``P(X > x) = (1 + xi (x - location) / beta) ** (-1 / xi)`` above the location. For tail estimation from draws, see :class:`prospicio.risk.Gpd`; this class is the same distribution as a pricing severity. Parameters ---------- xi : float Shape. beta : float Scale; finite and positive. location : float, default 0.0 Examples -------- >>> from prospicio.distributions import GeneralizedPareto >>> g = GeneralizedPareto.riegel(1000.0, 2.0, 1.5) >>> round(g.survival(2000.0), 12) == round((7 / 3) ** -1.5, 12) True local_pareto_to_piecewise(t, alpha, rel_tolerance=0.0001, stop_survival=1e-09, stop_at=Ellipsis) Converts the local Pareto distribution with local alpha ``alpha(x)`` above ``t`` to a piecewise Pareto that matches its survival function exactly at the thresholds and within ``rel_tolerance`` between them. Parameters ---------- t : float Threshold; ``P(X > x) = 1`` below it. alpha : callable ``alpha(x) -> float``, finite and non-negative, positive where the conversion stops. rel_tolerance : float, default 1e-4 stop_survival : float, default 1e-9 Stop once the survival function falls below this. stop_at : float, default inf Stop at this amount. Returns ------- tuple of (PiecewisePareto, float, float) The approximation, the largest relative error found, and where the approximated range ends (the last alpha continues above it). Examples -------- >>> import math >>> from prospicio.distributions import local_pareto_to_piecewise >>> pp, err, end = local_pareto_to_piecewise(1000.0, lambda x: 1.5 + 0.3 * math.log(x / 1000.0)) >>> err <= 1e-4 True ## Claim counts Frequency distributions for frequency-severity aggregation. Poisson(lam) Poisson claim counts with mean ``lam``. Parameters ---------- lam : float Mean number of claims; must be finite and non-negative. Raises ------ ValueError If ``lam`` is negative or not finite. Examples -------- >>> from prospicio.distributions import Poisson >>> n = Poisson(3.0) >>> round(n.pmf(0), 6) 0.049787 NegativeBinomial(r, beta) Negative binomial claim counts: mean ``r * beta``, variance ``r * beta * (1 + beta)`` (Klugman, Panjer & Willmot). SciPy's ``nbinom(n=r, p=1/(1+beta))`` is the same distribution. Parameters ---------- r : float Shape; must be positive. beta : float Scale; must be positive. Raises ------ ValueError If ``r`` or ``beta`` is not positive and finite. Examples -------- >>> from prospicio.distributions import NegativeBinomial >>> n = NegativeBinomial.from_mean_variance(10.0, 30.0) >>> round(n.variance(), 9) 30.0 Binomial(n, p) Binomial claim counts: ``n`` risks, each claiming with probability ``p``. Parameters ---------- n : int Number of trials. p : float Claim probability in ``[0, 1)``. Examples -------- >>> from prospicio.distributions import Binomial >>> Binomial(10, 0.3).mean() 3.0 claim_count(mean, dispersion) The claim count with this mean and dispersion ``Var[N] / E[N]``: binomial below 1, Poisson at 1, negative binomial above 1. A binomial needs a whole number of trials, so below 1 the trials are ``mean / (1 - dispersion)`` rounded up: the mean is kept and the dispersion moves up to the nearest attainable value. Parameters ---------- mean : float dispersion : float Positive. Returns ------- Binomial or Poisson or NegativeBinomial Examples -------- >>> from prospicio.distributions import claim_count >>> claim_count(4.0, 2.5) NegativeBinomial(r=2.6666666666666665, beta=1.5) ## Predictive distributions The joint result every model returns. PredictiveDistribution(dims, components, draws) The joint result every model returns: draws for each simulation (row) and component (column), keyed by dimension values. The ``mean``, ``quantile``, ``var`` and ``tvar`` methods describe the total over all components, computed from row sums. Parameters ---------- dims : list of str Dimension names, e.g. ``["lob", "origin"]``. components : list of tuple One key per component, with one ``int`` or ``str`` per dimension. draws : list of list of float One row per simulation, one value per component. Raises ------ ValueError If the keys or the draws do not fit together. Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> pd = PredictiveDistribution(["line"], [("A",), ("B",)], ... [[0.0, 0.0], [0.0, 0.0], [0.0, 100.0], [100.0, 0.0]]) >>> pd.var(0.75) 100.0 >>> pd.marginal(("A",)).var(0.75) 0.0 ## Aggregate loss Compound distributions on a grid, and simulated years of losses. panjer(frequency, severity, points) Aggregate loss ``S = X_1 + ... + X_N`` by Panjer's recursion. Parameters ---------- frequency : Poisson, NegativeBinomial or Binomial severity : Grid Severity on a grid; the result uses its step. points : int Points in the aggregate grid. Returns ------- tuple of (Grid, CompoundReport) Raises ------ ValueError If ``points`` is 0 or ``P(S = 0)`` underflows (use ``fft``). Examples -------- >>> from prospicio.aggregate import panjer >>> from prospicio.distributions import Grid, Poisson >>> sev = Grid(1.0, [0.1, 0.3, 0.25, 0.2, 0.1, 0.05]) >>> agg, report = panjer(Poisson(3.0), sev, 100) >>> round(agg.mean(), 6) 6.15 fft(frequency, severity, points) Aggregate loss ``S = X_1 + ... + X_N`` by fast Fourier transform. Works for large claim counts that make ``panjer`` underflow. Check ``report.aliasing_error`` before using the result. Parameters ---------- frequency : Poisson, NegativeBinomial or Binomial severity : Grid points : int Returns ------- tuple of (Grid, CompoundReport) Raises ------ ValueError If ``points`` is 0. CompoundReport What a compound calculation produced and the error it introduced. Returned with the aggregate grid by ``panjer`` and ``fft``. Read ``aliasing_error`` first for FFT results: when it is not negligible the grid is unreliable, including ``tail_mass``. simulate_events(frequency, severity, n_sims, seed) Simulates ``n_sims`` years of claims: a count from ``frequency``, then that many independent losses from ``severity``. Parameters ---------- frequency : Poisson, NegativeBinomial or Binomial severity : Lognormal, Grid, Pareto, PiecewisePareto, LogAffinePareto or GeneralizedPareto n_sims : int seed : int Returns ------- EventSet Raises ------ ValueError If ``n_sims`` is 0. Examples -------- >>> from prospicio.aggregate import simulate_events >>> from prospicio.distributions import Lognormal, Poisson >>> events = simulate_events(Poisson(5.0), Lognormal.from_mean_cv(1000.0, 1.0), 20_000, 42) >>> abs(events.totals().mean() - 5000.0) < 75.0 True EventSet Simulated years of individual losses, for applying per-loss terms such as reinsurance layers. Created by ``simulate_events``. Year ``i`` was drawn from stream ``i`` of the generator keyed by ``seed``, so results do not depend on the number of threads. ## Reinsurance Excess-of-loss layers and towers, applied to simulated losses or exactly on the grid. Layer(name, limit, attachment, share=1.0, aggregate_deductible=0.0, aggregate_limit=None, reinstatements=None, premium=0.0, reinstatement_rates=None, pro_rata_time=False) A per-occurrence excess-of-loss layer: ``limit`` xs ``attachment`` on each loss, then annual terms. For one year, ``ceded = share * min(max(sum of per-loss recoveries - aggregate_deductible, 0), aggregate_limit)``, less any loss corridor (``with_loss_corridor``) before the annual limit. Parameters ---------- name : str limit : float Per-occurrence limit; may be ``inf``. attachment : float share : float, default 1.0 Placed share, in ``(0, 1]``. aggregate_deductible : float, default 0.0 aggregate_limit : float, default inf reinstatements : int, optional Free reinstatements: sets ``aggregate_limit`` to ``limit * (reinstatements + 1)``; cannot be combined with ``aggregate_limit``. premium : float, default 0.0 Upfront premium for the placed share; used only by ``reinstatement_rates``. reinstatement_rates : list of float, optional Paid reinstatements, one rate per reinstatement as a fraction of ``premium`` (1.0 is 100%), pro rata as to amount. Sets ``aggregate_limit`` to ``limit * (len(reinstatement_rates) + 1)``; cannot be combined with ``aggregate_limit`` or ``reinstatements``. pro_rata_time : bool, default False Paid reinstatements also pro rata as to time: the limit a loss at time ``t`` (the fraction of the year elapsed) uses up is charged at ``1 - t``. Needs ``reinstatement_rates``, and events with times (``EventSet.with_uniform_times``, or ``times=`` in ``EventSet.from_years``). Raises ------ ValueError If a term is out of range, or more than one of ``aggregate_limit``, ``reinstatements`` and ``reinstatement_rates`` is given. Examples -------- >>> from prospicio.reinsurance import Layer >>> layer = Layer("5x5", 5e6, 5e6, reinstatements=1) >>> layer.ceded([7e6]) 2000000.0 >>> layer.ceded([12e6, 20e6, 30e6]) 10000000.0 >>> paid = Layer("10x10", 10.0, 10.0, premium=2.0, reinstatement_rates=[1.0, 0.5]) >>> paid.reinstatement_premium([22.0, 12.0]) 2.2 >>> timed = Layer("10x10", 10.0, 10.0, premium=2.0, reinstatement_rates=[1.0, 0.5], ... pro_rata_time=True) >>> round(timed.reinstatement_premium([22.0, 12.0], times=[0.25, 0.5]), 12) 1.6 Tower(layers) A reinsurance programme: layers in inuring stages. ``Tower(layers)`` is one stage: every layer sees the gross losses. ``Tower.inuring(stages)`` applies stages in order, each seeing the losses net of all earlier stages, event by event. Parameters ---------- layers : list of Layer At least one; names must be unique. Raises ------ ValueError If there are no layers or two share a name. Examples -------- >>> from prospicio.aggregate import simulate_events >>> from prospicio.reinsurance import Layer, Tower >>> from prospicio.distributions import Lognormal, Poisson >>> events = simulate_events(Poisson(2.0), Lognormal.from_mean_cv(3e6, 1.5), 1_000, 7) >>> tower = Tower([Layer("5x5", 5e6, 5e6), Layer("15x10", 15e6, 10e6)]) >>> result = tower.apply(events) >>> [k[0] for k in result.aggregate(["kind"]).components()] ['gross', 'ceded', 'net'] TowerGrids A tower's annual distributions on the grid, from ``Tower.on_grid``. Every grid is a marginal distribution. Ceded grids are at the placed share and after annual terms; a layer with share ``c`` has step ``c * h``. ## Models Terms and design matrices, GLMs and GAMs, metrics, resampling and MCMC diagnostics. Terms() The terms of a model: an intercept, numeric columns and factors. Build them up, then ``fit`` them to training data to learn the factor levels; the result builds the same design matrix on any data. Examples -------- >>> from prospicio.models import Terms >>> data = {"age": [30.0, 45.0, 60.0], "region": ["N", "S", "W"]} >>> coding = Terms().intercept().numeric("age").factor("region").fit(data) >>> coding.names ['(Intercept)', 'age', 'region[S]', 'region[W]'] Coding Terms with factor levels learned from training data, from ``Terms.fit``. Design(columns, names, offset=None, weights=None) A design matrix with an offset and prior weights. Parameters ---------- columns : list of list of float One list per column. names : list of str offset : list of float, optional weights : list of float, optional Examples -------- >>> from prospicio.models import Design >>> d = Design([[1.0, 1.0], [0.0, 2.0]], ["(Intercept)", "x"]) >>> d.n_rows, d.names (2, ['(Intercept)', 'x']) Glm(family, link=None, dispersion=None, theta=None, power=None, link_power=None) A generalized linear model, fitted by IRLS. Parameters ---------- family : str ``"gaussian"``, ``"poisson"``, ``"gamma"``, ``"inverse_gaussian"``, ``"binomial"``, ``"negative_binomial"`` (needs ``theta``) or ``"tweedie"`` (needs ``power``). link : str, optional ``"identity"``, ``"log"``, ``"logit"``, ``"probit"``, ``"cloglog"``, ``"inverse"``, ``"inverse_squared"`` or ``"power"`` (needs ``link_power``); the family's canonical link by default. dispersion : str or float, optional ``"pearson"``, ``"deviance"`` or a fixed value. By default 1 for the Poisson, binomial and negative binomial and Pearson's estimate otherwise; ``"pearson"`` with the Poisson is the over-dispersed (quasi-) Poisson, which accepts negative responses as long as the fitted means stay positive. theta : float, optional power : float, optional link_power : float, optional Examples -------- >>> from prospicio.models import Design, Glm >>> d = Design([[1.0] * 4, [0.0, 0.0, 1.0, 1.0]], ["(Intercept)", "young"], ... offset=[0.0, 0.0, 0.0, 0.0]) >>> fit = Glm("poisson", "log").fit(d, [1.0, 3.0, 4.0, 6.0]) >>> round(fit.coefficients[1], 10) == round(__import__("math").log(5 / 2), 10) True GlmFit A fitted GLM, from ``Glm.fit``. BayesGlm(family, link=None, prior_sd=2.5, intercept_sd=10.0, dispersion=None, dispersion_scale=10.0, chains=4, tune=1000, draws=1000, seed=0, target_accept=0.8, max_depth=10, theta=None, power=None, link_power=None) A Bayesian GLM sampled with NUTS (nuts-rs, the Rust core of nutpie). Normal priors with mean 0 on the coefficients: standard deviation ``intercept_sd`` for an all-ones column, ``prior_sd`` for the others (on the link scale; standardize covariates). For the Gaussian, gamma and inverse Gaussian the dispersion is sampled too, with a half-normal prior of scale ``dispersion_scale``, unless ``dispersion`` fixes it. Chains run in parallel, start near the maximum-likelihood fit, and replay exactly from ``seed``. Parameters ---------- family : str link : str, optional prior_sd : float, default 2.5 intercept_sd : float, default 10.0 dispersion : float, optional A fixed dispersion; 1 by default for the Poisson, binomial and negative binomial. A Tweedie needs one. dispersion_scale : float, default 10.0 chains, tune, draws : int, default 4, 1000, 1000 seed : int, default 0 target_accept : float, default 0.8 max_depth : int, default 10 theta, power, link_power : float, optional Examples -------- >>> from prospicio.models import BayesGlm, Design >>> x = [(i % 4) - 1.5 for i in range(40)] >>> y = [[1.0, 2.0, 3.0, 5.0][i % 4] for i in range(40)] >>> d = Design([[1.0] * 40, x], ["(Intercept)", "x"]) >>> fit = BayesGlm("poisson", chains=2, tune=300, draws=300).fit(d, y) >>> all(s["rhat"] < 1.05 for s in fit.summary()) True BayesGlmFit A sampled Bayesian GLM, from ``BayesGlm.fit``. ElasticNet(family, link=None, alpha=1.0, lam=0.0, standardize=True, penalty_factor=None, theta=None, power=None, link_power=None) An elastic-net GLM: the lasso (``alpha=1``), ridge (``alpha=0``) and everything between, minimizing glmnet's objective ``sum(w * d) / (2 * sum(w)) + lam * sum(pf * ((1 - alpha) / 2 * b**2 + alpha * |b|))`` over coefficients ``b`` of standardized columns. The design's first all-ones column is the unpenalized intercept; coefficients are reported on the design's scale. Parameters ---------- family : str As ``Glm``. link : str, optional As ``Glm``; the canonical link by default. alpha : float, default 1.0 Mixing between ridge (0) and the lasso (1). lam : float, default 0.0 Penalty strength (``lambda`` in glmnet). standardize : bool, default True Penalize the coefficients of columns scaled to unit standard deviation. penalty_factor : list of float, optional One factor per design column (the intercept's is ignored); 0 leaves a column unpenalized. theta : float, optional power : float, optional link_power : float, optional Examples -------- >>> from prospicio.models import Design, ElasticNet >>> x = [float(i) for i in range(6)] >>> d = Design([[1.0] * 6, x, [1.0, 0.0] * 3], ["(Intercept)", "x1", "x2"]) >>> y = [1.0, 3.1, 4.9, 7.2, 9.0, 10.8] >>> net = ElasticNet("gaussian", alpha=1.0) >>> top = net.lambda_max(d, y) >>> net.with_lam(1.01 * top).fit(d, y).coefficients[1:] [0.0, 0.0] ElasticNetFit A fitted elastic net, from ``ElasticNet.fit`` or ``ElasticNet.path``. CvPath Cross-validated scores along an elastic-net path, from ``ElasticNet.cross_validate``. Gam(glm, smooths, smoothing=None) A generalized additive model: a ``Glm`` plus P-spline smooths of numeric design columns, with smoothing chosen by GCV or UBRE. Parameters ---------- glm : Glm Family, link and dispersion. smooths : list of str or (str, int) The design columns to smooth, optionally with the number of basis functions (10 by default). smoothing : str or list of float, default "auto" ``"auto"`` (UBRE for a fixed dispersion, GCV otherwise), ``"gcv"``, ``"ubre"``, or fixed smoothing parameters, one per smooth. Examples -------- >>> import math >>> from prospicio.models import Design, Gam, Glm >>> x = [i / 99 for i in range(100)] >>> y = [math.sin(6 * v) for v in x] >>> d = Design([[1.0] * 100, x], ["(Intercept)", "x"]) >>> fit = Gam(Glm("gaussian"), ["x"]).fit(d, y) >>> abs(fit.predict(d)[50] - y[50]) < 0.01 True GamFit A fitted GAM, from ``Gam.fit``. deviance(family, y, mu, weights=None, theta=None, power=None) Deviance ``sum w d(y, mu)`` of a family. Parameters ---------- family : str y : list of float mu : list of float weights : list of float, optional theta : float, optional power : float, optional Returns ------- float gini(y, pred, exposure=None) Gini index of the ordered Lorenz curve. Parameters ---------- y : list of float pred : list of float exposure : list of float, optional Returns ------- float lift(y, pred, exposure=None, bands=10) Lift table: rows sorted by predicted rate, cut into bands of about equal exposure. Parameters ---------- y : list of float pred : list of float exposure : list of float, optional bands : int, default 10 Returns ------- list of dict ``exposure``, ``expected`` and ``actual`` per band. crps(draws, y) Continuous ranked probability score of equally likely draws for an outcome; lower is better. Parameters ---------- draws : list of float y : float Returns ------- float log_score(family, y, mu, dispersion=1.0, weights=None, theta=None, power=None) Mean log score ``-(1/n) sum log f(y_i)`` of the outcomes under the family's predictive distribution; lower is better. Parameters ---------- family : str y : list of float mu : list of float dispersion : float, default 1.0 weights : list of float, optional theta : float, optional power : float, optional Returns ------- float Examples -------- >>> from prospicio.models import log_score >>> round(log_score("poisson", [0.0], [1.0]), 12) 1.0 pit(family, y, mu, dispersion=1.0, weights=None, seed=0, theta=None, power=None) Probability integral transform of each outcome under the family's predictive distribution, randomized where it has atoms (counts, a Tweedie's zero); uniform when the model is calibrated. Parameters ---------- family : str y : list of float mu : list of float dispersion : float, default 1.0 weights : list of float, optional seed : int, default 0 Seeds the randomization. theta : float, optional power : float, optional Returns ------- list of float pit_from_draws(draws, y, u=0.5) The PIT of ``y`` under the empirical distribution of ``draws``, randomized over ties by ``u``. Parameters ---------- draws : list of float y : float u : float, default 0.5 Returns ------- float pit_histogram(pit, bins=10) Counts of PIT values in ``bins`` equal-width bins of ``[0, 1]``. Parameters ---------- pit : list of float bins : int, default 10 Returns ------- list of int ks_uniform(values) Kolmogorov-Smirnov distance from the uniform on ``[0, 1]``; about ``1.36 / sqrt(n)`` or less 95% of the time under uniformity. Parameters ---------- values : list of float Returns ------- float k_fold(n, k, seed) ``k``-fold splits of ``n`` rows, shuffled with ``seed``. Parameters ---------- n : int k : int seed : int Returns ------- list of (list of int, list of int) ``(train, test)`` row indices per fold. group_k_fold(groups, k, seed) Grouped ``k``-fold splits: each group's rows stay in one fold. Parameters ---------- groups : list of str k : int seed : int Returns ------- list of (list of int, list of int) time_ordered(periods, n_test) Time-ordered splits: for each of the last ``n_test`` periods, train on earlier periods and test on that one (for a triangle, the calendar diagonal backtest). Parameters ---------- periods : list of int n_test : int Returns ------- list of (list of int, list of int) mcmc_diagnostics(chains) MCMC diagnostics of chains of draws (Vehtari et al. 2021, as R's ``posterior``): rank-normalized split R-hat, bulk and tail effective sample sizes, the effective sample size of the mean and its Monte Carlo standard error. Parameters ---------- chains : list of list of float Equal-length chains, at least 4 draws each. Returns ------- dict ``rhat``, ``ess_bulk``, ``ess_tail``, ``ess_mean``, ``mcse_mean``. Examples -------- >>> from prospicio.models import mcmc_diagnostics >>> a = [float((i * 37) % 101) for i in range(400)] >>> b = [float((i * 53 + 7) % 101) for i in range(400)] >>> mcmc_diagnostics([a, b])["rhat"] < 1.01 True elpd_loo(log_lik, r_eff=None) Leave-one-out cross-validation by Pareto-smoothed importance sampling (PSIS-LOO), from one fit's pointwise log-likelihood draws. Matches the R package ``loo``. Parameters ---------- log_lik : list of list of float One row per posterior draw, one column per observation: ``log p(y_i | theta_s)``. r_eff : list of float, optional Relative efficiency of the draws per observation (1 for independent draws). Returns ------- Elpd With ``pareto_k`` and ``k_threshold``. stacking_weights(lpd) Stacking weights from pointwise held-out log predictive densities (Yao et al., 2018): the weights on the simplex that maximize the log score of the mixture of the models' predictive distributions. Works for any model: pass PSIS-LOO pointwise values (``Elpd.pointwise``) for a Bayesian fit, or cross-validated log densities for any other. A model that adds nothing gets weight exactly 0. Parameters ---------- lpd : list of list of float One list per model, each with one log density per observation. Returns ------- list of float One weight per model, summing to 1. Examples -------- >>> from prospicio.models import stacking_weights >>> w = stacking_weights([[-0.1, -0.1, -3.0, -3.0], [-3.0, -3.0, -0.1, -0.1]]) >>> [round(x, 9) for x in w] [0.5, 0.5] pseudo_bma_weights(lpd, bootstrap=True, n_draws=1000, seed=0) Pseudo-BMA weights, ``w_k`` proportional to ``exp(elpd_k)``; with ``bootstrap=True``, pseudo-BMA+ weights averaged over Bayesian-bootstrap replicates of the observations, which keeps a model that is only slightly better from taking all the weight. Parameters ---------- lpd : list of list of float One list per model, as for ``stacking_weights``. bootstrap : bool, default True n_draws : int, default 1000 Bootstrap replicates. seed : int, default 0 Replicate ``b`` uses stream ``b`` of ``seed``. Returns ------- list of float BayesStacking(concentration=None, chains=4, tune=1000, draws=1000, seed=0) Bayesian stacking: a posterior for the stacking weights, with a Dirichlet prior, sampled by NUTS from pointwise held-out log densities (Yao et al., 2018). ``stacking_weights`` gives the optimum alone. Parameters ---------- concentration : list of float, optional Dirichlet concentration, one per model (default 1, uniform). chains, tune, draws : int, default 4, 1000, 1000 seed : int, default 0 HierarchicalStacking(discrete=0, alpha_loc=0.0, alpha_scale=1.0, beta_loc=0.0, beta_scale=1.0, partial_pooling=False, tau_mu_global=1.0, tau_mu_discrete=1.0, tau_mu_continuous=1.0, tau_sigma_discrete=1.0, tau_sigma_continuous=1.0, adaptive=None, chains=4, tune=1000, draws=1000, seed=0) Hierarchical stacking (Yao, Pirš, Vehtari and Gelman, 2022): model weights that vary with covariates, ``w = softmax(alpha + B x)`` against the last model as reference, so a model can be trusted in one part of the portfolio and not another. The priors are those of BayesBlend's ``HierarchicalBayesStacking``; sampled by NUTS. Scale continuous covariates (BayesBlend divides by twice the standard deviation) and dummy-code discrete ones before fitting, with the dummies first. With ``partial_pooling``, each model's slopes on the discrete covariates, and separately on the continuous ones, are drawn around a model-level mean, itself drawn around a global mean. A scale of 0 removes a level: ``tau_mu_global=0`` fixes the global mean at 0, ``tau_mu_*=0`` pools completely and ``tau_sigma_*=0`` sets every slope to its model's mean. BayesBlend warns that pooling needs at least three covariates. ``adaptive`` multiplies the prior scales by ``N**lambda`` with ``lambda ~ Exponential(adaptive)``, weakening them as the data grow. Parameters ---------- discrete : int, default 0 Number of leading covariates that are dummy codes. alpha_loc, alpha_scale : float, default 0.0, 1.0 beta_loc, beta_scale : float, default 0.0, 1.0 Slope prior without pooling. partial_pooling : bool, default False tau_mu_global, tau_mu_discrete, tau_mu_continuous : float, default 1.0 tau_sigma_discrete, tau_sigma_continuous : float, default 1.0 Pooling scales, BayesBlend's defaults. adaptive : float, optional Rate of the exponential prior on ``lambda`` (BayesBlend uses 4). chains, tune, draws : int, default 4, 1000, 1000 seed : int, default 0 Examples -------- >>> from prospicio.models import HierarchicalStacking >>> x = [i / 99 - 0.5 for i in range(100)] >>> a = [-0.5 if v < 0 else -2.0 for v in x] >>> b = [-2.0 if v < 0 else -0.5 for v in x] >>> fit = HierarchicalStacking(chains=2, tune=300, draws=300).fit([a, b], [x]) >>> w = fit.weights([[-0.4, 0.4]]) >>> w[0][0] > 0.7 and w[1][0] < 0.3 True StackingFit Posterior stacking weights, from ``BayesStacking.fit`` or ``HierarchicalStacking.fit``. elpd_waic(log_lik) WAIC from pointwise log-likelihood draws: ``lppd`` less the variance of each observation's log-likelihood. Parameters ---------- log_lik : list of list of float One row per posterior draw, one column per observation. Returns ------- Elpd lppd(log_lik) In-sample log pointwise predictive density, ``sum_i log(mean_s p(y_i | theta_s))``. Parameters ---------- log_lik : list of list of float One row per posterior draw, one column per observation. Returns ------- float Examples -------- >>> import math >>> from prospicio.models import lppd >>> round(lppd([[math.log(0.5)], [math.log(0.25)]]), 12) == round(math.log(0.375), 12) True Elpd An ELPD estimate from ``elpd_loo`` or ``elpd_waic``. deviance_score(family, theta=None, power=None) A score for ``cross_validate``: the family's mean deviance on the test rows, ``sum(w * d) / sum(w)`` with the test design's weights. Parameters ---------- family : str As ``Glm``. theta : float, optional power : float, optional Returns ------- callable ``score(y_test, predicted, test_design) -> float``. Examples -------- >>> from prospicio.models import Design, deviance_score >>> d = Design([[1.0, 1.0]], ["(Intercept)"]) >>> deviance_score("gaussian")([1.0, 3.0], [2.0, 2.0], d) 1.0 cross_validate(model, design, y, splits, score, n_jobs=None) Fits ``model`` on each split's training rows and scores it on the test rows. Folds run on threads; the Rust fits release the GIL, so they run in parallel. Parameters ---------- model : object Anything with ``fit(design, y)`` returning an object with ``predict(design)``: ``Glm``, ``ElasticNet``, ``Gam``. design : Design y : list of float splits : list of (list of int, list of int) ``(train, test)`` rows, as ``k_fold`` returns. score : callable ``score(y_test, predicted, test_design) -> float``, a loss (lower is better), such as ``deviance_score("poisson")``. n_jobs : int, optional Threads; one per split by default. Returns ------- list of float One score per split. Examples -------- >>> from prospicio.models import Design, Glm, cross_validate, deviance_score, k_fold >>> d = Design([[1.0] * 6, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]], ["(Intercept)", "x"]) >>> y = [1.0, 2.9, 5.1, 7.0, 8.9, 11.2] >>> scores = cross_validate(Glm("gaussian"), d, y, k_fold(6, 3, 1), deviance_score("gaussian")) >>> len(scores) 3 grid_search(candidates, make, design, y, splits, score, n_jobs=None) Scores every candidate by ``cross_validate`` on the same splits and picks the lowest mean score. Parameters ---------- candidates : list Hyperparameter values, in any form ``make`` accepts. make : callable ``make(candidate) -> model``. design : Design y : list of float splits : list of (list of int, list of int) score : callable As ``cross_validate``. n_jobs : int, optional Returns ------- SearchResult Examples -------- >>> from prospicio.models import Design, ElasticNet, deviance_score, grid_search, k_fold >>> x = [i / 10 for i in range(40)] >>> d = Design([[1.0] * 40, x], ["(Intercept)", "x"]) >>> y = [1.0 + 2.0 * v for v in x] >>> found = grid_search([0.0, 1.0, 10.0], lambda lam: ElasticNet("gaussian", lam=lam), ... d, y, k_fold(40, 4, 1), deviance_score("gaussian")) >>> found.best_candidate 0.0 random_search(n, seed, draw, make, design, y, splits, score, n_jobs=None) Draws ``n`` candidates with ``draw(rng)`` from ``random.Random(seed)`` and scores them as ``grid_search`` does: with several hyperparameters, random candidates cover each range better than a grid of the same size. Parameters ---------- n : int seed : int draw : callable ``draw(rng) -> candidate``, for example ``lambda rng: log_uniform(rng, 1e-4, 1.0)``. make, design, y, splits, score, n_jobs As ``grid_search``. Returns ------- SearchResult log_uniform(rng, low, high) A draw log-uniform between ``low`` and ``high``, for learning rates and penalties. Parameters ---------- rng : random.Random low, high : float Positive bounds. Returns ------- float SearchResult(scores: list, best: int) -> None The result of ``grid_search`` or ``random_search``. Attributes ---------- scores : list of (object, float) Each candidate with its mean score over the splits. best : int Index of the lowest mean score. compare(models, design, y, splits, scores, n_jobs=None) Fits every model on each split's training rows and scores its test predictions on every metric: one table across engines. Parameters ---------- models : dict Name to model; each model has ``fit(design, y)`` returning an object with ``predict(design)``, as for ``cross_validate``. design : Design y : list of float splits : list of (list of int, list of int) scores : dict Name to ``score(y_test, predicted, test_design) -> float``, losses (lower is better). n_jobs : int, optional Threads; one per split by default. Returns ------- Comparison Examples -------- >>> from prospicio.models import Design, ElasticNet, Glm, compare, deviance_score, k_fold >>> x = [i / 10 for i in range(40)] >>> d = Design([[1.0] * 40, x], ["(Intercept)", "x"]) >>> y = [1.0 + 2.0 * v + (0.3 if i % 3 else -0.6) for i, v in enumerate(x)] >>> c = compare({"glm": Glm("gaussian"), "ridge": ElasticNet("gaussian", alpha=0.0, lam=5.0)}, ... d, y, k_fold(40, 4, 1), {"deviance": deviance_score("gaussian")}) >>> c.best("deviance") 'glm' >>> [row["model"] for row in c.table()] ['glm', 'ridge'] Comparison(models: list, metrics: list, split_scores: dict) -> None The result of ``compare``: every model's score on every metric and split. Attributes ---------- models : list of str metrics : list of str split_scores : dict ``split_scores[(model, metric)]`` is the list of per-split scores. actual_vs_expected(periods, family, y, mu, weights=None, dispersion=1.0, theta=None, power=None) Actual against expected by period for a stored model's predictions on new data, with each period's z-score under the model and a test for drift. ``A = sum(w * y)`` and ``E = sum(w * mu)`` per period; the z-score is ``(A - E) / sqrt(dispersion * sum(w * V(mu)))`` with the family's variance function ``V``, about standard normal while the model holds. ``trend`` is the slope of ``A / E - 1`` per period step (periods in sorted order), weighted by each period's precision. Parameters ---------- periods : list of int or str One period label per row. family : str y : list of float Actuals. mu : list of float The model's predicted means. weights : list of float, optional Prior weights as fitted (exposure for a rate); leave out for counts with exposure in the offset. dispersion : float, default 1.0 theta, power : float, optional Negative binomial ``theta``, Tweedie ``power``. Returns ------- dict ``periods`` (a list of dicts with ``period``, ``n``, ``weight``, ``actual``, ``expected``, ``ratio``, ``std_dev`` and ``z``), ``total`` (the same without ``period``), ``trend``, ``trend_std_error`` and ``trend_z``. Examples -------- >>> from prospicio.models import actual_vs_expected >>> m = actual_vs_expected([2023, 2023, 2024, 2024], "poisson", ... [1.0, 3.0, 2.0, 6.0], [2.0, 2.0, 2.0, 2.0]) >>> m["periods"][1]["ratio"], m["periods"][1]["z"] (2.0, 2.0) simulate_from_means(family, means, n_sims, seed, dispersion=None, weights=None, theta=None, power=None) Joint predictive draws from fitted means, for engines that give only a mean per row (the boosting adapters): the family adds process noise and several mean vectors (bootstrap refits) add parameter uncertainty. Simulation ``i`` uses stream ``i`` of ``seed``: it picks one mean vector uniformly, then draws each row's response from the family with that mean, the row's dispersion and the row's weight. Components are keyed ``row = 0, 1, ...``, as ``GlmFit.predict_distribution`` keys them. Parameters ---------- family : str As in ``Glm``. means : list of list of float One or more mean vectors, one value per row each. n_sims : int seed : int dispersion : float or list of float, optional One value for every row, or one per row (from a dispersion model); 1 by default. weights : list of float, optional Prior weights; 1 by default. theta, power : float, optional Negative binomial ``theta``, Tweedie ``power``. Returns ------- PredictiveDistribution Examples -------- >>> from prospicio.models import simulate_from_means >>> pd = simulate_from_means("poisson", [[0.1, 0.4]], 20_000, 7) >>> round(pd.mean(), 1) 0.5 ## Gradient boosting LightGBM and XGBoost behind the model protocol, with joint predictive draws. Booster(family='poisson', engine='lightgbm', power=None, n_rounds=200, learning_rate=0.05, params=None, n_boot=0, seed=0, alpha=None, dispersion_model=False) Gradient-boosted trees for a response family, through LightGBM or XGBoost. Parameters ---------- family : {"poisson", "gamma", "tweedie", "gaussian", "quantile"}, default "poisson" The engine's objective; a log link for the first three. ``"quantile"`` fits the ``alpha`` quantile of the response (no link, no offset). engine : {"lightgbm", "xgboost"}, default "lightgbm" power : float, optional Tweedie variance power in (1, 2); required for ``"tweedie"``. alpha : float, optional Quantile level in (0, 1); required for ``"quantile"``. dispersion_model : bool, default False Fit a dispersion per row (gamma, Tweedie and Gaussian families). n_rounds : int, default 200 Boosting rounds. learning_rate : float, default 0.05 params : dict, optional Further engine parameters (``num_leaves``, ``max_depth``, ``monotone_constraints``, ...), passed as they are. n_boot : int, default 0 Bootstrap refits for parameter uncertainty in ``predict_distribution``; 0 gives process noise only. seed : int, default 0 Seeds the engine and the bootstrap resamples. Examples -------- >>> from prospicio.boosting import Booster >>> from prospicio.models import Design >>> x = [i % 10 / 10 for i in range(200)] >>> d = Design([x], ["x"], offset=[0.0] * 200) >>> y = [float(i % 3 == 0) + v for i, v in enumerate(x)] >>> fit = Booster("poisson", n_rounds=20).fit(d, y) >>> len(fit.predict(d)) 200 BoosterFit(spec, names, model, boots, dispersion, base) A fitted ``Booster``: means and joint predictive draws. Attributes ---------- model The engine's booster (``lightgbm.Booster`` or ``xgboost.Booster``). base : float The constant added to the offset so the trees start from the mean. dispersion : float 1 for the Poisson; otherwise Pearson's estimate on the training rows, ``sum w (y - mu)^2 / V(mu) / n`` (no degrees-of-freedom correction: trees have no fixed parameter count). dispersion_fit With ``dispersion_model=True``, the dispersion booster and its starting log level; otherwise ``None``. ## Pricing The collective model, layer rating, reinsurance tower matching, and risk-loaded prices from simulated losses. CollectiveModel(frequency, severity) The collective risk model: a claim count and a severity, with layer moments in closed form. For the layer ``limit`` xs ``attachment`` applied to each loss, with ``Y`` the loss to the layer from one claim, the aggregate has mean ``E[N] E[Y]`` and variance ``E[N] Var[Y] + Var[N] E[Y]**2``. Parameters ---------- frequency : Poisson, NegativeBinomial or Binomial severity : Lognormal, Grid, Pareto, PiecewisePareto, LogAffinePareto or GeneralizedPareto Examples -------- >>> from prospicio.distributions import Pareto, claim_count >>> from prospicio.pricing import CollectiveModel >>> m = CollectiveModel(claim_count(2.0, 1.5), Pareto(1e6, 2.0)) >>> round(m.layer_mean(4e6, 1e6)) 1600000 >>> m.excess_frequency(2e6) 0.5 ilf(severity, limit, basic_limit) Increased limit factor ``LEV(limit) / LEV(basic_limit)``. Parameters ---------- severity : a severity limit : float basic_limit : float Returns ------- float Examples -------- >>> from prospicio.distributions import Pareto >>> from prospicio.pricing import ilf >>> round(ilf(Pareto(100.0, 2.0), 1000.0, 200.0), 12) 1.266666666667 loss_elimination_ratio(severity, deductible) Loss elimination ratio of a deductible, ``LEV(deductible) / E[X]``. Parameters ---------- severity : a severity deductible : float Returns ------- float pareto_extrapolation(from_, to, alpha, truncation=None) Expected loss of layer ``to`` per unit of expected loss of layer ``from_``, under a Pareto with this alpha (and truncation). Parameters ---------- from_ : tuple of float ``(limit, attachment)``. to : tuple of float ``(limit, attachment)``. alpha : float truncation : float, optional Returns ------- float Examples -------- >>> from prospicio.pricing import pareto_extrapolation >>> round(pareto_extrapolation((1e6, 1e6), (2e6, 2e6), 2.0), 12) 0.5 alpha_between_layers(a, b, truncation=None) The Pareto alpha at which two layers have the given expected losses. Parameters ---------- a : tuple of float ``(limit, attachment, expected_loss)``. b : tuple of float ``(limit, attachment, expected_loss)``; one layer must lie above the other. truncation : float, optional Returns ------- float alpha_between_frequency_and_layer(threshold, frequency, limit, attachment, expected_loss, truncation=None) The Pareto alpha at which ``frequency`` losses a year above ``threshold`` give the layer an expected loss of ``expected_loss``. Parameters ---------- threshold : float frequency : float limit : float attachment : float expected_loss : float truncation : float, optional Returns ------- float alpha_between_frequencies(threshold_1, frequency_1, threshold_2, frequency_2, truncation=None) The Pareto alpha between two excess frequencies. Parameters ---------- threshold_1 : float frequency_1 : float threshold_2 : float frequency_2 : float truncation : float, optional Returns ------- float Examples -------- >>> from prospicio.pricing import alpha_between_frequencies >>> round(alpha_between_frequencies(1e6, 4.0, 2e6, 1.0), 12) 2.0 match_tower(attachments, layer_losses, frequencies=None, rule='minimize') Matches a tower of contiguous layers, the last unlimited, with one frequency and a piecewise Pareto severity (Riegel 2018). Parameters ---------- attachments : list of float Increasing attachment points; layer ``i`` runs to the next one, the last is unlimited. layer_losses : list of float Expected loss a year of each layer. frequencies : list of float or None, optional Expected losses a year above each attachment point; ``None`` (or a ``None`` entry) to derive them. rule : {"minimize", "midpoint"}, default "minimize" How the free threshold inside each layer is chosen. Returns ------- TowerModel Examples -------- >>> from prospicio.pricing import match_tower >>> m = match_tower([1000.0, 1500.0, 2000.0], [100.0, 90.0, 120.0], [0.25, None, None]) >>> round(m.layer_loss(500.0, 1500.0), 9) 90.0 fit_pml_curve(return_periods, amounts, tail_alpha=2.0, truncation=None) The model through the points of a PML curve: ``amounts[j]`` is exceeded once in ``return_periods[j]`` years. Parameters ---------- return_periods : list of float amounts : list of float tail_alpha : float, default 2.0 Alpha above the largest amount. truncation : float, optional Truncation of the last piece. Returns ------- TowerModel fit_references(layers=Ellipsis, frequencies=Ellipsis, default_alpha=2.0, rule='minimize') A model that reproduces every reference: expected layer losses (which may overlap or leave gaps) and excess frequencies. Parameters ---------- layers : list of tuple of float, optional ``(limit, attachment, expected_loss)`` per layer. frequencies : list of tuple of float, optional ``(threshold, frequency)`` per excess frequency. default_alpha : float, default 2.0 Alpha above the highest point, unless an unlimited layer sets it. rule : {"minimize", "midpoint"}, default "minimize" Returns ------- TowerModel Examples -------- >>> from prospicio.pricing import fit_references >>> m = fit_references([(1000.0, 1000.0, 150.0), (3000.0, 1500.0, 160.0)], [(1000.0, 0.3)]) >>> round(m.layer_loss(3000.0, 1500.0), 6) 160.0 TowerModel A frequency and a piecewise Pareto severity that reproduce a tower, a PML curve or a set of references. price(losses, assets, *, cost_of_capital=None, distortion=None) Risk-loaded price of a cover from its simulated losses. The assets backing the loss are a distortion risk measure of it. The premium is either a pricing distortion of the loss, or set by a constant cost of capital ``r`` on the capital ``a - P``, which gives ``P = (E[X] + r a) / (1 + r)``. Parameters ---------- losses : Sampled or PredictiveDistribution Loss draws; for a ``PredictiveDistribution``, its total. assets : Distortion The measure that sets the assets, for example ``Distortion.tvar(0.99)``. cost_of_capital : float, optional Positive rate. Give this or ``distortion``. distortion : Distortion, optional Pricing distortion; it must load less than ``assets``. Returns ------- Price Raises ------ ValueError Unless exactly one rule is given, or if the premium exceeds the assets. Examples -------- >>> from prospicio.distributions import Sampled >>> from prospicio.pricing import price >>> from prospicio.risk import Distortion >>> p = price(Sampled([0.0, 0.0, 2.0, 6.0]), Distortion.tvar(0.5), cost_of_capital=0.25) >>> p.premium, p.capital (2.4, 1.6) price_portfolio(pd, assets, *, cost_of_capital=None, distortion=None) Prices a portfolio and allocates the price to its components. Premium and assets are each allocated by co-measure (the natural allocation): component prices add up to the portfolio's, and a component that diversifies the portfolio is priced below its standalone price. With a cost of capital, every component earns the rate on its allocated capital. Parameters ---------- pd : PredictiveDistribution Components that add up to the portfolio: segments or covers, not gross, ceded and net side by side. assets : Distortion cost_of_capital : float, optional distortion : Distortion, optional Exactly one of ``cost_of_capital`` and ``distortion``, as in ``price``. Returns ------- PortfolioPrice Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.pricing import price_portfolio >>> from prospicio.risk import Distortion >>> pd = PredictiveDistribution(["cover"], [("a",), ("b",)], ... [[0.0, 2.0], [1.0, 1.0], [4.0, 0.0], [8.0, 0.0]]) >>> p = price_portfolio(pd, Distortion.tvar(0.5), cost_of_capital=0.1) >>> [round(c.premium, 6) for c in p.allocated] [3.5, 0.681818] >>> p.allocated[1].margin < 0 # the second cover hedges the first True Price The risk-loaded price of a cover, or of one component's share of a portfolio: expected loss, premium and the assets backing the loss. Returned by ``price`` and ``price_portfolio``. PortfolioPrice Prices of a portfolio's components and of the portfolio as a whole. Returned by ``price_portfolio``. Mbbefd(b, g) The MBBEFD exposure curve and destruction-rate distribution (Bernegger, 1997), with ``b >= 0`` and ``g >= 1``; ``1/g`` is the probability of a total loss. ``G(x)`` is the share of a risk's expected loss below the fraction ``x`` of its maximum possible loss (MPL). ``Mbbefd.swiss_re(c)`` gives Bernegger's one-parameter family: ``c = 1.5, 2, 3, 4`` are the Swiss Re curves and ``c = 5`` the Lloyd's curve. Parameters ---------- b : float g : float Examples -------- >>> from prospicio.pricing import Mbbefd >>> c3 = Mbbefd.swiss_re(3.0) >>> top = c3.layer_share(5e6, 5e6, 10e6) >>> bottom = c3.layer_share(5e6, 0.0, 10e6) >>> round(top + bottom, 12), top < bottom (1.0, True) TabulatedCurve(x, g) A tabulated exposure curve: points ``(x, G(x))`` from ``(0, 0)`` to ``(1, 1)``, interpolated linearly, as published curves are given (Salzmann's homeowners scale, Ludwig's curves, ISO PSOLD tables, a reinsurer's own). The table must be concave (its slopes never increase). Its destruction rate is discrete: the points' ``x`` with probabilities from the drops in slope, and a total loss with probability last slope over first. Its mean rate is the first chord's, ``x1 / G(x1)``, so a table needs fine first points for the expected loss to be right. Parameters ---------- x : list of float Increasing from 0 to 1. g : list of float ``G(x)``, from 0 to 1. Raises ------ ValueError If the points do not run from ``(0, 0)`` to ``(1, 1)``, or are not increasing and concave. Examples -------- >>> from prospicio.pricing import TabulatedCurve >>> t = TabulatedCurve([0.0, 0.1, 0.5, 1.0], [0.0, 0.4, 0.8, 1.0]) >>> round(t.curve([0.3])[0], 12), t.mean_rate() (0.6, 0.25) RiskProfile(sums_insured, risks, curves, expected_losses=None, premiums=None, loss_ratio=None, lower=None, upper=None, spread='uniform') A risk profile for property per-risk business: bands of sum insured, each with an expected loss (given, or premium times a loss ratio) and its own exposure curve. Each band's representative risk has sum insured ``SI`` (its total sum insured over its number of risks, say), taken as its MPL. The band expects ``EL / (SI * curve.mean_rate)`` losses a year; each simulated loss is the band's ``SI`` times a destruction rate from the band's curve, and carries that ``SI``, so a surplus treaty (``Layer.surplus``) and the per-risk excess of loss it inures to apply to the events. The exposure-rated expectations (``expected_layer_loss``, ``expected_surplus_loss``) check the simulation. Parameters ---------- sums_insured : list of float One per band. risks : list of float Number of risks per band (for reference). curves : Mbbefd or TabulatedCurve, or a list of them One curve for every band, or one per band. expected_losses : list of float, optional Expected annual loss per band. Give this, or ``premiums``. premiums : list of float, optional Premium per band, with ``loss_ratio``. loss_ratio : float or list of float, optional Expected loss ratio, one for all bands or one per band. lower, upper : list of float or None, optional Bounds of each band's sums insured (``None`` for a band without). A band with bounds spreads its risks' sums insured uniformly between them: its mean ``SI`` is ``(lower + upper) / 2`` (in place of ``sums_insured``), each simulated loss draws its own ``SI`` between the bounds, and the exposure-rated expectations average over the band, weighted by sum insured. spread : {"uniform", "tilted"}, default "uniform" How a band with bounds spreads its sums insured. ``"tilted"`` keeps ``sums_insured`` as the mean, so the spread matches both the bounds and the band's total sum insured: the density ``∝ exp(θ s)`` on the bounds with ``θ`` solved for that mean (``sums_insured`` must lie strictly between the bounds). Examples -------- >>> from prospicio.pricing import Mbbefd, RiskProfile >>> p = RiskProfile([1e6, 10e6], [800, 50], Mbbefd.swiss_re(3.0), ... premiums=[2e6, 1e6], loss_ratio=0.6) >>> round(p.expected_loss()) 1800000 >>> events = p.simulate(1000, 7) >>> events.has_sums_insured True severity_exposure_curve(severity, mpl, x) The exposure curve of a severity capped at the maximum possible loss ``mpl``: ``G(x) = LEV(x mpl) / LEV(mpl)`` at each ``x``. Parameters ---------- severity : a severity mpl : float x : list of float Returns ------- list of float Examples -------- >>> from prospicio.distributions import Pareto >>> from prospicio.pricing import severity_exposure_curve >>> g = severity_exposure_curve(Pareto(1e5, 1.5), 1e7, [0.0, 0.5, 1.0]) >>> g[0], round(g[2], 12), g[1] > 0.5 (0.0, 1.0, True) ## Reserving Loss triangles, the chain ladder with tail factors, Mack's model with its one-year view, the expected-loss methods, the ODP bootstrap and Clark's growth curves. Triangle A loss triangle with four axes: index (segment), column (measure), origin and development age, in chainladder-python's order. Segments are named by key columns such as ``"lob"`` and ``"state"``: ``keys`` gives their names and ``index`` one label per segment. A triangle without keys has one segment, ``"Total"``. Build one from a long table with ``from_long`` or ``from_frame``. Ages are whole months from the start of the origin period, so age 12 on a 2021 accident year is valued at December 2021. Cells that were not observed are ``nan`` in ``values``; an observed zero stays zero. Examples -------- >>> from prospicio.reserving import Triangle >>> tri = Triangle.from_long( ... origin=[2020, 2020, 2021], ... development=[12, 24, 12], ... values={"paid": [100.0, 150.0, 110.0]}, ... ) >>> tri.shape (1, 1, 2, 2) >>> tri.origins, tri.development, tri.valuation (['2020', '2021'], [12, 24], datetime.date(2021, 12, 31)) >>> tri.values[0][0] [[100.0, 150.0], [110.0, nan]] ChainLadder(average='volume', sigma_interpolation='log-linear', tail=None) The chain-ladder method: each origin's latest value projected to ultimate with age-to-age factors estimated from the triangle and a tail. Parameters ---------- average : {"volume", "simple", "regression"}, default "volume" How link ratios are averaged into one factor per age: volume weighted, their mean, or least squares through the origin (Mack's ``alpha`` of 1, 0 and 2). sigma_interpolation : {"log-linear", "mack"}, default "log-linear" How a variance parameter with a single link ratio is filled in. tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional Development past the oldest age: a number is a constant factor from the oldest age to ultimate. No tail (a factor of 1) by default. Examples -------- >>> from prospicio.reserving import ChainLadder, Triangle >>> tri = Triangle.from_long([2020, 2020, 2021], [12, 24, 12], {"paid": [100.0, 150.0, 200.0]}) >>> fit = ChainLadder().fit(tri, "paid") >>> fit.ldf, fit.ultimate, fit.total_reserve ([1.5], [150.0, 300.0], 100.0) ChainLadderFit A fitted chain-ladder projection of every segment of a triangle column. Per-origin lists (``origins``, ``latest``, ``ultimate``, ``reserve``) run over the origins of each segment in turn, like the rows of ``to_frame()``, so a single-segment fit has one value per origin. Per-age lists (``ldf``, ``cdf``, ``sigma``, ``std_err``) and the tail need a single-segment fit; for several segments use ``development_frame()`` (per age), ``totals_frame()`` (``tail``, ``tail_sigma``, ``tail_std_err``) or ``segment(...)``. Examples -------- >>> from prospicio.reserving import ChainLadder, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021] * 2, ... [12, 24, 12] * 2, ... {"paid": [100.0, 150.0, 200.0, 10.0, 20.0, 30.0]}, ... keys={"lob": ["Auto"] * 3 + ["Home"] * 3}, ... ) >>> fit = ChainLadder().fit(tri, "paid") >>> fit.index, fit.reserve (['Auto', 'Home'], [0.0, 100.0, 0.0, 30.0]) >>> fit.segment(lob="Home").ldf [2.0] TailConstant(factor=1.0, decay=0.5, attachment_age=None) A given tail factor, as chainladder-python's ``TailConstant``. The factor applies from the attachment age to ultimate. Past the attachment it is spread over the following periods as ``1 + x * decay**k``, the last factor making up the difference; this shapes the factors past the attachment, not the factor to ultimate. An attachment before the oldest age replaces the estimated factors from there. Parameters ---------- factor : float, default 1.0 Factor from the attachment age to ultimate; finite and positive. decay : float, default 0.5 Share of each period's development kept in the next, from 0 to 1. attachment_age : int, optional Age in months the factor attaches at (the first age at or after it); the oldest age by default. An age at or before the youngest replaces every estimated factor (chainladder-python ignores such an attachment). Examples -------- >>> from prospicio.reserving import ChainLadder, TailConstant, Triangle >>> tri = Triangle.from_long([2020, 2020, 2021], [12, 24, 12], {"paid": [100.0, 150.0, 200.0]}) >>> fit = ChainLadder(tail=TailConstant(1.05)).fit(tri, "paid") >>> fit.tail, round(fit.ultimate[1], 6) (1.05, 315.0) TailCurve(curve='exponential', fit_period=Ellipsis, extrap_periods=100, attachment_age=None) A curve fitted to the estimated factors and extrapolated, as chainladder-python's ``TailCurve``. Factors above 1.00001 in the fit period are regressed by least squares: ``ln(f - 1)`` on the 1-based development index ``k`` (exponential) or on ``ln(k)`` (inverse power). The fitted curve replaces the factors from the attachment age on and runs ``extrap_periods`` periods past the oldest age. Parameters ---------- curve : {"exponential", "inverse_power"}, default "exponential" fit_period : tuple of (int or None, int or None), default (None, None) Ages in months whose factors enter the fit: from the last age at or before the first (inclusive) to the last age at or before the second (exclusive), as chainladder-python reads them; ``None`` is open-ended. extrap_periods : int, default 100 Number of periods past the oldest age the curve is extrapolated. attachment_age : int, optional Age in months the curve attaches at (the first age at or after it); the oldest age by default. Examples -------- >>> from prospicio.reserving import ChainLadder, TailCurve, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> fit = ChainLadder(tail=TailCurve()).fit(tri, "values") >>> 1.0 < fit.tail < 1.05 True TailBondy(earliest_age=None, attachment_age=None) The Bondy tail, as chainladder-python's ``TailBondy``. Each log factor from ``earliest_age`` on is taken as ``b`` times the one before it, ``b`` fitted by least squares. The fitted factors are ``f0 ** (b ** j)`` from the factor ``f0`` at ``earliest_age``, and those past the next one multiply to the last fitted factor raised to ``b / (1 - b)``. With the default ``earliest_age`` (the age of the last factor) ``b`` is 1/2 and the tail repeats the last factor. Parameters ---------- earliest_age : int, optional First age in months whose factor enters the fit (the last age at or before it, as chainladder-python reads it); the age of the last factor by default. attachment_age : int, optional The factor from this age (the last age at or before it) to the next is kept and the fitted ones replace those after it; the age of the last factor by default. Not before ``earliest_age``. Examples -------- >>> from prospicio.reserving import ChainLadder, TailBondy, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2020, 2021, 2021, 2022], ... [12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 110.0, 170.0, 120.0], ... ) >>> round(ChainLadder(tail=TailBondy()).fit(tri, "values").tail, 12) 1.1 TailLogLinear() R ChainLadder's ``tail = TRUE`` rule (its ``tailfactor`` function). When the third- and second-last factors multiply to more than 1.0001, ``ln(f - 1)`` is regressed on the development index over the factors above 1 and the next 100 extrapolated factors are multiplied; otherwise the tail is 1. A tail above 2 is reset to 1, as R does. Examples -------- >>> from prospicio.reserving import Mack, TailLogLinear, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> fit = Mack(tail=TailLogLinear()).fit(tri, "values") >>> fit.tail > 1.0 and fit.standard_error[0] > 0.0 True Mack(average='volume', sigma_interpolation='log-linear', tail=None, tail_sigma=None, tail_std_err=None) Mack's distribution-free chain ladder: the chain-ladder projection plus the standard error of each origin's reserve and of the total, split into process and parameter risk (Mack 1993, 1999). A tail other than 1 is one more development step, from the oldest age to ultimate, with its own sigma and standard error, as R ChainLadder's ``MackChainLadder(tail = ...)``; unless given, both are extrapolated log-linearly. Every origin, the oldest included, carries the tail's risk. A tail below 1 follows chainladder-python: it scales the ultimates and carries the risk read where a tail of 1.001 would be. R's ``MackChainLadder`` ignores a tail below 1 altogether. Parameters ---------- average : {"volume", "simple", "regression"}, default "volume" sigma_interpolation : {"log-linear", "mack"}, default "log-linear" tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional As ``ChainLadder``; no tail by default. tail_sigma : float, optional The tail's sigma (R's ``tail.sigma``); extrapolated if not given. Unused when the tail factor is 1. tail_std_err : float, optional The tail factor's standard error (R's ``tail.se``); extrapolated if not given. Unused when the tail factor is 1. Examples -------- >>> from prospicio.reserving import Mack, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> fit = Mack().fit(tri, "values") >>> fit.total_standard_error > 0 and fit.standard_error[0] == 0 True MackFit A fitted Mack model of every segment: the chain-ladder fields, plus standard errors of each origin's reserve and of each segment's total. Per-origin lists run over the origins of each segment in turn, like the rows of ``to_frame()``. Per-age lists, the tail and the totals' standard errors need a single-segment fit; for several segments use ``development_frame()``, ``totals_frame()`` (the totals' standard errors and the tail) or ``segment(...)``. ``total_ultimate`` and ``total_reserve`` sum over every segment. ClaimsDevelopmentResult Merz and Wüthrich's (2008) one-year view of a Mack fit: standard errors of the claims development result (CDR), the change in the chain-ladder ultimate over a calendar year, per origin and in total. The total includes the covariance between origins. Year ``k`` of the run-off is R ChainLadder's ``CDR(k)S.E.``; summed in square over the years, the run-off gives back Mack's standard error. Examples -------- >>> from prospicio.reserving import Mack, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> mack = Mack().fit(tri, "values") >>> cdr = mack.claims_development_result() >>> len(cdr.by_calendar_year), cdr.one_year_standard_error[0] (3, 0.0) >>> abs(cdr.total_run_off_standard_error - mack.total_standard_error) < 1e-9 True ExpectedLoss(apriori=1.0, average='volume', sigma_interpolation='log-linear', tail=None) The expected loss ratio method: each origin's ultimate is ``apriori`` times its exposure, whatever has been observed. The chain ladder is still fitted for the development pattern the fit reports. The exposure is a measure column of the same triangle (premium, say): each origin's latest observed cumulative value in the segment fitted. Parameters ---------- apriori : float, default 1.0 Expected loss ratio: the ultimate per unit of exposure; positive. average : {"volume", "simple", "regression"}, default "volume" How link ratios are averaged, as in ``ChainLadder``. sigma_interpolation : {"log-linear", "mack"}, default "log-linear" tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional As ``ChainLadder``; no tail by default. Examples -------- >>> from prospicio.reserving import ExpectedLoss, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021], [12, 24, 12], ... {"paid": [100.0, 150.0, 200.0], "premium": [250.0, 250.0, 400.0]}, ... ) >>> fit = ExpectedLoss(apriori=0.5).fit(tri, "paid", "premium") >>> fit.ultimate, fit.reserve ([125.0, 200.0], [-25.0, 0.0]) BornhuetterFerguson(apriori=1.0, average='volume', sigma_interpolation='log-linear', tail=None) The Bornhuetter–Ferguson method: each origin's latest value plus the expected loss ``apriori * exposure`` times the share still to develop, ``1 - 1 / cdf``, as chainladder-python's ``BornhuetterFerguson``. The exposure is a measure column of the same triangle (premium, say): each origin's latest observed cumulative value in the segment fitted. Parameters ---------- apriori : float, default 1.0 Expected loss ratio: the expected ultimate per unit of exposure; positive. average : {"volume", "simple", "regression"}, default "volume" How link ratios are averaged, as in ``ChainLadder``. sigma_interpolation : {"log-linear", "mack"}, default "log-linear" tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional As ``ChainLadder``; no tail by default. Examples -------- >>> from prospicio.reserving import BornhuetterFerguson, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021], [12, 24, 12], ... {"paid": [100.0, 150.0, 200.0], "premium": [250.0, 250.0, 400.0]}, ... ) >>> fit = BornhuetterFerguson(apriori=0.5).fit(tri, "paid", "premium") >>> [round(u, 2) for u in fit.ultimate] [150.0, 266.67] Benktander(apriori=1.0, n_iters=1, average='volume', sigma_interpolation='log-linear', tail=None) The Benktander (iterated Bornhuetter–Ferguson) method: starting from ``U(0) = apriori * exposure``, ``U(k) = latest + (1 - 1 / cdf) * U(k-1)`` for ``n_iters`` steps, as chainladder-python's ``Benktander``. ``n_iters=0`` is the expected loss method, 1 is Bornhuetter–Ferguson, and many iterations approach the chain ladder. The steps are summed in closed form, so a large ``n_iters`` is cheap; where an origin's ``cdf`` is below 1/2 they diverge instead. Parameters ---------- apriori : float, default 1.0 Expected loss ratio of the starting ultimate; positive. n_iters : int, default 1 Number of Bornhuetter–Ferguson steps. average : {"volume", "simple", "regression"}, default "volume" How link ratios are averaged, as in ``ChainLadder``. sigma_interpolation : {"log-linear", "mack"}, default "log-linear" tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional As ``ChainLadder``; no tail by default. Examples -------- >>> from prospicio.reserving import Benktander, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021], [12, 24, 12], ... {"paid": [100.0, 150.0, 200.0], "premium": [250.0, 250.0, 400.0]}, ... ) >>> fit = Benktander(apriori=0.5, n_iters=2).fit(tri, "paid", "premium") >>> [round(u, 2) for u in fit.ultimate] [150.0, 288.89] CapeCod(trend=0.0, decay=1.0, average='volume', sigma_interpolation='log-linear', tail=None) The Cape Cod (Stanard–Bühlmann) method: Bornhuetter–Ferguson with each origin's apriori estimated from the triangle, as chainladder-python's ``CapeCod``. Origin ``j``'s used-up exposure is ``exposure[j] / cdf[j]`` and its latest value is trended to the triangle's valuation by ``(1 + trend) ** (months / 12)``, the months running from the end of the origin period. Origin ``i``'s trended apriori is the sum of the trended latest values weighted by ``decay ** abs(i - j)`` over the same weighted sum of used-up exposures; dividing by its own trend factor gives the apriori of its Bornhuetter–Ferguson ultimate. Parameters ---------- trend : float, default 0.0 Annual trend of the loss ratio; above -1. decay : float, default 1.0 Weight of an origin ``n`` periods away, ``decay ** n``; from 0 to 1. With 1 every origin shares one loss ratio. average : {"volume", "simple", "regression"}, default "volume" How link ratios are averaged, as in ``ChainLadder``. sigma_interpolation : {"log-linear", "mack"}, default "log-linear" tail : float, TailConstant, TailCurve, TailBondy or TailLogLinear, optional As ``ChainLadder``; no tail by default. Examples -------- >>> from prospicio.reserving import CapeCod, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021], [12, 24, 12], ... {"paid": [100.0, 150.0, 200.0], "premium": [250.0, 250.0, 400.0]}, ... ) >>> fit = CapeCod().fit(tri, "paid", "premium") >>> [round(a, 4) for a in fit.apriori], [round(u, 2) for u in fit.ultimate] ([0.6774, 0.6774], [150.0, 290.32]) ExpectedLossFit A fitted expected-loss method (``ExpectedLoss``, ``BornhuetterFerguson`` or ``Benktander``) of every segment of a triangle column. Per-origin lists (``origins``, ``latest``, ``exposure``, ``apriori``, ``ultimate``, ``reserve``) run over the origins of each segment in turn, like the rows of ``to_frame()``. ``ultimate`` and ``reserve`` are this method's; ``chain_ladder`` holds the chain ladder's. Per-age lists need a single-segment fit; for several segments use ``development_frame()`` or ``segment(...)``. Examples -------- >>> from prospicio.reserving import BornhuetterFerguson, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021] * 2, ... [12, 24, 12] * 2, ... {"paid": [100.0, 150.0, 200.0, 10.0, 20.0, 30.0], ... "premium": [250.0, 250.0, 400.0, 500.0, 500.0, 800.0]}, ... keys={"lob": ["Auto"] * 3 + ["Home"] * 3}, ... ) >>> fit = BornhuetterFerguson(apriori=0.5).fit(tri, "paid", "premium") >>> fit.exposure, fit.segment(lob="Home").ultimate ([250.0, 400.0, 500.0, 800.0], [20.0, 230.0]) CapeCodFit A fitted Cape Cod of every segment of a triangle column: the fields of ``ExpectedLossFit``, with ``apriori`` the detrended loss ratio applied to each origin (chainladder-python's ``detrended_apriori_``), plus ``trended_apriori`` before detrending (its ``apriori_``). Per-origin lists run over the origins of each segment in turn, like the rows of ``to_frame()``. Examples -------- >>> from prospicio.reserving import CapeCod, Triangle >>> tri = Triangle.from_long( ... [2020, 2020, 2021], [12, 24, 12], ... {"paid": [100.0, 150.0, 200.0], "premium": [250.0, 250.0, 400.0]}, ... ) >>> fit = CapeCod(trend=0.1).fit(tri, "paid", "premium") >>> round(fit.trended_apriori[0] / fit.apriori[0], 10) 1.1 OdpBootstrap(n_sims=10000, seed=0, process='gamma') Over-dispersed Poisson bootstrap of the chain ladder (England and Verrall 2002), as R ChainLadder's ``BootChainLadder``: adjusted Pearson residuals of the volume-weighted chain ladder are resampled into pseudo triangles, each is re-projected, and process error is added to every future incremental value. Simulation ``i`` uses random stream ``i`` of ``seed`` for every segment in turn, so results do not depend on the number of threads. Parameters ---------- n_sims : int, default 10000 Number of simulations; positive. seed : int, default 0 Seed of the simulation streams, from 0 to ``2**64 - 1``. R accepts seeds below ``2**53``; a seed in both ranges gives the same draws. process : {"gamma", "none"}, default "gamma" Process error on each simulated future incremental value: Gamma with the expected value as mean and variance ``scale * |mean|`` (R's ``process.distr = "gamma"``), or none for parameter error only. Raises ------ ValueError If ``n_sims`` is zero or ``process`` is unknown. OverflowError If ``n_sims`` or ``seed`` is negative or too large. Examples -------- >>> from prospicio.reserving import OdpBootstrap, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> fit = OdpBootstrap(n_sims=2000, seed=42).fit(tri, "values") >>> fit.reserves.components() [('2020',), ('2021',), ('2022',), ('2023',)] >>> fit.reserves.mean() > 0 True OdpBootstrapFit A fitted ODP bootstrap of every segment. ``reserves`` is one joint distribution with the triangle's keys and ``"origin"`` as dimensions, so ``reserves.aggregate(["lob"])`` keeps the dependence between segments. Per-origin lists run over the origins of each segment in turn, like the rows of ``to_frame()`` and the components of ``reserves``. ``fitted``, ``residuals`` and ``scale`` need a single-segment fit; for several segments use ``segment(...)`` or ``totals_frame()``. ``fitted`` and ``residuals`` are nested lists indexed ``[origin][development]``, like one segment of ``Triangle.values``, with ``nan`` where the triangle is not observed. MackBootstrap(n_sims=10000, seed=0, process='gamma', average='volume', sigma_interpolation='log-linear', centre_residuals=True) Mack's bootstrap for the lifetime and one-year views (England, Verrall and Wüthrich 2019, Appendix 1): the scaled bias-adjusted residuals of the link ratios are resampled into pseudo factors, and each future cumulative value, to the last age (``fit``) or over the coming year (``one_year``), is drawn from the one before ``C`` (the observed latest value for the first) with mean ``f* C`` and Mack's variance ``sigma**2 * abs(C)**(2 - alpha)``. The lifetime view's standard deviation approximates Mack's analytic standard error. Beside ``OdpBootstrap`` (variance ``scale`` times the mean increment), it gives the one-year view under Mack's process: with the volume-weighted chain ladder and no tail, its standard deviation is ``MackFit.claims_development_result()``'s (Merz and Wüthrich) within Monte Carlo error. The residuals are centred by default (``centre_residuals``), so the mean CDR is Merz and Wüthrich's zero and the lifetime mean reserve the chain ladder's, and EVW's Table 4 expected reserves agree. Uncentred, as EVW's Appendix 1 is written, the pool's non-zero mean biases the pseudo factors: the mean CDR is about -0.2 (RAA), -0.04 (GenIns) and +0.18 (ABC) times its standard deviation, the lifetime mean reserve about +17%, +0.7% and -0.8% off the chain ladder's, and the one-year standard deviation up to 1.3% wide on RAA. Simulation ``i`` uses random stream ``i`` of ``seed`` for every segment in turn. Parameters ---------- n_sims : int, default 10000 Number of simulations; positive. seed : int, default 0 Seed of the simulation streams, from 0 to ``2**64 - 1``. process : {"gamma", "lognormal", "residuals", "normal", "none"}, default "gamma" Process error on each next cumulative value: Gamma or lognormal (negated for a negative mean) or normal, with Mack's mean and variance; the mean plus a resampled residual times the standard deviation, which carries the residuals' mean and variance; or none for parameter error only. average : {"volume", "simple", "regression"}, default "volume" How Mack's model averages the link ratios (its ``alpha``). sigma_interpolation : {"log-linear", "mack"}, default "log-linear" How a sigma behind a single link ratio is filled in. centre_residuals : bool, default True Subtract the residuals' mean before resampling them, so that the pseudo factors are unbiased, the mean CDR is about zero and the lifetime mean reserve is the chain ladder's. ``False`` resamples them uncentred, as EVW's Appendix 1 is written. Raises ------ ValueError If ``n_sims`` is zero or a setting is unknown. OverflowError If ``n_sims`` or ``seed`` is negative or too large. Examples -------- >>> from prospicio.reserving import ChainLadder, Mack, MackBootstrap, Triangle >>> tri = Triangle.from_long( ... [2020] * 4 + [2021] * 3 + [2022] * 2 + [2023], ... [12, 24, 36, 48, 12, 24, 36, 12, 24, 12], ... [100.0, 150.0, 165.0, 170.0, 110.0, 170.0, 180.0, 120.0, 175.0, 130.0], ... ) >>> fit = MackBootstrap(n_sims=2000, seed=42).one_year(tri, "values", ChainLadder()) >>> fit.model 'mack' >>> fit.cdr.variance() ** 0.5 < Mack().fit(tri, "values").total_standard_error True MackBootstrapFit A fitted bootstrap of Mack's model, the lifetime view, of every segment (``MackBootstrap.fit``). ``reserves`` is one joint distribution with the triangle's keys and ``"origin"`` as dimensions, so ``reserves.aggregate(["lob"])`` keeps the dependence between segments. Per-origin lists run over the origins of each segment in turn, like the rows of ``to_frame()`` and the components of ``reserves``. ``residuals`` needs a single-segment fit; for several segments use ``segment(...)``. OneYearFit The simulated one-year view of every segment, from ``OdpBootstrap.one_year`` or ``MackBootstrap.one_year`` (``model`` says which). ``cdr`` is one joint distribution of the claims development result with the triangle's keys and ``"origin"`` as dimensions, so ``cdr.aggregate(["lob"])`` keeps the dependence between segments, and ``cdr.quantile(0.005)`` is minus the one-year value at risk at 99.5%. Per-origin lists run over the origins of each segment in turn, like the rows of ``to_frame()`` and the components of ``cdr``. ``fitted``, ``residuals`` and ``scale`` need a single-segment fit; for several segments use ``segment(...)`` or ``totals_frame()``. ``fitted`` and ``scale`` are the ODP bootstrap's (``OdpBootstrapFit``'s), ``mack`` is Mack's bootstrap's model (a ``MackFit`` of every segment, as ``chain_ladder``), and ``residuals`` are either's. ClarkLdf(curve='loglogistic', max_age=None) Clark's LDF method (Clark 2003), as R ChainLadder's ``ClarkLDF``: each origin's expected ultimate and a growth curve are fitted to the incremental losses by over-dispersed Poisson maximum likelihood, with ages measured from the average date of loss (the middle of the origin period, R's ``adol = TRUE``). The ultimate is the latest value developed by the fitted curve to ``max_age``. Process risk is the scale times the fitted reserve, and parameter risk the delta method on the parameters' covariance, the scale times the inverse Fisher information. Parameters ---------- curve : {"loglogistic", "weibull"}, default "loglogistic" The growth curve ``G``: ``x**omega / (x**omega + theta**omega)`` or ``1 - exp(-(x / theta)**omega)``. max_age : float, optional Age in months at which development stops; at least the triangle's last age. ``None`` develops to infinity. Raises ------ ValueError If ``curve`` is unknown. Examples -------- >>> from prospicio.reserving import ClarkLdf, Triangle >>> rows = [[110.0, 290.0, 370.0, 420.0, 440.0], [95.0, 300.0, 390.0, 425.0], ... [130.0, 320.0, 410.0], [105.0, 305.0], [120.0]] >>> tri = Triangle.from_long( ... [2020 + i for i, row in enumerate(rows) for _ in row], ... [12 * (d + 1) for row in rows for d in range(len(row))], ... [v for row in rows for v in row], ... ) >>> fit = ClarkLdf(curve="weibull", max_age=120).fit(tri, "values") >>> fit.omega > 0 and fit.total_standard_error > fit.total_process_risk True >>> round(fit.ultimate[2] * fit.growth(36) / fit.growth(120), 6) 410.0 ClarkCapeCod(curve='loglogistic', max_age=None) Clark's Cape Cod method (Clark 2003), as R ChainLadder's ``ClarkCapeCod``: one expected loss ratio times each origin's exposure and a growth curve are fitted to the incremental losses by over-dispersed Poisson maximum likelihood, with ages measured from the average date of loss. The reserve is the fitted ``elr * exposure * (G(max_age) - G(age))``; process and parameter risk are as in ``ClarkLdf``. Parameters ---------- curve : {"loglogistic", "weibull"}, default "loglogistic" The growth curve, as in ``ClarkLdf``. max_age : float, optional Age in months at which development stops; at least the triangle's last age. ``None`` develops to infinity. Raises ------ ValueError If ``curve`` is unknown. Examples -------- >>> from prospicio.reserving import ClarkCapeCod, Triangle >>> rows = [[110.0, 290.0, 370.0, 420.0, 440.0], [95.0, 300.0, 390.0, 425.0], ... [130.0, 320.0, 410.0], [105.0, 305.0], [120.0]] >>> tri = Triangle.from_long( ... [2020 + i for i, row in enumerate(rows) for _ in row], ... [12 * (d + 1) for row in rows for d in range(len(row))], ... {"paid": [v for row in rows for v in row], ... "premium": [800.0 for row in rows for _ in row]}, ... ) >>> fit = ClarkCapeCod().fit(tri, "paid", "premium") >>> 0 < fit.elr < 1 and fit.expected_ultimate == [fit.elr * 800.0] * 5 True ClarkFit A fitted Clark LDF or Cape Cod model of every segment of a triangle column. Per-origin lists (``origins``, ``latest``, ``expected_ultimate``, ``ultimate``, ``reserve`` and the standard errors) run over the origins of each segment in turn, like the rows of ``to_frame()``. The fitted parameters and the standard errors of the total need a single-segment fit; for several segments use ``totals_frame()`` or ``segment(...)``. ``total_ultimate`` and ``total_reserve`` sum over every segment. Examples -------- >>> from prospicio.reserving import ClarkLdf, Triangle >>> rows = [[110.0, 290.0, 370.0, 420.0, 440.0], [95.0, 300.0, 390.0, 425.0], ... [130.0, 320.0, 410.0], [105.0, 305.0], [120.0]] >>> tri = Triangle.from_long( ... [2020 + i for i, row in enumerate(rows) for _ in row], ... [12 * (d + 1) for row in rows for d in range(len(row))], ... [v for row in rows for v in row], ... ) >>> fit = ClarkLdf().fit(tri, "values") >>> len(fit.covariance), fit.elr, fit.growth(float("inf")) (7, None, 1.0) ## Risk measures Distortion risk measures, and capital allocation to components. Distortion A distortion risk measure: ``rho(X) = integral of g(S(x)) dx`` for a concave distortion ``g`` of the survival function. Make one with ``Distortion.tvar``, ``Distortion.wang``, ``Distortion.proportional_hazard``, ``Distortion.dual_power`` or ``Distortion.exponential``. Every one is coherent, and each has a parameter value that gives the mean (or a limit that does). Examples -------- >>> from prospicio.distributions import Sampled >>> from prospicio.risk import Distortion >>> x = Sampled([1.0, 2.0, 3.0, 4.0]) >>> Distortion.tvar(0.5).measure(x) 3.5 >>> Distortion.tvar(0.5).weights(4) [0.0, 0.0, 0.5, 0.5] allocate(pd, distortion) Allocates a distortion risk measure of the total to the components. Euler allocation by co-measure: simulations are ranked by their total and each component gets the distortion-weighted sum of its own draws. The contributions sum to ``distortion.measure(pd)``; for ``Distortion.tvar(p)`` they are the CoTVaRs. Components must add up to the portfolio being allocated. Parameters ---------- pd : PredictiveDistribution distortion : Distortion Returns ------- list of float One contribution per component, in ``pd.components()`` order. Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import Distortion, allocate >>> pd = PredictiveDistribution(["lob"], [("motor",), ("property",)], ... [[1.0, 2.0], [4.0, 1.0], [2.0, 5.0], [3.0, 6.0]]) >>> allocate(pd, Distortion.tvar(0.5)) [2.5, 5.5] capital(pd, distortion, method='euler') Allocates the distortion risk measure of a portfolio's total to its components. Methods: - ``"euler"``: co-measure (for TVaR, the CoTVaRs); consistent with marginal changes to the portfolio. - ``"covariance"``: ``rho(S) Cov(X_j, S) / Var(S)``. - ``"proportional"``: stand-alone measures scaled to ``rho(S)``. - ``"marginal"``: ``rho(S) - rho(S - X_j)`` (Merton-Perold); does not add up to ``rho(S)``. - ``"shapley"``: Shapley value of ``v(T) = rho(sum of T)``; at most 12 components. Parameters ---------- pd : PredictiveDistribution Components that add up to the portfolio. distortion : Distortion method : str, default "euler" Returns ------- Allocation Raises ------ ValueError For an unknown method, a constant total (``"covariance"``), stand-alone measures summing to 0 (``"proportional"``) or more than 12 components (``"shapley"``). Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import Distortion, capital >>> pd = PredictiveDistribution(["lob"], [("motor",), ("property",)], ... [[1.0, 2.0], [4.0, 1.0], [2.0, 5.0], [3.0, 6.0]]) >>> a = capital(pd, Distortion.tvar(0.5)) >>> a.total, a.standalone, a.allocated (8.0, [3.5, 5.5], [2.5, 5.5]) >>> a.diversification_benefit() 1.0 entropic(dist, theta) Entropic risk measure ``(1 / theta) log E[exp(theta X)]``: the certainty equivalent of a loss under exponential utility. It rises from the mean (``theta -> 0``) to the largest draw (``theta -> inf``); for a normal loss it is ``mu + theta sigma**2 / 2``. Parameters ---------- dist : Sampled, PredictiveDistribution or list of float A predictive distribution is measured on its total. theta : float Risk aversion, positive. Returns ------- float Examples -------- >>> import math >>> from prospicio.risk import entropic >>> round(entropic([0.0, 1.0], math.log(2.0)), 12) == round(math.log2(1.5), 12) True esscher(dist, h) Esscher premium ``E[X exp(h X)] / E[exp(h X)]``: the mean after tilting probability towards large losses. The mean at ``h = 0``; ``mu + h sigma**2`` for a normal loss. Parameters ---------- dist : Sampled, PredictiveDistribution or list of float A predictive distribution is measured on its total. h : float Returns ------- float Examples -------- >>> import math >>> from prospicio.risk import esscher >>> round(esscher([0.0, 1.0], math.log(3.0)), 12) 0.75 marginal_expected_shortfall(pd, p) Marginal expected shortfall of each component at level ``p``: its mean over the simulations where the total is in its worst ``1 - p``. The same as ``allocate(pd, Distortion.tvar(p))``; it sums to the total's TVaR. Parameters ---------- pd : PredictiveDistribution p : float Returns ------- list of float One per component, in ``pd.components()`` order. Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import marginal_expected_shortfall >>> pd = PredictiveDistribution(["lob"], [("motor",), ("property",)], ... [[1.0, 2.0], [4.0, 1.0], [2.0, 5.0], [3.0, 6.0]]) >>> marginal_expected_shortfall(pd, 0.5) [2.5, 5.5] covar(pd, key, p, q) CoVaR of a component: the total's VaR at level ``q`` over the simulations where the component is at or above its own VaR at ``p``. Compare it with the total's unconditional VaR at ``q`` to see how much one segment's bad years drag the portfolio. Parameters ---------- pd : PredictiveDistribution key : tuple The component's key. p : float The component's distress level. q : float The level of the total's VaR. Returns ------- float Raises ------ ValueError If there is no component ``key``. Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import covar >>> pd = PredictiveDistribution(["lob"], [("a",), ("b",)], ... [[1.0, 0.0], [2.0, 1.0], [3.0, 5.0], [4.0, 1.0]]) >>> covar(pd, ("a",), 0.75, 0.5) 5.0 esscher_allocation(pd, h) Esscher allocation: each component's mean under the Esscher transform of the total, ``E[X_j exp(h S)] / E[exp(h S)]``. The contributions sum to ``esscher(pd, h)``; at ``h = 0`` they are the means. Parameters ---------- pd : PredictiveDistribution h : float Returns ------- list of float One per component, in ``pd.components()`` order. Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import esscher_allocation >>> pd = PredictiveDistribution(["lob"], [("motor",), ("property",)], ... [[1.0, 2.0], [4.0, 1.0], [2.0, 5.0], [3.0, 6.0]]) >>> esscher_allocation(pd, 0.0) [2.5, 3.5] Allocation Capital allocation of a distortion risk measure, from ``capital``. ## Dependence Copulas, and reordering existing draws to a target correlation. GaussianCopula(correlation) The Gaussian copula with correlation matrix ``correlation``. Parameters ---------- correlation : list of list of float Symmetric, unit diagonal, positive definite. Raises ------ ValueError If the matrix is not a valid correlation matrix. Examples -------- >>> from prospicio.risk import GaussianCopula >>> c = GaussianCopula([[1.0, 0.5], [0.5, 1.0]]) >>> u = c.sample(3, seed=1) >>> len(u), all(0.0 < x < 1.0 for row in u for x in row) (3, True) StudentTCopula(correlation, nu) The Student t copula with correlation matrix ``correlation`` and ``nu`` degrees of freedom: Gaussian-like correlation with joint extremes. Parameters ---------- correlation : list of list of float Symmetric, unit diagonal, positive definite. nu : float Degrees of freedom, positive. Raises ------ ValueError If the matrix or ``nu`` is invalid. Examples -------- >>> from prospicio.risk import StudentTCopula >>> StudentTCopula([[1.0, 0.5], [0.5, 1.0]], 4.0).dim 2 ArchimedeanCopula(family, theta, dim) An exchangeable Archimedean copula: Clayton, Gumbel, Frank or Joe. Parameters ---------- family : {"clayton", "gumbel", "frank", "joe"} theta : float Positive for Clayton and Frank; at least 1 for Gumbel and Joe. dim : int Raises ------ ValueError If the family is unknown or ``theta`` is out of range. Examples -------- >>> from prospicio.risk import ArchimedeanCopula >>> c = ArchimedeanCopula("clayton", 2.0, 3) # Kendall's tau 0.5 >>> c.dim, c.family (3, 'clayton') simulate(copula, marginals, n_sims, seed, keys=None, dims=None) Simulates marginals joined by a copula. In simulation ``i``, draws uniforms from ``copula`` with stream ``i`` of ``seed`` and applies each marginal's quantile function. Parameters ---------- copula : GaussianCopula, StudentTCopula or ArchimedeanCopula marginals : list of Lognormal, Grid or Pareto-family severities One per copula dimension. n_sims : int seed : int keys : list of tuple, optional One component key per marginal; defaults to ``(0,), (1,), ...``. dims : list of str, default ["component"] Returns ------- PredictiveDistribution Examples -------- >>> from prospicio.distributions import Lognormal >>> from prospicio.risk import GaussianCopula, simulate >>> c = GaussianCopula([[1.0, 0.4], [0.4, 1.0]]) >>> pd = simulate(c, [Lognormal.from_mean_cv(100.0, 0.2), Lognormal.from_mean_cv(50.0, 1.0)], ... 10_000, 42, keys=[("motor",), ("property",)], dims=["lob"]) >>> abs(pd.mean() - 150.0) < 3.0 True iman_conover(pd, correlation, seed) Reorders each component's draws to a target correlation (Iman-Conover). Every component keeps exactly its own draws; only their pairing across simulations changes. The correlation of the result's normal scores is close to ``correlation``, and Spearman's rho close to ``(6 / pi) asin(correlation / 2)``. Parameters ---------- pd : PredictiveDistribution correlation : list of list of float One row and column per component. seed : int Returns ------- PredictiveDistribution Examples -------- >>> from prospicio.distributions import PredictiveDistribution >>> from prospicio.risk import iman_conover >>> rows = [[float(i), float((i * 7919) % 1000)] for i in range(1000)] >>> pd = PredictiveDistribution(["lob"], [(0,), (1,)], rows) >>> joined = iman_conover(pd, [[1.0, 0.7], [0.7, 1.0]], seed=3) >>> sorted(joined.marginal((1,)).draws) == sorted(pd.marginal((1,)).draws) True ## Extreme value tails Generalized Pareto tails over a threshold, beyond the draws. Gpd(xi, beta) The generalized Pareto distribution, as SciPy's ``genpareto(c=xi, scale=beta)``. ``P(X > x) = (1 + xi x / beta)**(-1 / xi)`` for ``x >= 0``. Parameters ---------- xi : float Shape; moments of order ``1 / xi`` and above are infinite. beta : float Scale, positive. Examples -------- >>> from prospicio.risk import Gpd >>> g = Gpd(0.5, 2.0) >>> g.mean() 4.0 >>> fit = Gpd.fit([g.quantile((i - 0.5) / 1000) for i in range(1, 1001)]) >>> round(fit.xi, 2), round(fit.beta, 2) (0.5, 2.0) PotTail A peaks-over-threshold tail: draws above a threshold modelled by a fitted generalized Pareto distribution, for VaR and TVaR beyond the draws. Make one with ``PotTail.fit(draws, level)``, which takes the threshold at the empirical ``level`` quantile. Examples -------- >>> from prospicio.distributions import Lognormal, Sampled >>> from prospicio.risk import PotTail >>> d = Lognormal(0.0, 1.0) >>> s = Sampled([d.quantile((i - 0.5) / 100_000) for i in range(1, 100_001)]) >>> tail = PotTail.fit(s, 0.95) >>> abs(tail.var(0.999) / d.quantile(0.999) - 1) < 0.02 True mean_excess(draws, thresholds) The empirical mean-excess function ``e(u) = E[X - u | X > u]`` at each threshold, linear above a threshold where a GPD fits. Parameters ---------- draws : list of float thresholds : list of float Returns ------- list of (float, float, int) ``(u, e(u), number of draws above u)``; ``e(u)`` is NaN when none are. Examples -------- >>> from prospicio.risk import mean_excess >>> mean_excess([1.0, 2.0, 3.0, 4.0], [2.0]) [(2.0, 1.5, 2)] hill(draws, ks) Hill estimates of the tail index ``xi`` (``1 / alpha``) from the ``k`` largest draws, for each ``k``. Parameters ---------- draws : list of float ks : list of int Returns ------- list of float Raises ------ ValueError If a ``k`` is 0 or not below the number of draws, or the ``k + 1`` largest draws are not all positive.