risksim

Portfolio Monte Carlo simulation and risk measures: simulate aggregate outcomes across a portfolio of contracts and summarize the distribution with standard risk measures.

The package builds a portfolio from contract definitions, simulates aggregate outcomes, and reports risk measures (such as value-at-risk and tail value-at-risk) over the simulated distribution — the capital view that sits at the end of the experience → pricing → loss → tail → capital pipeline.

var and tvar follow the ecosystem-wide empirical estimators (inverted-CDF order statistic; Acerbi–Tasche), so a portfolio’s risk measures here match the same quantities computed in lossmodels or extremeloss byte for byte, and every simulation accepts the shared rng argument (None, seed, or Generator) for bit-reproducible runs — see Conventions.

See the API reference below for the full surface; each object’s docstring carries its own usage.

Portfolio simulation

A portfolio is a list of named items, each wrapping anything with a .sample method — a fitted lossmodels.CollectiveRiskModel drops straight in:

import risksim as rs

port = rs.Portfolio([
    rs.PortfolioItem("commercial", crm),
    rs.PortfolioItem("surety", other_model, weight=0.4),
])
res = port.simulate(100_000, rng=7)
res.gross_losses, res.component_losses, res.component_names

simulate returns a SimulationResult carrying the gross, ceded, and retained loss vectors, the per-component draws, and — when a contract is applied — the per-layer recoveries. The rng argument (None, seed, or Generator) resolves to one generator threaded through the components, so the same seed reproduces every array bit for bit.

Contracts

AggregateLayer(attachment, limit, share) follows the ecosystem coverage semantics: limit is the layer width, so the layer exhausts at attachment + limit (see Conventions). ContractProgram stacks layers, and apply_contract applies either to any loss vector, returning (ceded, retained):

treaty = rs.AggregateLayer(attachment=3_200_000, limit=1_500_000,
                           name="agg_stop_loss")
res = port.simulate(100_000, contract=treaty, rng=7)
res.ceded_losses, res.retained_losses, res.layer_losses

ceded, retained = rs.apply_contract([50.0, 150.0, 400.0],
                                    rs.AggregateLayer(attachment=100.0, limit=200.0))
# ceded    -> [  0.,  50., 200.]
# retained -> [ 50., 100., 200.]

Risk measures

rs.metrics.var(res.gross_losses, 0.99)
rs.metrics.tvar(res.retained_losses, 0.99)

These are the ecosystem-wide empirical estimators — see the intro above and Conventions.

Dependence between components

Independent components overstate diversification in exactly the tail metrics above. risksim.dependence.impose_rank_correlation fixes the default without touching any sampler: simulate each component as usual, then reorder (Iman-Conover) to a target rank correlation – marginals preserved exactly:

from risksim.dependence import impose_rank_correlation

matrix = np.column_stack([item.sample(n, rng) for item in items])
total = impose_rank_correlation(matrix, corr, rng).sum(axis=1)

Rank correlation is not tail dependence: normal scores leave joint extremes asymptotically independent at any rho. When the question is “do the components blow up together”, pass scores="t" with a small df – same rank correlation, genuinely clustered joint tails. Example 9 measures both effects on a two-line portfolio — the diversification benefit halving at ρ = 0.5, and the joint-exceedance probability separating the two score choices that every sum-based metric cannot tell apart.

Monte Carlo error

A simulated VaR without an error estimate is a random number with confidence. risksim.uncertainty answers “how much of this is signal” — each metric with the interval its sampling theory supports: normal theory for the mean, distribution-free order statistics for VaR (on the same ceil(n*q) rank convention as metrics.var, so points match exactly), percentile bootstrap for TVaR:

from risksim import uncertainty

uncertainty.summary_with_error(result.losses, quantiles=(0.95, 0.99), rng=7)
# {"mean": {...}, "var_95": {...}, "tvar_99": {estimate, se, ci_low, ci_high}}

uncertainty.quantile_ci(result.losses, q=0.99)   # order-statistic interval
uncertainty.bootstrap_ci(result.losses, lambda a: metrics.tvar(a, 0.99), rng=7)

quantile_ci reports se = nan deliberately: a quantile has no distribution-free standard error — the interval is the uncertainty statement. If the bands are too wide, the answer is more simulations, and now you can see it.

API reference

class AggregateLayer(attachment: float = 0.0, limit: float | None = None, share: float = 1.0, name: str | None = None)

Bases: object

Aggregate annual layer applied to simulated aggregate losses.

For aggregate annual loss S, ceded loss is:

C = share * min((S - attachment)+, limit)

with no cap if limit is None.

attachment_probability(losses: ndarray | list[float]) → float

Probability that a loss reaches the layer, i.e. P(loss > attachment).

This depends only on the attachment point, not on share or limit, so it stays correct under degenerate parameters (e.g. share=0).

apply_contract(losses: ndarray | list[float], contract: AggregateLayer | ContractProgram) → tuple[ndarray, ndarray]

Return (ceded, retained) arrays for a single aggregate layer or a multi-layer contract program.

class ContractProgram(layers: Sequence[AggregateLayer], name: str = 'contract_program')

Bases: object

