import os
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ.setdefault(_v, "2")   # (later sections use LightGBM; cap threads — energy-notebook lesson)
import numpy as np, pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
plt.rcParams.update({"figure.dpi": 110, "font.size": 10, "axes.grid": True,
                     "grid.alpha": 0.25, "axes.spines.top": False, "axes.spines.right": False})
DATA = Path("data")

# transaction is 683 MB / 394 cols; load a focused, analysis-relevant subset and downcast to save memory
TCOLS = (["TransactionID","isFraud","TransactionDT","TransactionAmt","ProductCD",
          "card1","card2","card3","card4","card5","card6","addr1","addr2"]
         + [f"D{i}" for i in range(1,16)] + [f"C{i}" for i in range(1,15)] + [f"M{i}" for i in range(1,10)])
tx = pd.read_csv(DATA/"train_transaction.csv", usecols=TCOLS)
idf = pd.read_csv(DATA/"train_identity.csv")
df = tx.merge(idf, on="TransactionID", how="left")           # left join: not every txn has identity
for c in df.select_dtypes("float64").columns: df[c] = df[c].astype("float32")

df["day"] = (df["TransactionDT"]//86400).astype(int)          # day index (relative)
N=len(df); F=int(df.isFraud.sum())
print(f"transactions {N:,} | fraud {F:,} ({100*F/N:.2f}%) | columns loaded {df.shape[1]}")
print(f"time span: day {df.day.min()} .. {df.day.max()} ({df.day.max()-df.day.min()+1} days)")
print(f"with identity records: {idf.shape[0]:,} ({100*idf.shape[0]/N:.0f}% of transactions)")

transactions 590,540 | fraud 20,663 (3.50%) | columns loaded 92
time span: day 1 .. 182 (182 days)
with identity records: 144,233 (24% of transactions)

# EDA 1 — imbalance and fraud rate by product segment
fig, ax = plt.subplots(1, 2, figsize=(12, 3.8))
ax[0].bar(["legit","fraud"], [N-F, F], color=["#2b6cb0","#c53030"]); ax[0].set_yscale("log")
ax[0].set_ylabel("transactions (log)"); ax[0].set_title(f"Imbalance — fraud {100*F/N:.2f}%")
for i,v in enumerate([N-F,F]): ax[0].text(i,v,f"{v:,}",ha="center",va="bottom",fontsize=9)

pr = df.groupby("ProductCD")["isFraud"].agg(["mean","size"]).sort_values("mean")
ax[1].bar(pr.index, pr["mean"]*100, color="#c53030")
ax[1].axhline(100*F/N, color="grey", ls=":", label=f"overall {100*F/N:.2f}%")
ax[1].set_ylabel("fraud rate (%)"); ax[1].set_xlabel("ProductCD"); ax[1].set_title("Fraud rate varies sharply by product")
ax[1].legend(fontsize=8)
for i,(p,r) in enumerate(zip(pr.index, pr["mean"]*100)): ax[1].text(i, r, f"{r:.1f}%", ha="center", va="bottom", fontsize=8)
plt.tight_layout(); plt.show()
print("ProductCD fraud rates:", {k: f"{v*100:.1f}%" for k,v in pr["mean"].items()})
print("Product 'C' (and 'S') carry far higher fraud than the dominant 'W' — a real, usable segment, unlike")
print("anything ULB's anonymized features could give us.")

ProductCD fraud rates: {'W': '2.0%', 'R': '3.8%', 'H': '4.8%', 'S': '5.9%', 'C': '11.7%'}
Product 'C' (and 'S') carry far higher fraud than the dominant 'W' — a real, usable segment, unlike
anything ULB's anonymized features could give us.

# EDA 2 — the missingness landscape (structured, and informative)
groups = {"D (timedeltas)": [f"D{i}" for i in range(1,16)],
          "C (counts)":     [f"C{i}" for i in range(1,15)],
          "M (match flags)":[f"M{i}" for i in range(1,10)],
          "id_ (identity)": [c for c in df.columns if c.startswith("id_")]}
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
miss = df[[c for g in groups.values() for c in g]].isna().mean()*100
allc = [c for g in groups.values() for c in g]
colmap = {c: col for g,col in zip(groups, ["#c53030","#2b6cb0","#2f855a","#dd6b20"]) for c in groups[g]}
ax[0].bar(range(len(allc)), miss.values, color=[colmap[c] for c in allc])
ax[0].set_xticks(range(len(allc))); ax[0].set_xticklabels(allc, rotation=90, fontsize=5)
ax[0].set_ylabel("% missing"); ax[0].set_title("Missingness by column (color = group)")

# is missingness itself predictive? fraud rate when D1 present vs missing, etc.
sig = {}
for c in ["D1","D2","D6","D11","DeviceType","id_01"]:
    if c in df: sig[c] = [df.loc[df[c].notna(),"isFraud"].mean()*100, df.loc[df[c].isna(),"isFraud"].mean()*100]
sg = pd.DataFrame(sig, index=["present","missing"]).T
sg.plot.bar(ax=ax[1], color=["#2b6cb0","#c53030"])
ax[1].axhline(100*F/N, color="grey", ls=":", label=f"overall {100*F/N:.2f}%")
ax[1].set_ylabel("fraud rate (%)"); ax[1].set_title("Missingness is informative (fraud rate present vs missing)")
ax[1].legend(fontsize=8); ax[1].tick_params(axis="x", labelrotation=0)
plt.tight_layout(); plt.show()
print("Whether a field is even present shifts the fraud rate — so missingness is signal, not just noise")
print("(we will keep NaNs as-is for the tree model in §5, which handles them natively).")

Whether a field is even present shifts the fraud rate — so missingness is signal, not just noise
(we will keep NaNs as-is for the tree model in §5, which handles them natively).

# EDA 3 — fraud through the 182-day window, and amount
fig, ax = plt.subplots(1, 2, figsize=(12, 3.8))
by_day = df.groupby("day")["isFraud"].mean()*100
ax[0].plot(by_day.index, by_day.values, color="#c53030", lw=0.9)
ax[0].axhline(100*F/N, color="grey", ls=":", label=f"overall {100*F/N:.2f}%")
ax[0].set_xlabel("day"); ax[0].set_ylabel("fraud rate (%)"); ax[0].set_title("Fraud rate over 182 days")
ax[0].legend(fontsize=8)
bins = np.logspace(0, np.log10(df.TransactionAmt.max()+1), 50)
ax[1].hist(df.loc[df.isFraud==0,"TransactionAmt"].clip(lower=1), bins=bins, density=True, alpha=0.6, color="#2b6cb0", label="legit")
ax[1].hist(df.loc[df.isFraud==1,"TransactionAmt"].clip(lower=1), bins=bins, density=True, alpha=0.6, color="#c53030", label="fraud")
ax[1].set_xscale("log"); ax[1].set_xlabel("amount ($, log)"); ax[1].set_ylabel("density")
ax[1].set_title("Transaction amount: fraud vs legit"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

SPLIT_DAY = int(df.day.quantile(0.80))
print(f"time-ordered split at day {SPLIT_DAY}: train days {df.day.min()}-{SPLIT_DAY} "
      f"({int((df.day<=SPLIT_DAY).sum()):,} txns, {int(df.loc[df.day<=SPLIT_DAY,'isFraud'].sum())} fraud) | "
      f"test {SPLIT_DAY+1}-{df.day.max()} ({int((df.day>SPLIT_DAY).sum()):,} txns, "
      f"{int(df.loc[df.day>SPLIT_DAY,'isFraud'].sum())} fraud)")

time-ordered split at day 141: train days 1-141 (475,006 txns, 16721 fraud) | test 142-182 (115,534 txns, 3942 fraud)

# --- build the anchor and show the raw D1 drifts while the anchor does not ----------------------------
df["D1_anchor"] = df["day"] - df["D1"]          # calendar day the card was first seen (fixed per card)
tr_mask = df["day"] <= SPLIT_DAY

fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot(df.groupby("day")["D1"].mean(), color="#c53030", label="raw D1 (days-since-first-seen)")
ax[0].plot(df.groupby("day")["D1_anchor"].mean(), color="#2b6cb0", label="D1_anchor = day - D1")
ax[0].axvline(SPLIT_DAY, color="grey", ls=":", label="train/test split")
ax[0].set_xlabel("day"); ax[0].set_ylabel("mean value"); ax[0].set_title("Raw D1 drifts upward; the anchor is stable")
ax[0].legend(fontsize=8)

# train vs test distribution: raw D1 shifts, anchor overlaps
b1 = np.linspace(0, 640, 60)
ax[1].hist(df.loc[tr_mask,"D1"].dropna(), bins=b1, density=True, alpha=0.5, color="#2b6cb0", label="D1 train")
ax[1].hist(df.loc[~tr_mask,"D1"].dropna(), bins=b1, density=True, alpha=0.5, color="#c53030", label="D1 test")
ax[1].set_xlabel("raw D1 (days)"); ax[1].set_ylabel("density"); ax[1].set_title("Raw D1: train vs test distribution shift")
ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()

# quantify the drift with a simple mean shift
print(f"raw D1  mean: train {df.loc[tr_mask,'D1'].mean():.1f}  ->  test {df.loc[~tr_mask,'D1'].mean():.1f}  "
      f"(shift {df.loc[~tr_mask,'D1'].mean()-df.loc[tr_mask,'D1'].mean():+.1f} days)")
print("The raw feature literally means something different in the test period — the anchor removes that by")
print("re-expressing it as a fixed calendar date instead of a moving 'days ago'.")

raw D1  mean: train 90.8  ->  test 109.1  (shift +18.3 days)
The raw feature literally means something different in the test period — the anchor removes that by
re-expressing it as a fixed calendar date instead of a moving 'days ago'.

# --- the payoff: the anchor turns transactions into recognizable card entities -----------------------
# a lightweight UID from card identity + the anchor (full construction in §3)
uid0 = (df["card1"].astype("Int64").astype(str) + "_" + df["addr1"].astype("Int64").astype(str)
        + "_" + df["D1_anchor"].round().astype("Int64").astype(str))
uid0 = uid0.where(df[["card1","addr1","D1_anchor"]].notna().all(axis=1))     # only where keys present
sizes = uid0.value_counts()                       # already sorted descending
multi = sizes[sizes > 1]
n_txn_multi = int(uid0.isin(multi.index).sum())

# example: the busiest card entity, its repeated transactions across days
ex = multi.index[0]
cols_show = ["day","TransactionAmt","ProductCD","isFraud"]
exdf = df.loc[uid0==ex, cols_show].sort_values("day").head(8)

fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].hist(sizes.values, bins=np.arange(1, 30), color="#2b6cb0", align="left")
ax[0].set_yscale("log"); ax[0].set_xlabel("transactions per UID (card1+addr1+anchor)")
ax[0].set_ylabel("number of UIDs (log)"); ax[0].set_title("Anchor + keys resolve repeat cards")
n_txn_multi = int((uid0.isin(multi.index)).sum())
ax[0].text(0.5, 0.9, f"{100*n_txn_multi/len(df):.0f}% of transactions belong\nto a UID seen >1 time",
           transform=ax[0].transAxes, fontsize=9, va="top")

ax[1].axis("off"); ax[1].text(0.0, 1.0, f"Example UID (one card): {ex}", fontsize=9, family="monospace",
                              va="top", transform=ax[1].transAxes)
ax[1].text(0.0, 0.86, exdf.to_string(index=False), fontsize=9, family="monospace", va="top",
           transform=ax[1].transAxes)
plt.tight_layout(); plt.show()

print(f"unique cards by card1 alone: {df['card1'].nunique():,}")
print(f"unique UIDs by card1+addr1+anchor: {sizes.shape[0]:,}  |  {100*n_txn_multi/len(df):.0f}% of txns are repeat-UID")
print("The anchor is what makes the repeats VISIBLE: the same card, transacting across many days, now falls")
print("into one group — the entity §3 formalizes and §4 turns into velocity / recency / history features.")

unique cards by card1 alone: 13,553
unique UIDs by card1+addr1+anchor: 199,070  |  69% of txns are repeat-UID
The anchor is what makes the repeats VISIBLE: the same card, transacting across many days, now falls
into one group — the entity §3 formalizes and §4 turns into velocity / recency / history features.

# --- construct the UID and measure coverage ----------------------------------------------------------
def key(s): return s.astype("Int64").astype(str)
df["uid"] = (key(df["card1"]) + "_" + key(df["addr1"]) + "_" + df["D1_anchor"].round().astype("Int64").astype(str))
df.loc[df[["card1","addr1","D1_anchor"]].isna().any(axis=1), "uid"] = np.nan   # unresolved if any key missing

resolved = df["uid"].notna()
sizes = df["uid"].value_counts()
multi = sizes[sizes > 1].index
in_multi = df["uid"].isin(multi)
print(f"resolved (all keys present): {100*resolved.mean():.0f}% of transactions")
print(f"distinct UIDs: {sizes.shape[0]:,} | transactions in a repeat UID: {100*in_multi.mean():.0f}%")
print(f"UID size: median {int(sizes.median())}, 95th pct {int(sizes.quantile(.95))}, max {int(sizes.max())}")

# persistence across the train/test split: entities whose history a test txn could use
uid_tr = set(df.loc[df.day<=SPLIT_DAY, "uid"].dropna().unique())
test_res = df.loc[(df.day>SPLIT_DAY) & resolved]
seen_before = test_res["uid"].isin(uid_tr)
print(f"\nof resolved TEST transactions, {100*seen_before.mean():.0f}% are on a card already seen in TRAIN")
print(f"  fraud rate — card seen in train: {100*test_res.loc[seen_before,'isFraud'].mean():.2f}%  vs  "
      f"new card: {100*test_res.loc[~seen_before,'isFraud'].mean():.2f}%")

resolved (all keys present): 89% of transactions
distinct UIDs: 199,070 | transactions in a repeat UID: 69%
UID size: median 1, 95th pct 9, max 1414

of resolved TEST transactions, 46% are on a card already seen in TRAIN
  fraud rate — card seen in train: 1.59%  vs  new card: 2.87%

# --- validate the UID: do grouped transactions share consistent card attributes? ---------------------
g = df.loc[in_multi].groupby("uid")
consist = {c: (g[c].nunique(dropna=True) <= 1).mean()*100 for c in ["card4","card6","ProductCD","card3","card5"]}
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
ax[0].bar(list(consist), list(consist.values()), color="#2f855a")
ax[0].set_ylabel("% of repeat-UIDs with a single value"); ax[0].set_ylim(0,105)
ax[0].set_title("UID validity: attribute consistency within an entity")
for i,(k,v) in enumerate(consist.items()): ax[0].text(i, v, f"{v:.0f}%", ha="center", va="bottom", fontsize=9)

# fraud concentration: cumulative fraud captured by the most fraud-heavy UIDs
fr_by_uid = df.loc[resolved].groupby("uid")["isFraud"].sum().sort_values(ascending=False)
cum = np.cumsum(fr_by_uid.values) / fr_by_uid.values.sum()
ax[1].plot(np.arange(1, len(cum)+1)/len(cum)*100, cum*100, color="#c53030")
ax[1].set_xlabel("% of UIDs (ranked by fraud count)"); ax[1].set_ylabel("% of all fraud captured")
ax[1].set_title("Fraud concentrates in a small set of card entities")
frac_uids_for_80 = 100*np.searchsorted(cum, 0.80)/len(cum)
ax[1].axhline(80, color="grey", ls=":"); ax[1].annotate(f"{frac_uids_for_80:.1f}% of UIDs hold 80% of fraud",
             (frac_uids_for_80, 80), fontsize=8, xytext=(frac_uids_for_80+10, 60),
             arrowprops=dict(arrowstyle="->"))
plt.tight_layout(); plt.show()
print("Card attributes are highly consistent within a UID (brand/type rarely vary) — evidence the UID")
print(f"really is grouping one card. And fraud is concentrated: ~{frac_uids_for_80:.0f}% of entities carry")
print("80% of all fraud — a strong signal that per-entity features (§4) should capture.")

Card attributes are highly consistent within a UID (brand/type rarely vary) — evidence the UID
really is grouping one card. And fraud is concentrated: ~1% of entities carry
80% of all fraud — a strong signal that per-entity features (§4) should capture.

# --- leak-safe per-UID aggregates (sorted by time; only PRIOR transactions of the same card) ---------
df = df.sort_values("TransactionDT").reset_index(drop=True)
g = df.groupby("uid")
df["uid_n_prior"]        = g.cumcount()                                   # 0 for a card's first txn
df["uid_amt_prior_sum"]  = g["TransactionAmt"].cumsum() - df["TransactionAmt"]
df["uid_amt_prior_mean"] = df["uid_amt_prior_sum"] / df["uid_n_prior"].where(df["uid_n_prior"] > 0)
df["uid_days_since_prev"]= df["day"] - g["day"].shift(1)                  # recency (NaN on first sighting)
df["amt_dev"]            = df["TransactionAmt"] - df["uid_amt_prior_mean"]# amount vs the card's own history
# rows with no UID (missing keys) keep NaN aggregates — fine for the tree model
df.loc[df["uid"].isna(), ["uid_n_prior","uid_amt_prior_mean","uid_days_since_prev","amt_dev"]] = np.nan
print("engineered features:", ["uid_n_prior","uid_amt_prior_mean","uid_days_since_prev","amt_dev"])
print(f"coverage (non-null uid_days_since_prev = repeat sightings): {100*df['uid_days_since_prev'].notna().mean():.0f}% of txns")

engineered features: ['uid_n_prior', 'uid_amt_prior_mean', 'uid_days_since_prev', 'amt_dev']
coverage (non-null uid_days_since_prev = repeat sightings): 55% of txns

# --- do these behavioral features separate fraud? ----------------------------------------------------
base = df["isFraud"].mean()*100
fig, ax = plt.subplots(1, 3, figsize=(14, 4))

# (1) history: fraud rate by number of prior sightings
nb_bucket = pd.cut(df["uid_n_prior"], [-1,0,1,2,5,10,1e9], labels=["0","1","2","3-5","6-10","10+"])
fr1 = df.groupby(nb_bucket, observed=True)["isFraud"].mean()*100
ax[0].bar(range(len(fr1)), fr1.values, color="#2b6cb0"); ax[0].set_xticks(range(len(fr1))); ax[0].set_xticklabels(fr1.index)
ax[0].axhline(base, color="grey", ls=":"); ax[0].set_xlabel("prior sightings of this card")
ax[0].set_ylabel("fraud rate (%)"); ax[0].set_title("History: fewer priors -> more fraud")

# (2) recency: fraud rate by days since previous txn
rb = pd.cut(df["uid_days_since_prev"], [-1,0,1,3,7,30,1e9], labels=["same day","1","2-3","4-7","8-30","30+"])
fr2 = df.groupby(rb, observed=True)["isFraud"].mean()*100
ax[1].bar(range(len(fr2)), fr2.values, color="#dd6b20"); ax[1].set_xticks(range(len(fr2))); ax[1].set_xticklabels(fr2.index, fontsize=8)
ax[1].axhline(base, color="grey", ls=":"); ax[1].set_xlabel("days since card's previous txn")
ax[1].set_ylabel("fraud rate (%)"); ax[1].set_title("Recency: bursts of activity")

# (3) deviation: fraud rate when the amount is out of character
dev_bucket = pd.qcut(df["amt_dev"].dropna(), 6, duplicates="drop")
fr3 = df.loc[df["amt_dev"].notna()].groupby(dev_bucket, observed=True)["isFraud"].mean()*100
ax[2].plot(range(len(fr3)), fr3.values, "o-", color="#c53030"); ax[2].axhline(base, color="grey", ls=":")
ax[2].set_xticks(range(len(fr3))); ax[2].set_xticklabels([f"{iv.mid:.0f}" for iv in fr3.index], rotation=45, fontsize=7)
ax[2].set_xlabel("amount - card's prior mean ($)"); ax[2].set_ylabel("fraud rate (%)")
ax[2].set_title("Deviation: charges out of character")
plt.tight_layout(); plt.show()

print(f"overall fraud rate {base:.2f}%. All three signals are leak-safe (prior transactions only):")
print(f"  history  — fraud RISES with prior sightings: {fr1.iloc[0]:.1f}% at first sighting -> {fr1.iloc[-1]:.1f}% at 10+ "
      f"(fraud arrives in bursts on a compromised card).")
print(f"  recency  — the sharpest signal: a repeat within ~1 day runs {fr2.max():.1f}% vs {fr2.min():.1f}% for a "
      f"long-dormant card (~{fr2.max()/fr2.min():.0f}x).")
print(f"  deviation— charges far from a card's own average, in EITHER direction, are elevated vs its norm.")
print("None of these is computable from a single transaction in isolation.")

overall fraud rate 3.50%. All three signals are leak-safe (prior transactions only):
  history  — fraud RISES with prior sightings: 2.0% at first sighting -> 3.8% at 10+ (fraud arrives in bursts on a compromised card).
  recency  — the sharpest signal: a repeat within ~1 day runs 6.1% vs 0.7% for a long-dormant card (~9x).
  deviation— charges far from a card's own average, in EITHER direction, are elevated vs its norm.
None of these is computable from a single transaction in isolation.

import lightgbm as lgb
from sklearn.metrics import average_precision_score, roc_auc_score

ENG = ["D1_anchor", "uid_n_prior", "uid_amt_prior_mean", "uid_days_since_prev", "amt_dev"]
DROP = {"TransactionID","isFraud","TransactionDT","day","uid","uid_amt_prior_sum", *ENG}
RAW = [c for c in df.columns if c not in DROP]
VENDOR = [c for c in RAW if len(c)>1 and c[0] in "CD" and c[1:].isdigit()]   # C1-14, D1-15 = Vesta's pre-built aggregates
MINIMAL = [c for c in RAW if c not in VENDOR]                                 # strip them to expose our features' value

X = df.copy()
for c in X.columns:
    if X[c].dtype == "object": X[c] = X[c].astype("category")
tr = (df.day <= SPLIT_DAY).values; te = (df.day > SPLIT_DAY).values; y = df["isFraud"].values

def fit_eval(feats, seed=0):
    catf = [c for c in feats if str(X[c].dtype) == "category"]
    m = lgb.LGBMClassifier(n_estimators=300, learning_rate=0.05, num_leaves=64, min_child_samples=50,
                           subsample=0.8, colsample_bytree=0.8, reg_lambda=5.0, random_state=seed, verbose=-1)
    m.fit(X.loc[tr, feats], y[tr], categorical_feature=catf)
    s = m.predict_proba(X.loc[te, feats])[:, 1]
    return m, average_precision_score(y[te], s), roc_auc_score(y[te], s)

runs = {"minimal (no vendor C/D)": MINIMAL, "minimal + ours": MINIMAL+ENG,
        "full raw (with vendor C/D)": RAW, "full raw + ours": RAW+ENG}
res = {}
for name, feats in runs.items():
    m, ap, auc = fit_eval(feats); res[name] = (m, ap, auc)
    print(f"  {name:28s} PR-AUC {ap:.4f}  ROC {auc:.4f}  ({len(feats)} feats)")
m_eng = res["full raw + ours"][0]

lift_min  = res["minimal + ours"][1]  - res["minimal (no vendor C/D)"][1]
lift_full = res["full raw + ours"][1] - res["full raw (with vendor C/D)"][1]
print(f"\nlift of our features on MINIMAL data:  {lift_min:+.4f} PR-AUC  (they stand in for the missing aggregates)")
print(f"lift of our features on FULL raw data: {lift_full:+.4f} PR-AUC  (mostly redundant with Vesta's C/D)")

fig, ax = plt.subplots(figsize=(9, 4))
names = list(runs); aps = [res[n][1] for n in names]
ax.bar(names, aps, color=["#a0aec0","#2f855a","#2b6cb0","#c53030"])
for i,v in enumerate(aps): ax.text(i, v, f"{v:.3f}", ha="center", va="bottom", fontsize=9)
ax.set_ylabel("out-of-time PR-AUC"); ax.set_title("Feature-engineering lift, isolated by baseline")
ax.tick_params(axis="x", labelsize=8, labelrotation=10)
plt.tight_layout(); plt.show()

  minimal (no vendor C/D)      PR-AUC 0.3202  ROC 0.8376  (59 feats)

  minimal + ours               PR-AUC 0.3505  ROC 0.8661  (64 feats)

  full raw (with vendor C/D)   PR-AUC 0.5375  ROC 0.9106  (88 feats)

  full raw + ours              PR-AUC 0.5397  ROC 0.9113  (93 feats)

lift of our features on MINIMAL data:  +0.0303 PR-AUC  (they stand in for the missing aggregates)
lift of our features on FULL raw data: +0.0022 PR-AUC  (mostly redundant with Vesta's C/D)

# --- SHAP (TreeSHAP, exact via LightGBM pred_contrib): which features drive the model? ----------------
rng = np.random.default_rng(0)
samp = rng.choice(np.where(te)[0], size=min(20000, int(te.sum())), replace=False)
feats = RAW + ENG
contrib = m_eng.booster_.predict(X.iloc[samp][feats], pred_contrib=True)   # (n, n_feat+1); last col = base
shap_vals = contrib[:, :-1]
imp = pd.Series(np.abs(shap_vals).mean(0), index=feats).sort_values(ascending=False)
top = imp.head(16)

fig, ax = plt.subplots(1, 2, figsize=(13, 5))
colors = ["#c53030" if f in ENG else "#2b6cb0" for f in top.index[::-1]]
ax[0].barh(range(len(top)), top.values[::-1], color=colors)
ax[0].set_yticks(range(len(top))); ax[0].set_yticklabels(top.index[::-1], fontsize=8)
ax[0].set_xlabel("mean |SHAP| (impact on log-odds)"); ax[0].set_title("Top features by SHAP (red = engineered)")

# SHAP dependence for the top engineered feature: how its value drives risk
best_eng = [f for f in imp.index if f in ENG][0]
j = feats.index(best_eng); xv = X.iloc[samp][best_eng].astype(float).values
ax[1].scatter(xv, shap_vals[:, j], s=5, alpha=0.2, color="#c53030")
ax[1].axhline(0, color="grey", lw=0.8); ax[1].set_xlabel(f"{best_eng} (value)")
ax[1].set_ylabel("SHAP value (→ pushes toward fraud)"); ax[1].set_title(f"SHAP dependence: {best_eng}")
if best_eng in ("amt_dev","uid_amt_prior_mean"): ax[1].set_xscale("symlog")
plt.tight_layout(); plt.show()

eng_ranks = {f: list(imp.index).index(f)+1 for f in ENG}
print(f"engineered feature SHAP ranks (of {len(imp)}): {eng_ranks}")
print(f"D1_anchor ranks #{eng_ranks['D1_anchor']} even with raw D1 present — the model reaches for the anchor")
print("transform directly. The UID aggregates rank lower HERE only because Vesta's C/D columns already encode")
print("similar entity signal (hence the small full-data lift); strip those and our features carry the load.")

engineered feature SHAP ranks (of 93): {'D1_anchor': 2, 'uid_n_prior': 23, 'uid_amt_prior_mean': 39, 'uid_days_since_prev': 46, 'amt_dev': 52}
D1_anchor ranks #2 even with raw D1 present — the model reaches for the anchor
transform directly. The UID aggregates rank lower HERE only because Vesta's C/D columns already encode
similar entity signal (hence the small full-data lift); strip those and our features carry the load.

import subprocess, sys
from sklearn.isotonic import IsotonicRegression
from sklearn.calibration import calibration_curve
from sklearn.metrics import brier_score_loss

# 3-way time split; refit the engineered score leak-safely on train only
D_TR, D_CA = 120, 141
trn = (df.day <= D_TR).values; cal = ((df.day > D_TR) & (df.day <= D_CA)).values; tst = (df.day > D_CA).values
feats = RAW + ENG; catf = [c for c in feats if str(X[c].dtype) == "category"]
gbm7 = lgb.LGBMClassifier(n_estimators=300, learning_rate=0.05, num_leaves=64, min_child_samples=50,
                          subsample=0.8, colsample_bytree=0.8, reg_lambda=5.0, verbose=-1)
gbm7.fit(X.loc[trn, feats], y[trn], categorical_feature=catf)
s_ca = gbm7.predict_proba(X.loc[cal, feats])[:, 1]; s_te = gbm7.predict_proba(X.loc[tst, feats])[:, 1]
yca, yte = y[cal], y[tst]
amt_ca = df.loc[cal, "TransactionAmt"].values; amt_te = df.loc[tst, "TransactionAmt"].values
clip = lambda s: np.clip(s, 1e-6, 1-1e-6); z_ca = np.log(clip(s_ca)/(1-clip(s_ca))); z_te = np.log(clip(s_te)/(1-clip(s_te)))

# REAL segments: ProductCD x card type x device
segkey = (df["ProductCD"].astype(str)+"|"+df["card6"].astype(str)+"|"+df["DeviceType"].astype(str))
seg_all, uniq = pd.factorize(segkey); df["seg"] = seg_all; G = len(uniq)
seg_ca = df.loc[cal, "seg"].values; seg_te = df.loc[tst, "seg"].values
print(f"3-way split — train {trn.sum():,} / calib {cal.sum():,} ({yca.sum()} fraud) / test {tst.sum():,} ({yte.sum()} fraud)")
print(f"real segments (ProductCD x card6 x DeviceType): {G}")

def ece(yv, p, nb=10):
    q = np.quantile(p, np.linspace(0,1,nb+1)); q[0]-=1e-9; b = np.digitize(p, q[1:-1]); e=0.0
    for k in range(nb):
        m=b==k
        if m.sum(): e += abs(p[m].mean()-yv[m].mean())*m.sum()
    return e/len(yv)

iso = IsotonicRegression(out_of_bounds="clip").fit(s_ca, yca); p_iso = iso.transform(s_te)
np.savez("calib_input.npz", z_ca=z_ca, seg_ca=seg_ca, yca=yca, G=G)   # reuse the ULB companion fitter
r = subprocess.run([sys.executable, "calib_fit.py"], capture_output=True, text=True)
print(r.stdout.strip() or r.stderr[-300:])
P = np.load("calib_post.npz"); a_d, b_d = P["a"], P["b"]
p_draws = 1/(1+np.exp(-(a_d[:, seg_te] + b_d[:, seg_te]*z_te[None,:])))
p_bayes = p_draws.mean(0); p_lo, p_hi = np.percentile(p_draws,5,0), np.percentile(p_draws,95,0)

print("\nCalibration on TEST (x1e-3, lower=better):")
for nm,p in [("raw score", s_te), ("isotonic", p_iso), ("Bayes-hier (real segments)", p_bayes)]:
    print(f"  {nm:28s} Brier {brier_score_loss(yte, np.clip(p,0,1))*1e3:.3f}   ECE {ece(yte, np.clip(p,0,1))*1e3:.3f}")

3-way split — train 414,542 / calib 60,464 (2121 fraud) / test 115,534 (3942 fraud)
real segments (ProductCD x card6 x DeviceType): 39

calib_post.npz saved | final ELBO loss 4911.52001953125

Calibration on TEST (x1e-3, lower=better):
  raw score                    Brier 22.566   ECE 3.039
  isotonic                     Brier 22.765   ECE 4.081
  Bayes-hier (real segments)   Brier 22.861   ECE 3.587

# reliability + does pooling help the THIN segments this time?
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot([0,1],[0,1],"k:",label="perfect")
for nm,p,col in [("raw",s_te,"#a0aec0"),("isotonic",p_iso,"#2f855a"),("Bayes-hier",p_bayes,"#2b6cb0")]:
    fp,mp = calibration_curve(yte, np.clip(p,0,1), n_bins=8, strategy="quantile"); ax[0].plot(mp,fp,"o-",color=col,label=nm)
ax[0].set_xlabel("predicted P(fraud)"); ax[0].set_ylabel("observed"); ax[0].set_title("Reliability (real-segment calibration)")
ax[0].legend(fontsize=8)

fps = pd.Series(yca).groupby(seg_ca).sum()
thin = set(fps[fps<=5].index); mask = np.isin(seg_te, list(thin))
bars = {"isotonic": ece(yte[mask], np.clip(p_iso[mask],0,1))*1e3, "Bayes-hier": ece(yte[mask], np.clip(p_bayes[mask],0,1))*1e3}
ax[1].bar(list(bars), list(bars.values()), color=["#2f855a","#2b6cb0"])
ax[1].set_ylabel("ECE on thin segments (x1e-3)"); ax[1].set_title(f"Thin segments (<=5 calib frauds): {mask.sum():,} test txns")
for i,v in enumerate(bars.values()): ax[1].text(i,v,f"{v:.3f}",ha="center",va="bottom",fontsize=9)
plt.tight_layout(); plt.show()
print(f"thin segments: {len(thin)} of {G}, {mask.sum():,} test txns, {int(yte[mask].sum())} frauds")
print("Whether the hierarchy beats isotonic on these real, data-poorer segments is the honest question ULB")
print("could not pose (its features were anonymized) — the bars above answer it for this data.")

thin segments: 8 of 39, 1,257 test txns, 49 frauds
Whether the hierarchy beats isotonic on these real, data-poorer segments is the honest question ULB
could not pose (its features were anonymized) — the bars above answer it for this data.

C_FP, K_FRAUD = 10.0, 10.0
CFN = lambda a: a + K_FRAUD
tau_star = lambda a: C_FP / (C_FP + CFN(a))
def rcost(flag, yv, amt): return float((flag*(1-yv)*C_FP + (1-flag)*yv*CFN(amt)).sum())

# calibrated posterior mean on calib (to tune the global-threshold baseline without leakage)
p_ca = (1/(1+np.exp(-(a_d[:, seg_ca] + b_d[:, seg_ca]*z_ca[None,:])))).mean(0)
grid = np.linspace(0.001, 0.9, 400)
tau_g = grid[np.argmin([rcost((p_ca>t).astype(int), yca, amt_ca) for t in grid])]

policies = {"flat 0.50": (p_bayes>0.5).astype(int),
            f"global tau={tau_g:.3f}": (p_bayes>tau_g).astype(int),
            "amount-varying tau*": (p_bayes>tau_star(amt_te)).astype(int)}
print(f"cost model: false decline ${C_FP:.0f}; approved fraud = amount + ${K_FRAUD:.0f}\n")
print(f"{'policy':22s}{'test $ cost':>12}{'fraud $ exposed':>16}{'caught':>9}{'false decl':>12}")
for nm, flag in policies.items():
    caught=int(((flag==1)&(yte==1)).sum()); fd=int(((flag==1)&(yte==0)).sum())
    exp=float(((1-flag)*yte*amt_te).sum())
    print(f"{nm:22s}{rcost(flag,yte,amt_te):>12,.0f}{exp:>16,.0f}{caught:>6}/{int(yte.sum())}{fd:>12,}")

cost model: false decline $10; approved fraud = amount + $10

policy                 test $ cost fraud $ exposed   caught  false decl
flat 0.50                  444,385         413,155  1372/3942         553
global tau=0.048           265,395         140,035  2934/3942      11,528
amount-varying tau*        218,217         111,497  2540/3942       9,270

# visualize the amount-varying threshold and the decision regions
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ag = np.logspace(0, np.log10(amt_te.max()+1), 200)
ax[0].plot(ag, tau_star(ag), color="#2b6cb0", lw=2, label=r"$\tau^*$(amount)")
ax[0].axhline(0.5, color="#c53030", ls="--", label="flat 0.50"); ax[0].axhline(tau_g, color="#dd6b20", ls=":", label=f"global {tau_g:.3f}")
ax[0].set_xscale("log"); ax[0].set_yscale("log"); ax[0].set_xlabel("amount ($, log)")
ax[0].set_ylabel("decline if P(fraud) >"); ax[0].set_title("The bar to decline falls as the stake rises"); ax[0].legend(fontsize=8)

rng = np.random.default_rng(0)
leg = np.where(yte==0)[0]; samp = rng.choice(leg, 5000, replace=False); idx = np.concatenate([np.where(yte==1)[0], samp])
ax[1].scatter(amt_te[idx][yte[idx]==0], p_bayes[idx][yte[idx]==0], s=5, alpha=0.2, color="#a0aec0", label="legit")
ax[1].scatter(amt_te[idx][yte[idx]==1], p_bayes[idx][yte[idx]==1], s=14, color="#c53030", label="fraud")
ax[1].plot(ag, tau_star(ag), color="#2b6cb0", lw=2, label=r"$\tau^*$ — decline above")
ax[1].set_xscale("log"); ax[1].set_yscale("log"); ax[1].set_xlabel("amount ($, log)"); ax[1].set_ylabel("calibrated P(fraud)")
ax[1].set_title("Decision regions"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
c05=rcost(policies["flat 0.50"],yte,amt_te); cg=rcost(policies[f"global tau={tau_g:.3f}"],yte,amt_te)
print(f"cost-driven thresholding vs the 0.5 default: ${c05:,.0f} -> ${cg:,.0f} ({100*(c05-cg)/c05:.0f}% lower) on the same score.")
print("As in the ULB notebook, tau*(amount) is the provably expected-cost-optimal per-transaction rule; the")
print("amount here is real dollars, so its stake-dependent bar is fully authentic.")

cost-driven thresholding vs the 0.5 default: $444,385 -> $265,395 (40% lower) on the same score.
As in the ULB notebook, tau*(amount) is the provably expected-cost-optimal per-transaction rule; the
amount here is real dollars, so its stake-dependent bar is fully authentic.

C_REV = 3.0
tau_te = tau_star(amt_te)
auto_decline = p_lo > tau_te; auto_approve = p_hi < tau_te; review = ~(auto_decline | auto_approve)
for nm, m in [("auto-approve", auto_approve), ("REVIEW", review), ("auto-decline", auto_decline)]:
    print(f"  {nm:13s} {int(m.sum()):>7,} ({100*m.mean():5.2f}%)  containing {int(yte[m].sum()):>4} of {int(yte.sum())} frauds")
cost_triage = (float((auto_decline&(yte==0)).sum())*C_FP + float(((auto_approve&(yte==1))*CFN(amt_te)).sum())
               + float(review.sum())*C_REV)
print(f"\nreview queue {int(review.sum()):,} ({100*review.mean():.2f}%), surfacing {int(yte[review].sum())} frauds; "
      f"triage cost ${cost_triage:,.0f} (review @ ${C_REV:.0f}/case)")

  auto-approve  101,921 (88.22%)  containing 1280 of 3942 frauds
  REVIEW          3,314 ( 2.87%)  containing  235 of 3942 frauds
  auto-decline   10,299 ( 8.91%)  containing 2427 of 3942 frauds

review queue 3,314 (2.87%), surfacing 235 frauds; triage cost $198,030 (review @ $3/case)

# the dollar ladder, decomposed
flag_g = policies[f"global tau={tau_g:.3f}"]
def decomp(appr_fraud, decl_legit, n_rev=0):
    return (float((appr_fraud*CFN(amt_te)).sum()), float(decl_legit.sum())*C_FP, float(n_rev)*C_REV)
ladder = {"approve all\n(no model)": decomp(yte==1, np.zeros(len(yte),bool)),
          "model +\nflat 0.5":      decomp((policies['flat 0.50']==0)&(yte==1), (policies['flat 0.50']==1)&(yte==0)),
          "model +\ncost threshold":decomp((flag_g==0)&(yte==1), (flag_g==1)&(yte==0)),
          "+ uncertainty\ntriage":  decomp(auto_approve&(yte==1), auto_decline&(yte==0), review.sum())}
tot = {k: sum(v) for k,v in ladder.items()}
labels=list(ladder); miss=[ladder[k][0] for k in labels]; fdc=[ladder[k][1] for k in labels]; rev=[ladder[k][2] for k in labels]
fig, ax = plt.subplots(figsize=(10, 4.6))
ax.bar(labels, miss, color="#c53030", label="missed-fraud loss")
ax.bar(labels, fdc, bottom=miss, color="#dd6b20", label="false-decline friction")
ax.bar(labels, rev, bottom=np.array(miss)+np.array(fdc), color="#2b6cb0", label="analyst review")
for i,k in enumerate(labels): ax.text(i, tot[k], f"${tot[k]:,.0f}", ha="center", va="bottom", fontsize=9, fontweight="bold")
ax.set_ylabel("realized test cost ($)"); ax.set_title("Dollar ladder on the engineered IEEE-CIS score (out-of-time)"); ax.legend(fontsize=9)
plt.tight_layout(); plt.show()
base=tot["approve all\n(no model)"]; final=tot["+ uncertainty\ntriage"]; naive=tot["model +\nflat 0.5"]
print("realized out-of-time test cost by system:")
for k in labels: print(f"  {k.replace(chr(10),' '):26s} ${tot[k]:>10,.0f}")
print(f"\nno-model -> full system: ${base:,.0f} -> ${final:,.0f} ({100*(base-final)/base:.0f}% lower); "
      f"naive-0.5 -> full: {100*(naive-final)/naive:.0f}% lower, from the decision layer on top of the engineered score.")

realized out-of-time test cost by system:
  approve all (no model)     $   626,072
  model + flat 0.5           $   444,385
  model + cost threshold     $   265,394
  + uncertainty triage       $   198,030

no-model -> full system: $626,072 -> $198,030 (68% lower); naive-0.5 -> full: 55% lower, from the decision layer on top of the engineered score.

§	Content
1	The data: transaction ⋈ identity, imbalance, missingness, the 182-day time-split
2	D-column normalization — relative timedeltas → a fixed calendar anchor
3	UID construction — stitch transactions into card/account entities
4	Behavioral aggregates — per-entity velocity / recency / history
5	The lift, measured — GBM PR-AUC raw vs engineered, out-of-time, with SHAP
6	Completing the cycle — calibrate the score on real segments
7	Cost-sensitive decision — amount-varying threshold on real amounts
8	Uncertainty triage + dollar ladder — the realized-$ scoreboard
9	The two-notebook system, and what it demonstrates

Group	Fields	Meaning
Target / timing	`isFraud`, `TransactionDT`, `TransactionAmt`	label; seconds since a reference; amount
Product	`ProductCD`	product category (W/C/R/H/S) — a real segment
Entity keys	`card1…card6`, `addr1`, `addr2`	card & address attributes — the raw material for a UID (§3)
Timedeltas	`D1…D15`	days since prior events (e.g. `D1` ≈ days since the card was first seen) — relative, the target of §2
Counts	`C1…C14`	anonymized counting features (per entity/time)
Match flags	`M1…M9`	whether fields match (e.g. name-on-card vs address)
Identity	`id_01…id_38`, `DeviceType`, `DeviceInfo`	device / network attributes, present for a minority of txns

IEEE-CIS Fraud Detection — the feature-engineering front end¶

1 · The data — real structure, real mess¶

Data dictionary (the groups that matter here)¶

What §1 establishes¶

2 · The D-columns — from relative timedelta to fixed anchor¶

What §2 establishes¶

3 · Building the UID — heuristic entity resolution¶

What §3 establishes¶

4 · Behavioral aggregates — velocity, recency, history¶

What §4 establishes¶

5 · The lift, measured — raw vs engineered, with SHAP¶

What §5 establishes¶

6 · Completing the cycle — calibrating the engineered score on real segments¶

What §6 establishes¶

7 · Cost-sensitive decision — on real amounts¶

What §7 establishes¶

8 · Uncertainty triage and the dollar ladder¶

What §8 establishes¶

9 · The two-notebook system, and what it demonstrates¶

What the pair demonstrates¶

Honest limitations & next steps¶