$$y = \mathbf{X}\beta + \varepsilon \qquad \hat\beta = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top y$$

OLS picks the β that minimizes the sum of squared residuals $\sum_i (y_i - \mathbf{x}_i^\top\beta)^2$.

Gauss–Markov assumptions

1. Linearity — E[y] is linear in the parameters
2. Full rank — no feature is an exact combination of others
3. Exogeneity — E[ε | X] = 0
4. Homoscedasticity — Var(ε | X) = σ², the same for every row
5. No autocorrelation — errors are uncorrelated

• 180,472 encounters from 2024, one row each
• C(col, Treatment('level')) makes dummy variables and names the base level
• Coefficients are additive dollars: an inpatient stay adds $37,722 to the base; a Metro member adds $94.50 to every encounter, from a $55 lab test to a $30,000 admission

A normal distribution puts no probability on an exact zero, and is symmetric. Aggregate loss per policy is usually zero — and occasionally enormous.

The standard deviation is roughly proportional to the mean (a constant coefficient of variation) — nothing like OLS’s single σ².

$$\ln y = \beta_0 + \beta_1 x_1 + \dots + \varepsilon \;\Longrightarrow\; y = e^{\beta_0}\cdot e^{\beta_1 x_1}\cdots e^{\varepsilon}$$

e["log_allowed"] = np.log(e["allowed"])
loglin = smf.ols("log_allowed ~ " + rhs, data=e).fit()

rel = np.exp(loglin.params)     # relativities
loglin.rsquared                 # 0.397
np.sqrt(loglin.scale)           # residual sd 1.019

• e^β is a relativity: Metro costs 1.163× Suburban, for every kind of encounter
• This is exactly how rating plans work: base rate × factor × factor × …

Category relativities have no single “truth”: heavy users’ encounters are more intensive, and the mix differs by category.

$$e^{\,\mathrm{E}[\ln y]} \neq \mathrm{E}[y] \qquad \text{for lognormal errors: } \mathrm{E}[y] = e^{\mathbf{x}^\top\beta + \sigma^2/2}$$

np.exp(fitted) estimates the median, not the mean. Duan’s smearing estimator multiplies by the average of e^residual:

naive = np.exp(loglin.fittedvalues)
smear = np.mean(np.exp(loglin.resid))          # 1.755

# heteroscedastic? smear within each category
cat = (e.assign(r=np.exp(loglin.resid))
        .groupby("service_category")["r"].mean())
fixed = naive * e["service_category"].map(cat)

$$\widehat{\mathrm{Var}}(\hat\beta) = (\mathbf{X}^\top\mathbf{X})^{-1}\,\Big(\sum_i \omega_i\,\hat e_i^2\,\mathbf{x}_i\mathbf{x}_i^\top\Big)\,(\mathbf{X}^\top\mathbf{X})^{-1}$$

• HC0 (White): ω_i = 1 — each row’s own squared residual
• HC3: ω_i = 1 / (1 − h_ii)² — inflates high-leverage rows; better in small samples

hc3 = smf.ols("log_allowed ~ " + rhs, data=e).fit(cov_type="HC3")
np.exp(hc3.conf_int())    # robust 95% CIs for the relativities