Collection of aggregate layers applied to the same gross loss.

This first version assumes the layers are intended to work together without overlap. For a standard non-overlapping tower, total ceded loss is the row-wise sum of ceded loss by layer.

class Portfolio(items: Sequence[PortfolioItem], name: str = 'portfolio')

Bases: object

Portfolio of aggregate-loss components.

sample() and sample_components() draw the components independently – the natural default, and the right one when the components genuinely do not move together. Independence is not imposed on the workflow, only on this sampler: to introduce dependence, draw the component matrix here and reorder it with risksim.dependence.impose_rank_correlation(), which preserves each component’s marginal exactly while imposing a target rank correlation (see that function for the rank-vs-tail-dependence caveat):

from risksim.dependence import impose_rank_correlation

matrix = portfolio.sample_components(n, rng)
dependent = impose_rank_correlation(matrix, corr, rng)
total = dependent.sum(axis=1)

The analytic variance(), std(), and summary() methods assume independence and are unaffected by that post-processing; compute dependent risk measures from the reordered sample via risksim.metrics.

sample_components(size: int = 1, rng: Generator | int | None = None) → ndarray

Draw a (size, n_components) matrix, one column per component.

Components are drawn independently. This matrix is exactly the input risksim.dependence.impose_rank_correlation() expects, so imposing dependence is a one-line post-processing step on the return value (see the class docstring).

variance() → float

Analytic portfolio variance under the independence assumption.

This is the closed-form sum of component variances and does not reflect any dependence imposed downstream via risksim.dependence.impose_rank_correlation(); for a dependent portfolio, take the variance of the reordered sample.

class PortfolioItem(name: 'str', model: 'SupportsSample', weight: 'float' = 1.0)

Bases: object

class SimulationResult(gross_losses: ndarray, ceded_losses: ndarray | None = None, retained_losses: ndarray | None = None, component_losses: ndarray | None = None, component_names: Sequence[str] | None = None, layer_losses: ndarray | None = None, layer_names: Sequence[str] | None = None, contract_name: str | None = None)

Bases: object

Container for portfolio simulation outputs.

If retained_losses is present, the primary losses view is retained/net loss. Otherwise, the primary losses view is gross loss.

risksim.dependence

Dependence between simulated components, without touching the samplers.

Independence across portfolio components is the classically dangerous default: it overstates diversification in exactly the tail metrics this package exists to compute. This module adds dependence by reordering (Iman & Conover, 1982): simulate every component with whatever machinery already exists, then permute each column of the results so their ranks follow a target correlation. Marginals are preserved exactly – each column is a permutation of itself – so nothing about any component’s own distribution changes; only which scenarios coincide.

Two honest limits, stated loudly. First, the target is a rank correlation: with scores="normal" the induced Spearman correlation matches the requested matrix to within the usual (6/pi) * asin(rho/2) distortion (under 0.02 absolute) – but rank correlation is not tail dependence, and normal scores produce joint extremes that are asymptotically independent no matter how high rho is. If the risk question is “do the components blow up together”, use scores="t" with a small df: t scores put genuine mass on joint tail events at the same rank correlation. Second, this imposes the dependence you assert; it does not estimate dependence from data.

The portfolio recipe is two lines:

matrix = np.column_stack([item.sample(n, rng) for item in items])
total = impose_rank_correlation(matrix, corr, rng).sum(axis=1)

impose_rank_correlation(samples: ndarray, target_corr: ndarray, rng: Any = None, scores: str = 'normal', df: float = 5.0) → ndarray

Reorder simulated columns to a target rank correlation (Iman-Conover).

Parameters:

• samples (ndarray, shape (n_sims, n_components)) – Independently simulated component outcomes. Not modified.

• target_corr (ndarray, shape (k, k)) – Desired correlation matrix: symmetric, unit diagonal, positive semidefinite.

• rng (optional) – Seed or numpy.random.Generator for the latent scores.

• scores ({"normal", "t"}) – Latent score family. "normal" gives rank correlation with no tail dependence; "t" adds joint-tail clustering, stronger for smaller df.

• df (float) – Degrees of freedom for scores="t"; must exceed 2 (the scores need a finite variance for the correlation target to mean anything).

Returns:

Same shape as samples; each column an exact permutation of the corresponding input column.

Return type:

ndarray

risksim.metrics

mean(losses: ndarray | list[float]) → float

Empirical mean of a simulated loss vector.

variance(losses: ndarray | list[float], ddof: int = 0) → float

Empirical variance of a simulated loss vector.

ddof follows numpy’s convention: 0 (default) is the population estimator dividing by n; 1 divides by n - 1.

std(losses: ndarray | list[float], ddof: int = 0) → float

Empirical standard deviation of a simulated loss vector.

ddof follows numpy’s convention: 0 (default) is the population estimator dividing by n; 1 divides by n - 1.

var(losses: ndarray | list[float], q: float | ndarray) → float | ndarray

Empirical Value-at-Risk.

Uses the actuarial (lower-quantile) definition

VaR_q(X) = inf{ x : F(x) >= q },

whose empirical plug-in is the order statistic x_(ceil(n*q)). Equivalent to np.quantile(losses, q, method="inverted_cdf").

