claims = pd.read_csv("data/motor/claims.csv")
claims["incurred"] = claims["paid"] + claims["outstanding"]
liab = claims[claims["coverage"] == "LIAB"].merge(
    train[["vehicle_id", "vehicle_style", "age_band",
           "urban", "log_miles", "channel", "driver_age"]],
    on="vehicle_id", how="inner", validate="many_to_one")
liab["young"] = liab["driver_age"] < 25
len(liab)       # 4,041 claims from policy years 2021–2024

$$Y \sim \text{Gamma}(\mu, \phi) \qquad \mathrm{Var}(Y) = \phi\,\mu^2 \qquad \ln\mu = \mathbf{x}^\top\beta$$

gamma_log = sm.families.Gamma(sm.families.links.Log())
sev = smf.glm("incurred ~ C(vehicle_style, Treatment('Sedan'))"
              " + young + urban", data=liab,
              family=gamma_log).fit()

np.exp(sev.params)
sev.scale              # dispersion 1.832 → CV = √φ = 1.35

• Standard deviation ∝ mean: a $10,000 claim has 10× the spread of a $1,000 claim
• Specify the log link explicitly: the Gamma’s default is the inverse link

The claims were simulated lognormal, and the likelihood finds it; Gamma residuals bend away at both ends.

$$\text{Pure Premium} = \widehat{\text{Frequency}} \times \widehat{\text{Severity}}$$

Valid when claim count and claim size are independent given the features, the standard collective-risk assumption.

# predicted rate per car-year: offset of zero = one year
train["freq_hat"] = freq.predict(train, offset=np.zeros(len(train)))
train["sev_hat"] = sev.predict(train)
train["pp_hat"] = train["freq_hat"] * train["sev_hat"]

(train["pp_hat"] * train["exposure"]).sum()   # 51,773,072
train["inc_liab"].sum()                       # 51,794,591

$$\text{PP} = \underbrace{447.03}_{\text{base}} \times f_{\text{style}} \times f_{\text{age}} \times f_{\text{young}} \times f_{\text{urban}} \times \Big(\tfrac{\text{miles}}{12{,}000}\Big)^{0.488}$$

Three risks, priced

for d in (train, hold):
    d["pp_hat"] = (freq.predict(d, offset=np.zeros(len(d)))
                   * sev.predict(d))

expected = (hold["pp_hat"] * hold["exposure"]).sum()
hold["inc_liab"].sum() / expected       # A/E = 1.022

On policy year 2025, which the models never saw, actual losses were 2.2% above expected: 7,261,043 against 7,107,968.

$$S = \sum_{j=1}^{N} X_j,\quad N \sim \text{Poisson},\; X_j \sim \text{Gamma}$$

S is exactly zero when N = 0, and continuous otherwise. It belongs to the EDF with

$$\mathrm{Var}(S) = \phi\,\mu^p, \qquad 1 < p < 2$$

• p → 1: Poisson-like (counts)
• p → 2: Gamma-like (sizes)
• Model pure premium directly, with exposure as a weight

train["pp"] = train["inc_liab"] / train["exposure"]

tw = smf.glm("pp ~ C(vehicle_style, Treatment('Sedan'))"
             " + C(age_band, Treatment('40-49'))"
             " + urban + log_miles",
             data=train, var_weights=train["exposure"],
             family=sm.families.Tweedie(
                 var_power=1.6, link=sm.families.links.Log()),
            ).fit()

Choose p by profile likelihood: fit a grid, keep the best.

Holdout A/E: Tweedie 1.025, two-part 1.022.