From Competitor Quotes to Risk Parameters: Implementing Market-Based Ratemaking in Python

import numpy as np
from scipy.optimize import isotonic_regression

# ---------------------------------------------------------------------------
# 1. Synthetic UK pet insurance data — 12 risk classes (4 breeds × 3 ages)
#    Market quotes (£/year) loosely calibrated to UK PCW rates, 2024.
#    Excess: £99. Coverage: vet fees up to £4,000.
# ---------------------------------------------------------------------------

BREEDS = ["Australian Shepherd", "Golden Retriever", "German Shepherd", "French Bulldog"]
AGES   = ["1yr", "4yr", "8yr"]

# Risk class labels (12 total)
risk_classes = [(b, a) for b in BREEDS for a in AGES]

# Observed competitor quotes (£/yr) — median of 5 UK insurers, synthetic.
# Ordering reflects: breed risk (Frenchie > GSD > Golden > Aussie),
# age risk (older = higher), consistent with BSAVA morbidity data.
market_quotes = np.array([
    # Aussie:  1yr    4yr    8yr
    180.0, 220.0, 310.0,
    # Golden:  1yr    4yr    8yr
    195.0, 245.0, 350.0,
    # GSD:     1yr    4yr    8yr
    215.0, 270.0, 385.0,
    # Frenchie: 1yr   4yr    8yr
    310.0, 390.0, 540.0,
])

# Coverage function parameters per risk class (deductible d, limit l, rate r)
# Simplified: r=1, d=99, l=4000 for all classes in this example.
DEDUCTIBLE = 99.0
LIMIT      = 4_000.0

# ---------------------------------------------------------------------------
# 2. Pure premium simulator — Poisson-LogNormal, vectorised
# ---------------------------------------------------------------------------

def simulate_pure_premiums(lam, mu, sigma=1.0, n_classes=12, R=2_000, rng=None):
    """
    For each risk class (differentiated by lam multipliers), simulate R
    Poisson-LogNormal aggregate losses and apply the coverage function.

    lam: base Poisson frequency (float)
    mu:  log-mean of severity (float)
    Returns: shape (n_classes,) array of expected net covered losses
    """
    if rng is None:
        rng = np.random.default_rng()

    # Risk class multipliers — frequency loadings by breed × age.
    # In practice you'd estimate these; here they encode the prior ordering.
    freq_mult = np.array([
        1.00, 1.15, 1.40,  # Aussie: low risk, rising with age
        1.10, 1.30, 1.60,  # Golden
        1.20, 1.45, 1.75,  # GSD
        1.70, 2.00, 2.40,  # Frenchie: highest risk
    ])

    pure_premiums = np.empty(n_classes)
    for i in range(n_classes):
        lam_i = lam * freq_mult[i]
        # Draw claim counts: shape (R,)
        n_claims = rng.poisson(lam_i, size=R)
        # Draw severities for all claims in batch
        total_claims = int(n_claims.sum())
        if total_claims == 0:
            pure_premiums[i] = 0.0
            continue
        severities = rng.lognormal(mean=mu, sigma=sigma, size=total_claims)
        # Apply coverage: min(max(sev - d, 0), l)
        covered = np.clip(severities - DEDUCTIBLE, 0, LIMIT)
        # Assign back to simulations
        agg = np.zeros(R)
        idx = 0
        for j, nc in enumerate(n_claims):
            agg[j] = covered[idx:idx + nc].sum()
            idx += nc
        pure_premiums[i] = agg.mean()

    return pure_premiums


# ---------------------------------------------------------------------------
# 3. Isotonic link — PAVA maps pure premiums to commercial premiums
# ---------------------------------------------------------------------------

def isotonic_link(pure_premiums, market_quotes):
    """
    Fit a monotone mapping f: pure_premium -> commercial_premium via PAVA.
    Returns fitted commercial premiums (same shape as pure_premiums).
    """
    # Sort by pure premium to define the monotonicity direction
    order = np.argsort(pure_premiums)
    sorted_quotes = market_quotes[order]
    # scipy.optimize.isotonic_regression minimises weighted sum of squares
    # subject to non-decreasing constraint
    result = isotonic_regression(sorted_quotes, increasing=True)
    # Map back to original order
    fitted = np.empty_like(pure_premiums)
    fitted[order] = result.x
    return fitted


# ---------------------------------------------------------------------------
# 4. Loss ratio penalty — the corridor enforces actuarial plausibility
# ---------------------------------------------------------------------------

def loss_ratio_penalty(pure_premiums, market_quotes, lr_low=0.50, lr_high=0.75):
    """
    Penalise particles where implied loss ratios (pure/commercial) fall
    outside the corridor [lr_low, lr_high].

    UK calibration: pet/travel [50%, 75%], home [45%, 65%], motor [55%, 75%].
    """
    lr = pure_premiums / market_quotes
    below = np.maximum(lr_low - lr, 0).sum()
    above = np.maximum(lr - lr_high, 0).sum()
    return below + above


# ---------------------------------------------------------------------------
# 5. ABC-SMC main loop
# ---------------------------------------------------------------------------