q may be a scalar (returns float) or array-like (returns np.ndarray of the same length).

tvar(losses: ndarray | list[float], q: float | ndarray) → float | ndarray

Empirical Tail Value-at-Risk (expected shortfall).

Uses the average-quantile (coherent) definition

TVaR_q(X) = (1 / (1 - q)) * integral_q^1 VaR_u(X) du,

whose empirical plug-in (Acerbi-Tasche) is, with sorted losses x_(1) <= ... <= x_(n) and k = ceil(n*q),

TVaR_q = [ sum_{i>k} x_(i) + x_(k) * (k - n*q) ] / (n * (1 - q)).

This is exact for the empirical distribution (correct with ties/atoms) and reduces to the mean of the largest n*(1-q) observations when n*q is an integer. Always satisfies tvar >= var.

q may be a scalar (returns float) or array-like (returns np.ndarray of the same length).

prob_exceeding(losses: ndarray | list[float], threshold: float) → float

Empirical exceedance probability P(X > threshold).

The inequality is strict, matching the survival-function convention used across the ecosystem. Also available under the ecosystem-standard spelling exceedance_probability.

exceedance_probability(losses: ndarray | list[float], threshold: float) → float

Empirical exceedance probability P(X > threshold).

The inequality is strict, matching the survival-function convention used across the ecosystem. Also available under the ecosystem-standard spelling exceedance_probability.

summary(losses: ndarray | list[float], quantiles: tuple[float, ...] = (0.95, 0.99)) → dict[str, Any]

One-call risk summary of a simulated loss vector.

Returns a dict with n_sims, mean, std, min, max, and a var_{p} / tvar_{p} pair per requested quantile, where {p} is the percentile label (0.99 -> "99", 0.995 -> "99.5"). VaR and TVaR use the ecosystem-wide empirical estimators — see var() and tvar().

risksim.uncertainty

Monte Carlo error quantification for simulation output.

A simulated VaR without an error estimate is a random number with confidence. Every function here takes the loss vector a simulation produced and answers “how much of this is signal”: normal-theory intervals for the mean, distribution-free order-statistic intervals for quantiles, and bootstrap intervals for anything else (TVaR in particular). Each metric gets the standard tool for that metric – mixing them in one summary is deliberate, because pretending one method fits all is how tail estimates end up with body-sized error bars.

All estimates use the same conventions as risksim.metrics (the lower-quantile order statistic for VaR), so the point values in summary_with_error() match risksim.metrics.summary() exactly.

bootstrap_ci(losses: ndarray | list[float], statistic: Callable[[ndarray], float], n_boot: int = 1000, confidence: float = 0.95, rng: Any = None) → dict[str, float]

Percentile-bootstrap confidence interval for any statistic.

Resamples the loss vector with replacement n_boot times and takes the empirical quantiles of the replicated statistic. The workhorse for statistics with no clean sampling theory – TVaR above all.

Parameters:

• losses (array-like) – Simulated losses.

• statistic (callable) – Maps a 1-d array to a float, e.g. lambda a: tvar(a, 0.99).

• n_boot (int) – Bootstrap replicates.

• confidence (float) – Interval level.

• rng (optional) – Seed or numpy.random.Generator for reproducibility.

Returns:

estimate (the statistic on the full sample – not the replicate mean), se (replicate standard deviation), ci_low, ci_high (percentile bounds).

Return type:

dict

mean_ci(losses: ndarray | list[float], confidence: float = 0.95) → dict[str, float]

Normal-theory confidence interval for the simulated mean.

The standard error is std(losses, ddof=1) / sqrt(n); with the sample sizes simulations run at, the normal interval is exact for all practical purposes.

Returns:

estimate, se, ci_low, ci_high.

Return type:

dict

quantile_ci(losses: ndarray | list[float], q: float, confidence: float = 0.95) → dict[str, float]

Distribution-free confidence interval for an empirical quantile.

The number of observations at or below the true q-quantile is Binomial(n, q); the interval takes the order statistics at ranks k -/+ z * sqrt(n q (1-q)) around the VaR rank k = ceil(n q) – no distributional assumption on the losses at all. Ranks clip to the sample: when n is too small for the requested tail, the bound honestly sits at the extreme order statistic rather than extrapolating.

se is nan by design: a quantile has no distribution-free standard error (its asymptotic variance involves the unknown density); the interval is the uncertainty statement.

Returns:

estimate, se (nan), ci_low, ci_high.

Return type:

dict

summary_with_error(losses: ndarray | list[float], quantiles: tuple[float, ...] = (0.95, 0.99), confidence: float = 0.95, n_boot: int = 1000, rng: Any = None) → dict[str, dict[str, float]]

risksim.metrics.summary(), with error bars on every metric.

Point estimates match metrics.summary exactly; each metric carries the interval its sampling theory supports: normal theory for the mean, order statistics for VaR (se is nan there – see quantile_ci()), bootstrap for TVaR.

Returns:

Keys like "mean", "var_95", "tvar_99"; each value has estimate, se, ci_low, ci_high.

Return type:

dict of str -> dict