portfolio_risk.py

"""Portfolio credit risk: expected loss, Basel IRB capital and a Monte Carloloss distribution for the whole card book, using the PDs from pd_model.py.     Expected loss (EL)  = PD x LGD x EAD, per account    EAD                 = drawn balance + CCF x undrawn limit    IRB capital (K)     = Basel II/III retail formula for qualifying revolving                          retail exposures, asset correlation R = 0.04    Monte Carlo         = one-factor Gaussian copula (the model behind the IRB                          formula), simulated with idiosyncratic default noise LGD and CCF are not in the data; they are stated assumptions with stress cases.The PD here is the model's probability of a missed payment in the followingmonth, which is a shorter horizon than the one-year PD Basel uses, so capitalfigures illustrate the method rather than a regulatory number."""import numpy as npimport pandas as pdfrom scipy.stats import norm from pd_model import STATUS_ORDER LGD = 0.80            # unsecured card balances: assumptionCCF = 0.50            # share of the undrawn limit drawn by default: assumptionR_QRRE = 0.04         # Basel asset correlation for qualifying revolving retailCONFIDENCE = 0.999    # Basel IRB confidence levelPD_FLOOR = 0.001      # Basel III floor for revolving retail PDSIMULATIONS = 20_000SEED = 7  def ead(features: pd.DataFrame, ccf: float = CCF) -> np.ndarray:    return features["drawn_balance"].to_numpy() + ccf * features["undrawn_limit"].to_numpy()  def irb_capital(pd_: np.ndarray, lgd: float, rho: float = R_QRRE, q: float = CONFIDENCE) -> np.ndarray:    """Capital per unit of EAD: loss at the q quantile of the systematic factor, less EL."""    conditional = norm.cdf((norm.ppf(pd_) + np.sqrt(rho) * norm.ppf(q)) / np.sqrt(1 - rho))    return lgd * (conditional - pd_)  def simulate_losses(pd_: np.ndarray, exposure: np.ndarray, lgd: float, rho: float,                    n: int = SIMULATIONS, seed: int = SEED, chunk: int = 250) -> np.ndarray:    """Portfolio loss in each scenario: an account defaults when its latent value    sqrt(rho) * Z + sqrt(1 - rho) * e falls below the PD threshold."""    rng = np.random.default_rng(seed)    threshold = norm.ppf(pd_).astype(np.float32)    loss_if_default = (lgd * exposure).astype(np.float32)    a, b = np.float32(np.sqrt(rho)), np.float32(np.sqrt(1 - rho))    losses = np.empty(n)    for start in range(0, n, chunk):        m = min(chunk, n - start)        z = rng.standard_normal((m, 1), dtype=np.float32)        e = rng.standard_normal((m, pd_.size), dtype=np.float32)        defaults = (a * z + b * e) < threshold        losses[start:start + m] = defaults @ loss_if_default    return losses  def tail(losses: np.ndarray, q: float) -> tuple[float, float]:    var = float(np.quantile(losses, q))    return var, float(losses[losses >= var].mean())  def scenario(name: str, pd_: np.ndarray, exposure: np.ndarray, lgd: float, rho: float, sims: int) -> dict:    el = float(np.sum(pd_ * lgd * exposure))    k = float(np.sum(irb_capital(pd_, lgd, rho) * exposure))    losses = simulate_losses(pd_, exposure, lgd, rho, n=sims)    var999, es999 = tail(losses, CONFIDENCE)    return {"scenario": name, "expected_loss": el, "irb_capital": k, "var_999": var999,            "unexpected_loss_999": var999 - float(losses.mean()), "es_999": es999,            "simulated_mean": float(losses.mean()), "losses": losses}  def run(features: pd.DataFrame, pd_model_values: np.ndarray) -> dict:    pd_ = np.clip(pd_model_values, PD_FLOOR, 0.9999)    exposure = ead(features)    total_ead = float(exposure.sum())    el_account = pd_ * LGD * exposure     base = scenario("Base", pd_, exposure, LGD, R_QRRE, SIMULATIONS)    stresses = [        scenario("PD x 1.5", np.clip(pd_ * 1.5, PD_FLOOR, 0.9999), exposure, LGD, R_QRRE, 10_000),        scenario("LGD 90%", pd_, exposure, 0.90, R_QRRE, 10_000),        scenario("CCF 100% (full limit drawn)", pd_, ead(features, 1.0), LGD, R_QRRE, 10_000),        scenario("Correlation 0.15", pd_, exposure, LGD, 0.15, 10_000),        scenario("Combined: PD x 1.5, LGD 90%, correlation 0.15",                 np.clip(pd_ * 1.5, PD_FLOOR, 0.9999), exposure, 0.90, 0.15, 10_000),    ]     # Loss distribution for the chart, as a share of total EAD.    share = base["losses"] / total_ead    counts, edges = np.histogram(share, bins=40)    histogram = [{"from": float(edges[i]), "to": float(edges[i + 1]), "share_of_sims": float(counts[i] / share.size)}                 for i in range(len(counts))]     # Where the expected loss sits: by most recent repayment status.    by_status = (pd.DataFrame({"segment": features["status_sep_group"], "ead": exposure, "el": el_account})                 .groupby("segment").sum().reindex(STATUS_ORDER))    concentration = [{"segment": seg, "share_of_ead": float(r["ead"] / total_ead),                      "share_of_el": float(r["el"] / el_account.sum())} for seg, r in by_status.iterrows()]     def strip(s: dict) -> dict:        return {k: v for k, v in s.items() if k != "losses"}     return {        "assumptions": {"lgd": LGD, "ccf": CCF, "correlation": R_QRRE, "confidence": CONFIDENCE,                        "pd_floor": PD_FLOOR, "simulations": SIMULATIONS, "currency": "NT$"},        "book": {"accounts": int(len(features)), "credit_limits": float(features["credit_limit"].sum()),                 "drawn": float(features["drawn_balance"].sum()), "ead": total_ead,                 "average_pd": float(pd_.mean()), "ead_weighted_pd": float(np.sum(pd_ * exposure) / total_ead)},        "base": strip(base),        "stress": [strip(s) for s in stresses],        "histogram": histogram,        "concentration": concentration,        "account_level": pd.DataFrame({"pd": pd_, "ead": exposure, "lgd": LGD}),    }