r"""Stochastic reserving: the over-dispersed-Poisson bootstrap.
The chain ladder gives a point reserve and Mack (1993) gives its analytic
standard error (both in :mod:`actuarialpy.reserving`). This module adds the
full *predictive distribution* of the unpaid claims by the semiparametric
bootstrap of England & Verrall (1999, 2002): resample the fitted Pearson
residuals to capture estimation error, then draw each future cell from its
process distribution to capture process error.
Model
-----
Incremental claims :math:`C_{ij}` (origin :math:`i`, development :math:`j`)
follow an over-dispersed Poisson: :math:`\mathrm{E}[C_{ij}] = m_{ij}`,
:math:`\mathrm{Var}[C_{ij}] = \phi\, m_{ij}`, with a multiplicative mean
:math:`m_{ij} = x_i y_j`. The maximum-likelihood :math:`\hat m_{ij}` equals the
volume-weighted chain-ladder fitted incrementals -- so this bootstrap is the
stochastic layer on top of the chain ladder, not a different model, and its
mean reproduces the chain-ladder reserve.
Pearson residuals on the observed cells, degrees-of-freedom adjusted,
.. math::
r_{ij} = \frac{C_{ij} - \hat m_{ij}}{\sqrt{\hat m_{ij}}}, \qquad
r^{\mathrm{adj}}_{ij} = \sqrt{\tfrac{n}{n-p}}\; r_{ij}, \qquad
\hat\phi = \frac{1}{n-p} \sum_{ij} r_{ij}^2,
with :math:`n = K(K+1)/2` observed cells and :math:`p = 2K-1` parameters.
Bootstrap
---------
For each replicate: (1) resample residuals and form a pseudo-triangle
:math:`C^*_{ij} = \hat m_{ij} + r^*_{ij}\sqrt{\hat m_{ij}}`; (2) refit the
chain ladder and project future means :math:`\hat m^*_{ij}` -- *estimation
error*; (3) draw each future cell :math:`\tilde C_{ij} \sim
\mathrm{Gamma}(\mathrm{shape}=\hat m^*_{ij}/\hat\phi,\ \mathrm{scale}=\hat\phi)`
-- *process error*; (4) the replicate reserve is :math:`\sum \tilde C_{ij}`.
The bootstrap mean carries the usual small upward bias of a product of ratio
estimators, so report :attr:`BootstrapODP.point_reserve` (the chain-ladder
estimate) as the central figure and use the bootstrap only for the spread and
the tail. Because the ODP variance assumption (:math:`\mathrm{Var} \propto`
mean) differs from Mack's (:math:`\mathrm{Var} \propto C_{ik}`), the bootstrap
prediction error need not match Mack's; the two are different variance models
of the same quantity, which is exactly why both are worth having.
"""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
import pandas as pd
from .utils.random import RNGLike, resolve_rng
def _as_cumulative_array(triangle: pd.DataFrame | np.ndarray) -> tuple[np.ndarray, list]:
"""Return (cumulative KxK float array, origin labels)."""
if isinstance(triangle, pd.DataFrame):
origins = list(triangle.index)
cum = triangle.to_numpy(dtype=float)
else:
cum = np.asarray(triangle, dtype=float)
origins = list(range(cum.shape[0]))
if cum.ndim != 2 or cum.shape[0] != cum.shape[1]:
raise ValueError(
f"triangle must be a square K x K cumulative triangle, got shape {cum.shape}"
)
return cum, origins
def _dev_factors(cum: np.ndarray, valid: np.ndarray) -> np.ndarray:
"""Volume-weighted age-to-age factors f_k (k = 0..K-2), batched over
leading dims. ``cum`` and ``valid`` share shape (..., K, K)."""
num = np.where(valid[..., :, 1:], cum[..., :, 1:], 0.0).sum(axis=-2)
den = np.where(valid[..., :, 1:], cum[..., :, :-1], 0.0).sum(axis=-2)
with np.errstate(divide="ignore", invalid="ignore"):
return np.where(den > 0, num / den, 1.0)
def _fitted_incrementals(cum: np.ndarray, obs: np.ndarray) -> np.ndarray:
"""ODP maximum-likelihood fitted incrementals for a cumulative triangle.
Equivalent to the chain-ladder fit: a single fitted ultimate per origin
backed down through the shared development factors, then differenced. This
is the mean surface the Pearson residuals and the bootstrap are built on.
"""
K = cum.shape[0]
latest_dev = np.array([int(np.max(np.where(obs[i])[0])) for i in range(K)])
latest_cum = cum[np.arange(K), latest_dev]
f = _dev_factors(cum, obs)
F = np.append(np.cumprod(f[::-1])[::-1], 1.0)
fitted_ult = latest_cum * F[latest_dev]
fit_cum = np.empty((K, K))
fit_cum[:, K - 1] = fitted_ult
for k in range(K - 2, -1, -1):
fit_cum[:, k] = fit_cum[:, k + 1] / f[k]
return np.column_stack([fit_cum[:, 0], np.diff(fit_cum, axis=1)])
[docs]
@dataclass
class ReserveDistribution:
"""A sample from the predictive distribution of the unpaid claims.
Wraps the simulated total reserves (and, optionally, the per-origin
reserves) with the chain-ladder point estimate and convenience summaries.
Risk measures use the same empirical estimators as the rest of the
ecosystem (inverted-CDF VaR, Acerbi-Tasche TVaR).
"""
reserves: np.ndarray # (n_sims,) total unpaid per replicate
point_reserve: float # chain-ladder point estimate
by_origin: np.ndarray | None = None # (n_sims, K) per-origin reserves
origins: list | None = None
@property
def n_sims(self) -> int:
return int(self.reserves.shape[0])
def mean(self) -> float:
return float(self.reserves.mean())
[docs]
def prediction_error(self) -> float:
"""Standard deviation of the predictive distribution (the reserve SE)."""
return float(self.reserves.std(ddof=1))
def cv(self) -> float:
return self.prediction_error() / self.mean()
[docs]
def quantile(self, q) -> np.ndarray | float:
"""Empirical VaR (lower-quantile / inverted-CDF convention)."""
return np.quantile(self.reserves, q, method="inverted_cdf")
var = quantile
[docs]
def tvar(self, q: float) -> float:
"""Empirical tail value-at-risk (Acerbi-Tasche)."""
x = np.sort(self.reserves)
n = x.size
k = int(np.ceil(n * q))
if k >= n:
return float(x[-1])
return float((x[k:].sum() + x[k - 1] * (k - n * q)) / (n * (1 - q)))
def summary(self, quantiles=(0.5, 0.75, 0.95, 0.99, 0.995)) -> dict:
out = {
"point_reserve": self.point_reserve,
"mean": self.mean(),
"prediction_error": self.prediction_error(),
"cv": self.cv(),
}
for q in quantiles:
out[f"var_{q}"] = float(self.quantile(q))
return out
[docs]
def to_frame(self, quantiles=(0.5, 0.75, 0.95, 0.99, 0.995)) -> pd.DataFrame:
"""Per-origin and total predictive summary as a tidy table."""
if self.by_origin is None:
raise ValueError("no per-origin sample; build with by_origin=True")
rows = []
labels = self.origins if self.origins is not None else range(self.by_origin.shape[1])
cols = np.column_stack([self.by_origin, self.reserves])
for j, name in enumerate(list(labels) + ["Total"]):
col = cols[:, j]
row = {
"origin": name,
"reserve": float(col.mean()),
"se": float(col.std(ddof=1)),
}
for q in quantiles:
row[f"q{int(q * 100)}"] = float(np.quantile(col, q, method="inverted_cdf"))
rows.append(row)
return pd.DataFrame(rows).set_index("origin")
[docs]
@dataclass
class BootstrapODP:
"""England-Verrall over-dispersed-Poisson bootstrap fitted to a triangle.
Fit with :meth:`fit` on a cumulative development triangle, then call
:meth:`sample` for reserve draws (the ``risksim`` seam) or
:meth:`reserve_distribution` for a summarized predictive distribution.
"""
fitted_incr: np.ndarray # (K, K) ODP fitted incrementals
residuals: np.ndarray # 1-D pool of adjusted Pearson residuals
dispersion: float # Pearson phi
obs_mask: np.ndarray # (K, K) True on observed (upper) cells
latest_dev: np.ndarray # (K,) latest observed development index per origin
point_ibnr: np.ndarray # (K,) chain-ladder IBNR per origin
origins: list
# ---- construction -----------------------------------------------------
[docs]
@classmethod
def fit(cls, triangle: pd.DataFrame | np.ndarray) -> "BootstrapODP":
"""Fit the ODP model to a cumulative development triangle.
``triangle`` is a square ``K x K`` cumulative triangle -- a DataFrame
(origins on the index, developments on the columns) or a 2-D array --
with the unobserved lower-right as ``NaN``. This is the same input
:class:`reservingmodels.ChainLadder` takes; build one from a claims
listing with :func:`reservingmodels.make_completion_triangle` or from
an ``actuarialpy.Experience`` via
:func:`reservingmodels.integrations.actuarialpy.triangle_from_experience`.
"""
cum, origins = _as_cumulative_array(triangle)
K = cum.shape[0]
obs = ~np.isnan(cum)
latest_dev = np.array([int(np.max(np.where(obs[i])[0])) for i in range(K)])
latest_cum = cum[np.arange(K), latest_dev]
f = _dev_factors(cum, obs) # (K-1,)
F = np.append(np.cumprod(f[::-1])[::-1], 1.0) # F_k = prod_{l>=k} f_l
fitted_ult = latest_cum * F[latest_dev]
fitted_incr = _fitted_incrementals(cum, obs)
good = obs & (fitted_incr > 0)
incr_obs = np.column_stack([cum[:, 0], np.diff(np.nan_to_num(cum), axis=1)])
r = (incr_obs[good] - fitted_incr[good]) / np.sqrt(fitted_incr[good])
n_obs = int(obs.sum())
dof = n_obs - (2 * K - 1)
if dof <= 0:
raise ValueError(
f"triangle too small for the ODP bootstrap: {n_obs} cells, "
f"{2 * K - 1} parameters (need n > 2K-1)"
)
dispersion = float((r ** 2).sum() / dof)
r_adj = r * np.sqrt(n_obs / dof)
return cls(
fitted_incr=fitted_incr,
residuals=r_adj,
dispersion=dispersion,
obs_mask=obs,
latest_dev=latest_dev,
point_ibnr=fitted_ult - latest_cum,
origins=origins,
)
# ---- point estimate ---------------------------------------------------
@property
def point_reserve(self) -> float:
"""The chain-ladder reserve (sum of per-origin IBNR)."""
return float(self.point_ibnr.sum())
# ---- the risksim seam -------------------------------------------------
[docs]
def sample(self, size: int = 1, rng: RNGLike = None) -> np.ndarray:
"""Draw ``size`` predictive total-reserve outcomes -> shape ``(size,)``.
One draw is one simulated total unpaid amount, so a fitted
``BootstrapODP`` drops straight into ``risksim.PortfolioItem`` as an
aggregate-loss component. ``rng`` is ``None`` (global state), an
``int`` seed, or a ``numpy.random.Generator`` threaded through a
portfolio simulation.
"""
return self.sample_by_origin(size, rng).sum(axis=1)
[docs]
def sample_by_origin(self, size: int = 1, rng: RNGLike = None) -> np.ndarray:
"""Predictive reserve per origin -> shape ``(size, K)``."""
gen = resolve_rng(rng)
obs = self.obs_mask
K = obs.shape[0]
ld = self.latest_dev
phi = self.dispersion
# 1. resample residuals onto observed cells -> pseudo incrementals
m_obs = self.fitted_incr[None, :, :] * obs[None, :, :]
eps = gen.choice(self.residuals, size=(size, K, K), replace=True)
pseudo_incr = np.where(
obs[None], m_obs + eps * np.sqrt(np.clip(m_obs, 0, None)), 0.0
)
# 2. cumulate and refit factors on each pseudo-triangle
pseudo_cum = np.cumsum(pseudo_incr, axis=2)
fstar = _dev_factors(pseudo_cum, np.broadcast_to(obs, (size, K, K)))
# 3-4. project future means from the pseudo leading edge, add process error
reserves = np.zeros((size, K))
for i in range(K):
if ld[i] >= K - 1:
continue # origin fully developed
lead = pseudo_cum[:, i, ld[i]] # (size,)
fac = np.cumprod(fstar[:, ld[i]:], axis=1) # (size, K-1-ld)
fut_cum = lead[:, None] * fac
prev = np.column_stack([lead, fut_cum[:, :-1]])
mstar = fut_cum - prev # future incremental means
shape = np.clip(mstar, 1e-9, None) / phi
draws = gen.gamma(shape=shape, scale=phi)
draws = np.where(mstar > 0, draws, 0.0)
reserves[:, i] = draws.sum(axis=1)
return reserves
# ---- convenience ------------------------------------------------------
[docs]
def reserve_distribution(
self, size: int = 100_000, rng: RNGLike = None
) -> ReserveDistribution:
"""Materialize a summarized predictive distribution of the reserve."""
by = self.sample_by_origin(size, rng)
return ReserveDistribution(
reserves=by.sum(axis=1),
point_reserve=self.point_reserve,
by_origin=by,
origins=self.origins,
)