def abc_smc(
    market_quotes,
    n_particles=500,
    n_sim=1_000,
    n_generations=7,
    tolerance_init=200.0,
    tolerance_decay=0.65,
    lr_low=0.50,
    lr_high=0.75,
    seed=42,
):
    """
    Population Monte Carlo ABC for market-based ratemaking.

    Returns: particles (n_particles, 2) array of (lam, mu) samples
             representing the approximate posterior.
    """
    rng = np.random.default_rng(seed)

    # Prior: Uniform over actuarially plausible ranges
    # lam ~ Uniform(0.05, 1.5),  mu ~ Uniform(4.5, 8.0)
    def sample_prior(n, rng):
        lam = rng.uniform(0.05, 1.5, size=n)
        mu  = rng.uniform(4.5,  8.0, size=n)
        return np.column_stack([lam, mu])

    # Generation 0: sample from prior, evaluate, keep best n_particles
    tolerance = tolerance_init
    candidates = sample_prior(n_particles * 10, rng)
    distances  = np.full(len(candidates), np.inf)

    for idx, (lam, mu) in enumerate(candidates):
        pp = simulate_pure_premiums(lam, mu, R=n_sim, rng=rng)
        fitted = isotonic_link(pp, market_quotes)
        rmse = np.sqrt(np.mean((fitted - market_quotes) ** 2))
        penalty = loss_ratio_penalty(pp, market_quotes, lr_low, lr_high)
        distances[idx] = rmse + 50.0 * penalty

    accepted_idx = np.argsort(distances)[:n_particles]
    particles    = candidates[accepted_idx]
    weights      = np.ones(n_particles) / n_particles

    # Generations 1..G: resample, perturb, accept/reject
    for gen in range(1, n_generations):
        tolerance *= tolerance_decay
        cov = 2.0 * np.cov(particles.T)  # Gaussian perturbation kernel

        new_particles = np.empty_like(particles)
        accepted = 0

        while accepted < n_particles:
            # Resample from current particles
            i = rng.choice(n_particles, p=weights)
            candidate = particles[i] + rng.multivariate_normal([0, 0], cov)
            lam, mu = candidate

            # Prior bounds check
            if not (0.05 < lam < 1.5 and 4.5 < mu < 8.0):
                continue

            pp = simulate_pure_premiums(lam, mu, R=n_sim, rng=rng)
            fitted = isotonic_link(pp, market_quotes)
            rmse = np.sqrt(np.mean((fitted - market_quotes) ** 2))
            penalty = loss_ratio_penalty(pp, market_quotes, lr_low, lr_high)
            dist = rmse + 50.0 * penalty

            if dist < tolerance:
                new_particles[accepted] = candidate
                accepted += 1

        particles = new_particles
        weights   = np.ones(n_particles) / n_particles
        print(f"  Gen {gen}: tolerance={tolerance:.1f}, accepted {n_particles}")

    return particles


# ---------------------------------------------------------------------------
# 6. Run it
# ---------------------------------------------------------------------------

if __name__ == "__main__":
    print("Running ABC-SMC...")
    posterior = abc_smc(market_quotes, n_particles=500, n_sim=1_000, n_generations=6)

    lam_map = np.median(posterior[:, 0])
    mu_map  = np.median(posterior[:, 1])
    print(f"\nMAP estimates:  lambda={lam_map:.3f},  mu={mu_map:.3f}")
    print(f"Posterior sd:   lambda={posterior[:, 0].std():.3f}, "
          f"mu={posterior[:, 1].std():.3f}")

    # Validate: pure premiums at MAP, loss ratios
    pp_map = simulate_pure_premiums(lam_map, mu_map, R=5_000)
    lr     = pp_map / market_quotes
    print(f"\nLoss ratio range: {lr.min():.1%} – {lr.max():.1%}")
    print(f"Risk ordering check (French Bulldog > GSD): "
          f"{pp_map[9] > pp_map[6]}")  # 4yr Frenchie vs 4yr GSD

from insurance_credibility import BuhlmannStraub

# ABC MAP estimates — feed directly as collective mean
lam_map = 0.28   # posterior median
mu_map  = 6.21   # posterior median

# Implied collective mean pure premium for a standard risk class
# (e.g. 4yr Australian Shepherd, the reference class)
collective_mean_pp = np.exp(mu_map + 0.5) * lam_map  # approx. £247/yr

# Convert to loss rate (pp / reference premium)
reference_premium = 220.0  # observed market quote, 4yr Aussie
collective_mean_lr = collective_mean_pp / reference_premium  # ~0.62

# Posterior variance informs the between-class variance a
# Wide posterior → high a → more credibility weight needed before data dominates
posterior_var_lam = posterior[:, 0].var()
posterior_var_mu  = posterior[:, 1].var()

# Bühlmann-Straub: at launch, Z=0, price = market prior.
# As your own claims arrive, Z rises and your data progressively dominates.
bs = BuhlmannStraub(
    collective_mean=collective_mean_lr,
    # k = v/a where v is within-class variance (claims noise)
    # and a is between-class variance. Set conservatively at launch.
    k=8.0,   # implies ~8 policy-years before Z reaches 0.5
)

# At launch: Z=0.0, credibility-weighted estimate = collective_mean entirely
print(bs.credibility_factor(n=0))    # 0.0 — pure market prior
print(bs.credibility_factor(n=50))   # ~0.86 — mostly own data after 50 years
print(bs.credibility_factor(n=8))    # 0.50 — equal weight at k policy-years